Group Abstract Group Abstract

Message Boards Message Boards

0
|
5.2K Views
|
5 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Integrate a function in a 4D region with singularities using NIntegrate?

Miguel Condesso
Posted 7 years ago

Hi everyone!

I need to integrate a function in a 4D region (x1,y1,x2,y2), which explodes whenever x1=x2&&y1=y2. The code is as following:

    {Lx, Ly} = {1/50, 1/50};

    f1[x_, y_, n_, m_] := Cos[n*Pi*(Lx + x)/(2*Lx)]*Sin[m*Pi*(Ly + y)/(2*Ly)];
    f2[x_, y_, n_, m_] := Sin[n*Pi*(Lx + x)/(2*Lx)]*Cos[m*Pi*(Ly + y)/(2*Ly)]*m*Lx/(n*Ly);
    dist[x1_, y1_, x2_, y2_] := Sqrt[(x1 - x2)^2 + (y1 - y2)^2];

    int[f_, n_, m_, i_, j_] := int[f, n, m, i, j] =
    NIntegrate[f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2],
    {y1, -Lx, Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx}]

Then, I want to evaluate int using either f1 or f2, for positive integer (n,m,i,j), for instance using AbsoluteTiming[int[f1,1,2,1,2]]. I may need to compute these integrals thousands of times, for the different possible combinations of (n,m,i,j) - in the most complex cases, each of them may reach up to 14, for instance. So I must be able to calculate the integral as fast as possible.

However, the method is simply not working so well, for some. When I compute, for instance: int[f1,4,1,1,1], the computation simply lasts for ever and I receive some scary error messages:

NIntegrate::errprec: Catastrophic loss of precision in the global error estimate due to insufficient WorkingPrecision or divergent integral.

I have tried multiple integration strategies and methods, but the thing that seems to work best is to simply define the integration boundaries using {y1,-Lx,Lx},{x1,-Lx,Lx},{y2,-Ly,y1,Ly},{x2,-Lx,x1,Lx}, taking advantage of the integration order to define the singularity locations.

The integration method "DuffyCoordinates" performed quite fast, but I didn't know whether I was correctly using the "Corners" option, e.g.:

 int2[f_,n_,m_,i_,j_] := NIntegrate[f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2],
    {y1, -Lx, Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx},Method -> {"DuffyCoordinates","Corners" -> {0, 0, 1, 1}}]

... But it also gets stuck for (n,m,i,j)=(4,1,1,1). So I don't know what else to do:

  • Am I doing something wrong with this couple of integration strategies, or missing something obvious?
  • I suspect much of the problems are due to the integral actually evaluating to 0, maybe because the functions are somehow odd (whatever that means in 4D?), and I'm integrating them in a symmetric interval. But I don't know how I could define my function int in a way that it would check whether the integrand is odd in this 4D symmetrical integration volume.
POSTED BY: Miguel Condesso
5 Replies
Sort By:

Miguel,

You definitely now have a singularity now!

The problem with your integrand for 

int[f1,1,4,1,1]

Is that it has a discontinuity but the limit is zero as it approaches from either direction. One way to deal with this is to use an integration method that does not use regular intervals I tried 

int[f_, n_, m_, i_, j_] := 
 int[f, n, m, i, j] = 
  NIntegrate[
   f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2], {y1, -Lx, 
    Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx}, 
   AccuracyGoal -> 7, Method -> "AdaptiveMonteCarlo"]

and it worked fine (and fast). I could not go to any more accuracy than 7 without having numerical issues. I hope this helps.

Regards,

Neil
POSTED BY: Neil Singer

Miguel,

There is no singularity in your equations -- you may have a bug --you divide by dist[], however, dist is defined to be 1/Sqrt[] so the sqrt is in the numerator and not the denominator so the timing issues relate to the integrand approaching 0 and not a singularity.

Is your expression correct?

Regards,

Neil
POSTED BY: Neil Singer

Oh, this is unfortunate, I've made a mistake defining dist[] as you've noticed. The sqrt should be in the denominator of the integrand. Thank you, I've corrected the definition of dist[] in the original post.

Now the calculations take generally longer and I guess this is finally due to the singularity, which now actually exists.
POSTED BY: Miguel Condesso

Miguel,

If you add AccuracyGoal to the NIntegrate it will fix your issues.

I ran

int[f_, n_, m_, i_, j_] := 
 int[f, n, m, i, j] = 
  NIntegrate[
   f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2], {y1, -Lx, 
    Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx}, 
   AccuracyGoal -> 8]

Without a problem:

In[10]:= AbsoluteTiming[int[f1, 1, 2, 1, 2]]

Out[10]= {0.006115, -1.9303*10^-9}

You were asking for more accuracy than practical because the integral goes to zero and is dominated by lack of precision. I chose AccuracyGoal->8 but I was able to go up to 12 (with a time penalty). If you need more precision than that you will need to specify the AccuracyGoal and the WorkingPrecision.

Regards,

Neil
POSTED BY: Neil Singer

Hello Neil,

Thank you for the advice, that option alone speeds up the computation in a lot of cases. That does make a lot of sense, and actually I don't need a very high precision, 3-4 significant digits should suffice. But for instance, when computing int[f1,1,4,1,1], the integration takes a long time, if it finishes at all, even setting the AccuracyGoal->8. When I specify the volume of integration as {y1,-Lx,Lx},{x1,-Lx,Lx},{y2,-Ly,y1,Ly},{x2,-Lx,x1,Lx} in order to exclude the singularity at x1==x2&&y1==y2, it takes around 1.5 minutes in my computer, but this slows down the calculation in other cases, for instance for int[f1,1,2,1,3], and actually Mathematica does not complain about the singularity.

I was thinking about the functions that Mathematica uses to preprocess the integrand. Since it is the product of sines and cosines, divided by a square root, it is probably easy to identify if it is odd/even (but I'm not sure how this works in a 4D space), and since the integration region is symmetric, maybe the calculations could be simplified this way.
POSTED BY: Miguel Condesso
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard