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

Inconsistent values when integrating [Python]

  1. Jun 15, 2015 #1
    I have a 2D Gaussian:

    ## f(x,y) = e^{-[(x-x_o)^2 + (y-y_o)^2]/(2*{sigma}^2)}##

    which I converted into polar coordinates and got:

    ## g(r,θ) = e^{-[r^2 + r_o^2 - 2*r*r_o(cos(θ)cos(θ_o) + sin(θ)sin(θ_o))]/({2*{sigma}^2})} ##

    The proof for how this was done is in the attached file, and it would be great if someone could verify my steps in case I messed up somewhere (although I have inputted values and it seems to work). I have plotted the functions, and they do seem comparable by visual inspection, too.

    The characteristics of the functions (including ## sigma_o##) are manually specified by me, so the values for ##x_o##, ##y_o## match up with ##r_o## and ##θ_o##. Now, when I integrate these functions over the same regions, I don't always get the same answers. For example, when ##r = 0## or ##θ = 0##, I do get the same answers, but when ##θ \neq 0##, then the answers are slightly off. I have been trying to search for where I'm going wrong, and it may be completely obvious (likely associated with the sine and cosine terms), but I'm just not seeing it. The fact that the values from completing the integration on Python (in polar coordinates) and WolframAlpha (in cartesian coordinates) are relatively similar for the cases where ##θ \neq 0## seems a little odd to me. If anyone has any thoughts, it would be greatly appreciated!

    Here is my code:

    Code (Text):

    import numpy as np
    from scipy import integrate

    RIT = [] # Region Intensity (integration)
    i = 0
    N = 1
    while i < N: #1 random beam generated
        i = i + 1
        sigma0 = 5 # just an example--arbitrary, I usually set it between 1 and 10
        r0 = #you can input any value here, I usually set it between 0 and 10
        theta0 = random.uniform(0,np.pi*2) #you can make it arbitrary instead of random if you wish
        def G(r,theta): #this is the Gaussian in polar coordinates
            return (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 includes dr and dtheta, NOT dx and dy any longer!!!
        RI = integrate.nquad(G, [[4,7],[0,0.5*np.pi]]) #this region is between r = 4 and 7 in the first quadrant
        RIT.append(RI)

    print RIT

     
    I've also been trying points in both f and g to see if they correspond, and they seem to work. For example, if:

    ##sigma = 1##
    ##x_o = 1##
    ##y_o = 2##

    then ## f(1,1) = e^{-0.5} ##

    This also corresponds to:

    ##r_o = \sqrt{26} ##
    ##θ_o = ~1.107##
    then ##g(\sqrt{2},pi/4) = e^{-0.5}##

    I am leaning towards there being an error in my actual integration, but do you see any?

    If I use these values for the code:
    ##sigma0 = 5##
    ##r0 = 2**0.5##
    ##theta0 = 1.75*np.pi##

    Then the value for the integral I get is: 1.33556147e+01
    When I do this in WolframAlpha, I get: 12.7035

    Any advice is welcome!
     

    Attached Files:

    Last edited: Jun 15, 2015
  2. jcsd
  3. Jun 15, 2015 #2
    Here is the WolframAlpha calculation in cartesian coordinates that should correspond to my polar coordinate integration in Python.
     

    Attached Files:

  4. Jun 15, 2015 #3

    RUber

    User Avatar
    Homework Helper

    I am not understanding why you have an imaginary part of your cartesian integral. You probably need to break that integral into two parts.
    for x = 0 to 4 and for x=4 to 7.
     
  5. Jun 15, 2015 #4
    You're right. Let me check that out. Complete oversight by me haha. Thank you!
     
  6. Jun 15, 2015 #5

    DEvens

    User Avatar
    Education Advisor
    Gold Member

    It's a little difficult to understand what exactly you are integrating. It seems like you are integrating from ##r## 4 to 7, and ##\theta## from 0 to pi/2. So this is some chunk of a 2-D Gaussian. Without checking if you have done that correctly, you are worried that when you change ##r_0## and ##\theta_0## that you get slightly different values.

    Well, yes. You are integrating a different part of the Gaussian when you do that. Or so it seems.

    Also, I want to pile on what RUber asked. What is the integral you fed to Wolfram? What's with the square roots in the integrand?
     
  7. Jun 15, 2015 #6
    And it turns out my code is correct, I just messed up my integral! Sweet! Thank you for catching my mistake!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Inconsistent values when integrating [Python]
  1. Strings in python (Replies: 15)

Loading...