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

Posted 5 months ago
650 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[HolonomicDifferentialRootReduceDumpy[NIntegrateLevinRuleDumpx]+(HolonomicDifferentialRootReduceDumpy^\[Prime]\[Prime])[NIntegrateLevinRuleDumpx],NIntegrateLevinRuleDumpFastLookupHolonomicDifferentialEquation] 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
Sort By:
Posted 5 months ago
 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 5 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.0005214754242472478922017323553271730151770785829665971637152618009066725650264883580. and 0.001751501250207990399788234152530489439832590360506388875897758779149405931549588380. 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 
 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?