Results 1 to 11 of 11

Math Help - Numerical integration over a hemisphere

  1. #1
    Newbie
    Joined
    Sep 2008
    Posts
    8

    Numerical integration over a hemisphere

    Dear all,

    I am trying to solve the following problem:
    A function is given on the surface of a unity hemisphere, and I need to find an efficient way to numerically integrate it over the surface.

    I started with an equidistant approach (both angles of the spherical coordinates subdivided into constant intervals). This gives good results but needs a large amount of computations.
    In the internet, I found an alternative (Gaussian) approach with equidistant points, but only for integration over spheres, not over hemispheres. This approach can be found here: Integration nodes for the sphere

    Can somebody help me on doing the numerical integration over the hemisphere? Any help is appreciated!

    Thank you,
    Jens
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Not sure what you want but if you choose to use Mathematica, here's the code I'd use to numerically integrate the function f(x,y,z) over the upper hemisphere of the unit sphere.

    Code:
    NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
       Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1}, 
      {y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]
    I realize not everyone has Mathematica or even wants to use it. Just sayin' if you're interested in a brute force approach to the answer using Mathematica, which I'd think would be rather quick (at 2GHz) for a reasonable function, that's what I'd start with. Here's a check of f(x,y,z)=x^2 over the upper hemisphere of the unit sphere, the answer of which is 2\pi/3, just to get some indication if the coding is correct:

    Code:
    f[x_, y_, z_] := x^2; 
    NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
       Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1}, 
      {y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]
    N[2*(Pi/3)]
    
    2.0944
    
    2.0944
    That took lest than 1/10 of a second. Can you post your function?
    Last edited by shawsend; September 18th 2008 at 09:01 AM. Reason: added example for check
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Sep 2008
    Posts
    8
    shawsend,

    Thank you for your answer!
    I actually can't use Matematica in this case, because I need to include the integration (even for one special function only) into my own code. I'd like to find some algorithm and then implement it into my program for this one special case only.

    Well, here is the function:
    you are given 4 angles (describing 2 points on a unit hemisphere):
    \gamma, \varphi \in [-\pi, \pi]
    \beta, \theta \in [0,\pi]
    \gamma and \theta are considered variables, while \varphi and \beta are given parameters.

    The function I would like to integrate is the cosine of the angle between the position vectors of the 2 points:
    f_1(\gamma, \theta) = \cos\theta \sin\beta \cos(\gamma-\varphi) + \sin\theta \cos\beta
    Let f_2(\gamma, \theta) = max[0,f_1(\gamma, \theta)], meaning that all angles >\pi/2 are discarded.

    \theta is the solar altitude angle. Since it is mathematically more common to start counting the angle from the z-axis, not from the x-y-plane, the following transformation is performed: \theta_z = \pi/2-\theta. This transforms the solar altitude angle \theta into the zenith angle \theta_z. All following calculations are based on the zenith angle. In terms of \theta_z, the functions f_1 and f_2 read:
    f_1(\gamma,\theta_z) = \sin\theta_z \sin\beta \cos(\gamma-\varphi) + \cos\theta_z \cos\beta
    f_2(\gamma, \theta_z) = max[0,f_1(\gamma, \theta_z)]

    Then, the integration I am trying to calculate numerically is:
    \phi = g(\beta,\varphi) = \frac{1}{\pi}\int_{\gamma=-\pi}^{\pi}{\int_{\theta_z=0}^{\pi/2}{f_2(\gamma, \theta_z)}\sin\theta_z\ d\theta_z d\gamma}
    I hope that you can now better see what I mean and what I want to do.
    Thanks in advance for your help!
    Jens
    Last edited by cool.bambus; September 24th 2008 at 01:40 PM. Reason: formula for f_1 corrected + simplified; another error corrected...
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Ok, I see what you mean. I've done numerical integration in C++. Just easier with Mathematica. Also, can link Mathematica to an external program but that's a programming project and something you're probably not interested in. I may still work with your integral to see what I get with Mathematica.
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Newbie
    Joined
    Sep 2008
    Posts
    8
    thanks for your effort!
    i am looking forward to getting your results.

    Jens
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Bamus, your notation is unusual to me (not the standard way to set up spherical coordinates and dot product of two vectors, at least to me). In your notation, the point p(\gamma,\beta), \gamma is measured in the counter-clockwise direction from the negative x-axis, and \beta is measured up from the x-y plane? Same with the other two symbols?

    Also, your integral looks a bit odd as your integrating f_2(\gamma,\theta) but using d\theta d\phi. I would gusss perhaps it would be something like:

    g(\beta,\varphi)=\frac{1}{\pi}\int_{-\pi}^{\pi}\int_{0}^{\pi} f_2(\gamma,\theta)sin(\theta)d\theta d\gamma

    Maybe no though. Anyway, generating a plot for g(\beta,\varphi) is problematic as the integration is taking a long time to complete. Probably due to the complicated integrand. So I switched to Monte Carlo method. This was better. Perhaps you might want to try this method in your approach to the problem
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Newbie
    Joined
    Sep 2008
    Posts
    8
    shawsend,

    thanks again.
    Sorry for the mistakes in the formulae and thank you for pointing them out. I corrected the errors, see my edited post above.

    Actually, I do not need to have the value of \phi, the integral, calculated as a general function \phi(\beta, \varphi). \beta and \varphi are just parameters, not variables in the proper sense. I am given values for \beta and \varphi at the start of the calculation, and I am supposed to do the calculation of \phi only once and only for this pair of \beta and \varphi.

    So, for the purpose of calculating the integral, we can treat \beta and \varphi as if they were real constants. Will we still need Monte Carlo? I have never worked with that method...

    Jens
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Hello Bamus,

    How about posting your values of \beta,\varphi. I'd like to generate a more accurate integral and see if it agrees with your value. Also Monte Carlo is a random method. Just like if you wanted to calculate the area under a parabola in the unit square, generate 100,000 random points in that square and count the number inside the parabola. The ratio inside/outside gives a pretty accurate value to the integral. The method is commonly used to calculate many-iterated integrals of many variables.

    That's just one method out of many available to Mathematica, but I would think that's unnecessary since we're only calculating one point and not 10,000 or so to make a decent plot of g(\beta,\varphi).
    Follow Math Help Forum on Facebook and Google+

  9. #9
    Newbie
    Joined
    Sep 2008
    Posts
    8
    -> function f_1(\gamma, \theta) corrected and simplified a little bit, see post above.

    Quote Originally Posted by shawsend View Post
    How about posting your values of \beta,\varphi. I'd like to generate a more accurate integral and see if it agrees with your value.
    For instance, typical values for \beta and \varphi would be:
    \beta = 35 = 0.61086524 rad
    \varphi = 15 = 0.26179939 rad

    You could check with
    \beta = 0 and \varphi = 0 which must give g(0,0) = 1
    \beta = 90 and \varphi = 0 which must give g(\pi/2,0) = 1/2


    Quote Originally Posted by shawsend View Post
    That's just one method out of many available to Mathematica, but I would think that's unnecessary since we're only calculating one point and not 10,000 or so to make a decent plot of g(\beta,\varphi).
    That's why in a recent post I pointed out that I don't need a general function g(\beta,\varphi).

    Thanks for your help!
    Jens
    Last edited by cool.bambus; September 21st 2008 at 07:48 PM.
    Follow Math Help Forum on Facebook and Google+

  10. #10
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Then I must be doing something wrong because I don't get those answers for g(0,0] and g(\pi/2,0). The code below is a little messy because of the greek letters but you can clearly see the last two command are the results for these two values.

    f[\[Gamma]_, \[Beta]_, \[CurlyPhi]_, \[Theta]_] :=
    Max[0, Cos[\[Theta]]*Cos[\[Gamma] - \[CurlyPhi]]*Sin[\[Beta]] +
    Cos[\[Beta]]*Sin[\[Theta]]]
    g[\[Beta]_, \[CurlyPhi]_] := (1/Pi)*NIntegrate[f[\[Gamma], \[Beta], \[CurlyPhi], \[Theta]]*
    Sin[\[Theta]], {\[Theta], 0, Pi}, {\[Gamma], -Pi, Pi},
    Method -> MonteCarlo];
    g[0, 0]
    g[Pi/2, 0]

    Out[3]=
    3.1486450543682474

    Out[4]=
    0.6402232798879435

    However I did generate a 15x15 array of points for g(\beta,\varphi) and plotted them. Is the plot below consistent with your expectations?
    Attached Thumbnails Attached Thumbnails Numerical integration over a hemisphere-bamusplot.jpg  
    Follow Math Help Forum on Facebook and Google+

  11. #11
    Newbie
    Joined
    Sep 2008
    Posts
    8
    I apologize again... I have been confusing about the use of the solar altitude angle \theta vs. the zenith angle \theta_z. This changes some \sin into \cos functions (and vice versa) in the functions. See the details and the corrected function in my post above, containing the definition of f_1(\gamma, \theta_z).

    The two tests yield:
    Special case 1:
    \beta = 0, \varphi = 0:
    In this case, f_1(\gamma,\theta_z)=\cos \theta_z, and as f_1\ge0 is always true for \theta_z \in [0,\pi/2], the value of the integral is: g(0,0) = 1 (checked this by hand).

    Special case 2:
    \beta = 90, \varphi = 0:
    In this case, f_1(\gamma,\theta_z)=\sin \theta_z \cos\gamma. Here, the condition f_1(\gamma,\theta_z)\ge0 is equivalent to \gamma \in [-\pi/2,\pi/2]. Thus, the value of the integral g(\pi/2,0) = 1/2 can be easily calculated by hand.
    Sorry again for the confusion.
    Now the problem should be well-defined...

    Jens
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. numerical integration
    Posted in the Calculus Forum
    Replies: 1
    Last Post: April 9th 2010, 02:49 AM
  2. Example of Numerical Integration?
    Posted in the Calculus Forum
    Replies: 1
    Last Post: November 7th 2009, 04:22 PM
  3. Numerical integration
    Posted in the Advanced Applied Math Forum
    Replies: 1
    Last Post: November 6th 2009, 09:25 AM
  4. multiple integration: hemisphere/cylinder
    Posted in the Calculus Forum
    Replies: 5
    Last Post: November 5th 2009, 03:39 PM
  5. Numerical Integration
    Posted in the Calculus Forum
    Replies: 2
    Last Post: August 30th 2007, 06:30 AM

Search Tags


/mathhelpforum @mathhelpforum