Numerical integration over a Green's function

[member 428835]

Hi PF!

I'm numerically integrating over a Green's function along with a few very odd functions. What I have looks like this

NIntegrate[-(1/((-1.` + x)^2 (1.` + x)^2 (1.` + y)^2))
     3.9787262092516675`*^14 (3.9999999999999907` +
     x (-14.99999999999903` +
        x (20.00000000000097` - 9.999999999999515` x +
           1.` x^3))) (-1.` + y)^2 (4.` + y) BesselJ[4,
    304.6888201459785` Sqrt[1 - x^2]] BesselJ[4,
    304.6888201459785` Sqrt[1 - y^2]] Cosh[
    310.00637327206255` - 304.6888201459785` x] Cosh[
    310.00637327206255` - 304.6888201459785` y] /. {x -> xx,
   y -> yy}, {xx, yy} \[Element]
  ImplicitRegion[
   Cos[Cos[\[Alpha]] Cos[\[Pi]/180]] < xx < yy < 1, {xx, yy}]]

but Mathematica throws a "numerical integration converging too slowly" error. How would you treat this?

I can send you the full notebook if you're interested, though that would have to be private. Thanks so much!

EDIT: For completeness the following technique worked in Mathematica and agrees with MCI method outlined in python below, though Mathematica is much much faster: redefine $G$ as a piecewise function, since after all it is (my technique above was to split the integration of $G$ over two separate domaines, one where $x>y$ and one where $y>x$). Then run NIntegrate over the full square domain. However, under Method specify LocalAdaptive. After much research I decided to use this technique, and no more errors or holdups, all is well.

 

[jambaugh]

If I were you I would code my own numerical integration. You could for example attempt Monte Carlo Integration. This involves generating random points uniformly distributed over a square region (with area A) containing your domain of integration.
For each point, test if it's in the region of integration and then evaluate your integrand function, then add that value to your cumulative sum (and increment your counter N for every randomly generated case, not just the ones in your domain of integration!).

 

[member 428835]

If I were you I would code my own numerical integration. You could for example attempt Monte Carlo Integration. This involves generating random points uniformly distributed over a square region (with area A) containing your domain of integration.
For each point, test if it's in the region of integration and then evaluate your integrand function, then add that value to your cumulative sum (and increment your counter N for every randomly generated case, not just the ones in your domain of integration!).

Can you explain the advantage to this approach over what NIntegrate currently does? Also, notice the region with which I am integrating: it's essentially $$\int_{\cos(\cos \pi/180)\sin(\pi/180)}^1\int_{\cos(\cos \pi/180)\sin(\pi/180)}^y f(x,y)\, dxdy$$ which I don't know how to perform numerically since a variable is in the integration limit.

 

[member 428835]

Also, the integrand when plotted over the integration region $[\cos(\cos \pi/180)\sin(\pi/180),1]\times [\cos(\cos \pi/180)\sin(\pi/180),1]$ looks very well behaved as you can see:

 

[jambaugh]

I was thinking about how to implement so, I spent about 30 minutes and wrote a quick example python script. Code included below. Of course, I used a much simpler integration problem but I think you can adapt the code easily enough. I later thought of an adaptation that would calculate the standard deviation of the estimator (sum*Area/N). I might fiddle around with it and pretty it up and will upload it here if I do.

(Basically, your sum*Area/N = E is your estimator. You can likewise estimate its standard deviation and construct a confidence interval of whatever confidence level and precision by letting it run until that is reached. I'm sure there are plenty of elegant packages out there which would handle it neatly.)

I can't speak to the advantages over NIntegrate because I simply don't know what it does. I can't speak to why Wolfram Alpha gave you the message it did.

As to the advantages or disadvantages of MCI, one advantage is that it is simple to implement, and it can reach lower precision results rather quickly. It also will handle higher dimensional integrals without much extra coding or resources. (I believe it is quite efficient memory-wise.) As with all mean type estimators, its standard error will decrease in proportion $\frac{1}{\sqrt{N}}$. I am not an expert and defer to those who are or what you find in searches on the topic. It also, due to the use of random sampling points, will avoid certain biases which might occur due to the regularity of the traditional meshes of direct numerical integration. Contrawise, if there are very localized features it could easily miss them until it's taken a very dense sample.

One point, the smaller your region of integration is relative to the rectangular envelope, the less efficient since you're region test is rejecting so many cases. Another point is that you can adapt it to use non-uniform sampling (offset with bias in the weights of the samples inversely proportional to their probability density.) This will allow you to, for example, do a MCI over an unbounded region or increase the density near pathological features. But I don't know how easy the error analysis would be.

Sorry if I'm rambling. Your question got me thinking about this all day. I think I might offer a project like this to students in my department for an undergraduate capstone.

Python Code Below: (It went quickly even with this slow interpreter. A c language implementation would fly.)
[Edit:] Note that the # symbol begins a comment on that line in Python.
from math import * from random import * # Monte Carlo method, integrating z>0 hemisphere of radius 2 over the unit circle in the x,y plane. def f(pt): # The integrand function, here z = sqrt( 4-x^2 - y^2) [python uses ** for powers like fortran not ^ like c and most other lang.] x= pt[0] y= pt[1] return (4.0 - x**2-y**2)**0.5 def testRegion(pt): return (pt[0]**2 + pt[1]**2 < 1.0) def genpoint(): # generating coordinates uniformly in the square -1<x<1, -1<y<1 (area = 4). x = 2.0*random()-1.0 y = 2.0*random()-1.0 return (x,y) #Intitialize Sum = 0.0 Area = 4.0 N = 0 #Main Loop. while(True): for i in range(10000): # Number of iterations before printing current result. N+=1 pt = genpoint() if testRegion(pt): Sum += f(pt) print(Sum*Area/N)

 

[pbuk]

Python code works much better posted like this:

from math import *
from random import *

#  Monte Carlo method, integrating z>0 hemisphere of radius 2 over the unit circle in the x,y plane.

def f(pt):
    # The integrand function, here z = sqrt( 4-x^2 - y^2)   [python uses ** for powers like fortran not ^ like c and most other lang.]
    x= pt[0]
    y= pt[1]
    return (4.0 - x**2-y**2)**0.5

def testRegion(pt):
    return (pt[0]**2 + pt[1]**2 < 1.0)

def genpoint():
    # generating coordinates uniformly in the square -1<x<1, -1<y<1 (area = 4).
    x = 2.0*random()-1.0
    y = 2.0*random()-1.0
    return (x,y)
 
#Intitialize

Sum = 0.0
Area = 4.0
N = 0

#Main Loop.

while(True):
    for i in range(10000):   # Number of iterations before printing current result.
        N+=1
        pt = genpoint()
        if testRegion(pt):
            Sum += f(pt)
    print(Sum*Area/N)

But in any case if you want to try MC then you can do this with NIntegrate simply by supplying Method -> "MonteCarlo" as an argument. You could also try some more sophisticated sampling-based algorithms.

Don't be surprised if this does not solve the convergence problem though, depending on what tolerance you are looking for.

 

[pbuk]

What are all those constants like 3.9999999999999907 for? This looks like it should be 4 but includes an error term due to some previous step using a numerical method.

Edit: and where has $\cos \left(\cos \left(\frac \pi {180}\right)\right)$ come from? Whilst it is possible to calculate this number it cannot have any physical or even mathematical significance.

 

[member 428835]

Wow, thank you both for the replies! Never used MCI before, but I'm checking it out. Very cool idea, and nice clear python script you wrote.

pbuk: I copy-pasted Mathematica output without really reading it. I didn't want to mess anything up by rounding so I thought this was easier. Those numbers stem from previous numerical calculations. And the cosine terms you ask about are a typo on my part: the integration limits in post 3 should read from $[\sin\alpha,1]\times[\sin\alpha,1]$

 

[member 428835]

Also, regarding Monte Carlo integration, is it preferred to compute 100 integrals each with 100 samples, or 1 integral with 10,000 samples? Seems like it should be the same, right?

 

[jambaugh]

Python code works much better posted like this: [...]

Thanks! I should have spent the time to look this up. The key-word highlighting makes it so much better. The error message mentioned by the OP makes me wonder if it actually was trying to use MCI. It would explain the "converging too slowly" type message since the algorithm, if MC would only detect the rate of convergence empirically given this method.

The question then is can you, using Nintegrate[] , change the loop limits, i.e. tell it to keep going further in trying to detect convergence? (As you might infer, my reply to OP was tempered by my ignorance of Mathematica's internals.

@joshmccraney I was "cludging" it together and chose my inner loop size to give appreciable convergence between the printed "spot checks". I added zeros until it fit on one page (giving about 3-4 decimals precision on the example problem within a single page's worth of output.)

 

[Orodruin]

Edit: and where has $\cos \left(\cos \left(\frac \pi {180}\right)\right)$ come from? Whilst it is possible to calculate this number it cannot have any physical or even mathematical significance.

Ok, three months old but I stumbled upon it just now.

I have come across this type of expression with relevant meaning exactly once that I can remember. I think it involved the components of a vector being parallel transported. The angle that the vector made to some reference set of curves varied as a cosine with the curve parameter and therefore the components took the functional form of nested trigonometric functions. I am not saying it is common but expressions like this do turn up sometimes.

 

[pbuk]

The cosine of the cosine of one degree?

 

[Orodruin]

Well, the argument would depend on the exact curve. I am sure it could be made one degree.

 

[member 428835]

Edit: and where has $\cos \left(\cos \left(\frac \pi {180}\right)\right)$ come from? Whilst it is possible to calculate this number it cannot have any physical or even mathematical significance.

Tough to explain briefly, but it's the edge of a curve that is ultimately revolved around an axis. In fluid dynamics this is called the contact-line, in this case of an equilibrium spherical-capped interface for liquid in a cylinder.