Message Boards Message Boards

GROUPS:

Integrate a function in a region around the origin with a singularity?

Posted 3 months ago
436 Views
|
4 Replies
|
1 Total Likes
|

I need to integrate a function in a region around the origin, the problem is that it has a singularity there. Having tried Integrate[] without any sort of response after waiting for a long time, I tried to perform the integration numerically. The code is something like this:

{Lx,Ly}={1/50, 1/50};
Aelem[x_, y_, n_, m_] := Aelem[x, y, n, m] = 
NIntegrate[current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2], {y1, -Ly, Ly}, {x1, -Lx, Lx}, Exclusions -> {0, 0}, Method -> {"LocalAdaptive"}]

Helem[i_, j_, n_, m_] := Helem[i, j, n, m] =
NIntegrate[current[x, y, i, j]*Aelem[x, y, n, m], {y, -Ly, Ly}, {x, -Lx, Lx}, Method -> {"LocalAdaptive"}]

And then these current functions have the shape

g[x_,y_,n_,m_]:=(2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]
current[x_,y_,n_,m_]:=D[g[x, y, n, m], y]
Helem[1,1,1,1]//AbsoluteTiming

I'm trying out some different integration strategies, and different ways to define the current function, to try to get the best solution. Running Helem[1,1,1,1]//AbsoluteTiming has some problems: 1. Sometimes the evaluation just takes for ever 2. I get the " NIntegrate::inumr error message: "The integrand has evaluated to non-numerical values (...)"" message, so I use ?NumericQ, but then I have problems with defining the current as the derivative of g; 3. Sometimes I get a really weird error message:

Throw::nocatch: Uncaught Throw[Holonomic`DifferentialRootReduceDump`y[NIntegrate`LevinRuleDump`x]+(Holonomic`DifferentialRootReduceDump`y^\[Prime]\[Prime])[NIntegrate`LevinRuleDump`x],NIntegrate`LevinRuleDump`FastLookupHolonomicDifferentialEquation] returned to top level.

I really need to be able to perform the integral as fast as possible, I need the function H for many combinations of (i,j,n,m), and I feel like I am shooting in the dark about the integration strategies and all. Does anyone have any input on how to improve this and make the calculations faster? Any input is appreciated!

Edit (25/06/2018): I am now trying to integrate all at once, so something like:

{Lx,Ly}={1/50, 1/50};
gnm[x_, y_, n_, m_] := (2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]
current[x_, y_, n_, m_] = D[gnm[x, y, n, m], y]; 
atonce[i_, j_, n_, m_] := atonce[i, j, n, m] =
NIntegrate[current[x, y, i, j]*current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2]
, {y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Lx, Lx}, {x1, -Lx, Lx}]
atonce[1,2,1,2]//AbsoluteTiming

And this is a little insane, whenever I try a different integration strategy (GlobalAdaptive, DuffyCoordinates, DoubleExponential), increase the MaxRecursion or the WorkingPrecision, the result becomes very different. I've tried to put both these parameters as high as 300, but I keep getting the error message "NIntegrate failed to converge to prescribed accuracy after X recursive bissections", or equivalent. I am considering trying something like integrating first on dx1dy1, somehow defining "x" and "y" as singularities, and only then integrate on dxdy, but I am unsure how to do this.

Edit 26/06/2018: Corrected the gnm function

4 Replies

Try:

{Lx, Ly} = {1/50, 1/50};
gnm[x_, y_, n_, m_] := (2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/2]
current[x_, y_, n_, m_] = D[gnm[x, y, n, m], y];
atonce[i_, j_, n_, m_] := 
 NIntegrate[
  current[x, y, i, j]*
   current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2], {y, -Lx, 
   Lx}, {x, -Lx, Lx}, {y1, -Lx, Lx}, {x1, -Lx, Lx}, 
  Exclusions -> (Sqrt[(x - x1)^2 + (y - y1)^2] == 0)]
atonce[1, 2, 1, 1] // AbsoluteTiming

(* {0.0139544, 0.} *)
Posted 3 months ago

I'm very sorry had mistyped the gnm function, now it is corrected:

gnm[x_, y_, n_, m_] := (2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]

Adding the Exclusions option, the Integration becomes a little slower, but the result is the same and the error message is still there:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 36 recursive bisections in x near {y,x,y1,x1} = {-0.0115588,-0.0025,-0.0115588,-0.000779382}. NIntegrate obtained -4.16285*10^-6 and 0.000030333190157773638` for the integral and error estimates.

Then, I impose WorkingPrecision -> 30, MaxRecursion -> 100 and the message becomes

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.00052147542424724789220173235532717301517707858296659716371526180090667256502648835`80. and 0.0017515012502079903997882341525304894398325903605063888758977587791494059315495883`80. for the integral and error estimates.

After this, I increase the suggested (by the error message) MaxErrorIncreases option: Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 3000}, and the error message stays the same, just now complaining that the error increased 3000 times, and my NIntegrate command already looks like this:

atonce[i_, j_, n_, m_] :=
NIntegrate[current[x, y, i, j]*current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2],{y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Lx, Lx}, {x1, -Lx, Lx},
Exclusions -> (Sqrt[(x - x1)^2 + (y - y1)^2] == 0), Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 3000}, WorkingPrecision -> 30, MaxRecursion -> 100]
atonce1[1, 2, 1, 1] // AbsoluteTiming

Multi-dimensional integral they are always hard to integrate. The answer is not simple. In your case we have problems by singularity.

I do not have an answer what going on inside NIntegrate, because I do not have such an experience related to this function.

Good Luck.

Regards Mariusz I.

Useful websites:

site1,site2

Posted 3 months ago

Thanks for the suggestions anyway! I have finally arrived at something promising! Basically I changed the method of handling the singularities, and the way I declare the exclusions. Any further suggestions that may speed things up are more than welcome! Here's the code:

Clear[Lx, Ly, gnm, current, try1]
{Lx, Ly} = {1/50, 1/50};

gnm[x_, y_, n_, m_] := (2*Lx/(n*\[Pi]))*Sin[n*\[Pi]*(x + Lx)/(2*Lx)]*Sin[m*\[Pi]*(y + Ly)/(2*Ly)];
current[x_, y_, n_, m_] = D[gnm[x, y, n, m], y];

try[i_, j_, n_, m_] := NIntegrate[jxnm[x, y, i, j]*jxnm[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2],
{y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Ly, y, Ly}, {x1, -Lx, x, Lx},Method -> {"GlobalAdaptive", "SingularityHandler" -> "DuffyCoordinates"}];

try[1, 1, 1, 1] // AbsoluteTiming

Edit I think it's especially slow when the value tends to 0, how can one deal with this?

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract