Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python Python: Integrating a function but getting large errors

  1. May 20, 2015 #1
    My code:
    Code (Python):

    import numpy as np
    import matplotlib.pylab as plt
    import math
    import random
    from scipy import integrate

    R1 = .001
    R2 = 7

    def G(r,theta):
        sigma = random.randint(4000., 7000.)/1000. # width of beam is 4 - 7mm
        r0 = random.randint(0, R1*1000.)/1000. #random centroid
        theta0 = random.uniform(0, np.pi*2) #angle
        return (100/(np.pi*2*sigma**2))*(np.e**(-((r**2 + r0**2 - 2*r*r0*(np.cos(theta)*\
        np.cos(theta0) + np.sin(theta)*np.sin(theta0)))/2*sigma**2)))*r # this r is here because the integrand is dr and dtheta, NOT dx and dy any longer!!!
    #amplitude is 100

    RTI = integrate.nquad(G, [[0,R2],[0,2*np.pi]]) #Integration

    print RTI
    I keep running this code, and everything seems to be working, but I keep getting results like: (0.1406360958428128, 0.012357271987783056)

    Also, the code itself takes quite a while to run and I'm sure this is an indication of something...just not sure where it's going on. As you can see, the error is almost the value of integration itself, and this error is to large. I've been looking at the code but am not sure why such a large error is arising. If anyone has any ideas, it would be greatly appreciated!
    Last edited by a moderator: May 20, 2015
  2. jcsd
  3. May 20, 2015 #2


    User Avatar
    Science Advisor

    I'm not sure what exactly you are trying to do. The integrate.nquad function will call the function G many times as it evaluates the integral. Each time it calls G, the parameters sigma, r0 and theta0 will have different values. So the value of the integral isn't really defined, because the function G is a different function each time it is called. I suspect that this isn't what you intend. You should probably move the random.randint calls outside of the G function and call them once to define sigma, r0 and theta0 before evaluating the integral. Does this make sense?
  4. May 21, 2015 #3

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    It's an indication that the integrator is struggling.

    A lousy numerical integrator might arbitrarily chop the space to be integrated into fixed size chunks, call the derivative function at each grid point, and apply some mathematical formulae to generate an estimate of the integral over than space. A better integrator will subdivide the space according to sensitivities revealed by calls to the derivative function.

    Your derivative function involves random. A good integrator will detect that. What a good integrator won't do is detect that your derivative function (in the context of computer science) is not a "function" (in the context of mathematics). Any integrator must necessarily assume that the function to be integrated is indeed a function. You don't have a function. You have randomness.

    What, exactly, are you trying to integrate?
  5. May 21, 2015 #4
    I have attached an image of what I am trying to do. Yes, you seem to get what I was attempting. I want to generate a single particular gaussian with the variables r0, theta0 and sigma0 randomly generated but held constant for that one run, then I want to integrate that one particular function over the circle. I would be creating numerous random gaussian functions like this.

    I have modified the code a little bit based on your suggestions, and it seems to be working, but any input with regards to thingsto modify would be great!

    Code (Text):
    RTIT = []

    R1 = 4. #
    R2 = 7. #

    i = 0

        while i < 10:

        i = i + 1

        sigma0 = random.randint(4000., 7000.)/1000. # width is 4 - 7mm
        r0 = random.randint(0, R1*1000.)/1000. #random centroid
        theta0 = random.uniform(0, np.pi*2) #angle centroid makes with polar axis

        def G(r,theta): #this is the Gaussian
            return (100/(np.pi*2*sigma0**2))*(np.e**(-((r**2 + r0**2 \
            - 2*r*r0*(np.cos(theta)*np.cos(theta0) + \
            np.sin(theta)*np.sin(theta0)))/2*sigma0**2)))*r # this r is here because
    #the integrand is dr and dtheta, NOT dx and dy any longer!!!
    #intensity is 100

        RTI = integrate.nquad(G, [[0,R2],[0,2*np.pi]]) #Total

    print RTIT

    Attached Files:

    Last edited by a moderator: May 21, 2015
  6. May 21, 2015 #5


    User Avatar
    Science Advisor

    A couple of things:

    If you enclose your code in
    Code (Text):
    tags, then the code will be more readable.

    Your code with the modifications looks good. Did it fix the large error problem you were initially having?
  7. May 22, 2015 #6
    Yep, it seems to be working great now. Thanks! I'll definitely use that from now on.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook