2
$\begingroup$

I need to integrate a function with a singularity at the origin. I need this integration to happen quite fast, and while Integrate[] simply keeps on going forever, using NIntegrate with LocalAdaptive does the job (although not really fast). The thing is that I keep getting the NIntegrate::inumr error message: "The integrand has evaluated to non-numerical values (...)". The code is sort of like this:

Edit: So the current function within Aelem actually depends on (x1,y1) and should therefore also be integrated in the first step

    Aelem[x_, y_, current_, n_, m_] := Aelem[x, y, current, 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[current_, i_, j_, n_, m_] := Helem[current, i, j, n, m] =
    NIntegrate[current[x, y, i, j]*Aelem[x, y, current, n, m],
    {y, -Ly, Ly}, {x, -Lx, Lx}, Method -> {"LocalAdaptive"}];

    g[x_,y_,n_,m_]:=Sin[(n*Pi*(x + Lx))/(2*Lx)]*Sin[(m*Pi*(y + Ly))/(2*Ly)]
    curr[x_,y_,n_,m_]:=D[g[x, y, n, m], y];
    {Lx,Ly}={0.1,0.1};
    Helem[curr,1,1,1,1]//Timing

To avoid the message, I've read that I should define the functions using ?NumericQ. This way, the definitions for g, curr and A would turn into g[x_?NumericQ,y_?NumericQ,n_,m_], and so forth.

Using this and calling Helem[g,1,1,1,1] does get rid of the errors! But then, calling Helem[curr,1,1,1,1], the output is as if the functions were not defined. I suppose there is a problem with applying the Derivative D to functions defined using NumericQ, but I don't see a way around this, since I really need to define these functions using some derivative.

Also, any tips for speeding up the integration are more than welcome since it is pretty slow right now!

$\endgroup$
3
$\begingroup$

It seems to be faster to integrate the first integral analyticaly and only the second with NIntegrate.

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

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

curr[x_, y_, n_, m_] = D[g[x, y, n, m], y];

Edit: Answer to your comment, that you have different curr-like-functions.

Since the curr-like-functions do not depend on x1 and y1, integrate first over 1/Sqrt[...] only once and insert curr-functions later.

(int1[x_, y_] = 
    Integrate[
      1/Sqrt[(x - x1)^2 + (y - y1)^2], {y1, -Ly, Ly}, {x1, -Lx, Lx}, 
      Assumptions -> -Lx <= x <= Lx && -Ly <= y <= Ly] // 
     Simplify[#, 
       Assumptions -> -Lx <= x <= Lx && -Ly <= y <= Ly] &) // AbsoluteTiming

(*  {144.4387336, 
     x Log[-1 + 10 y + Sqrt[2] Sqrt[1 - 10 x + 50 x^2 - 10 y + 50 y^2]] - 
     y Log[(1 + 10 x + 
     Sqrt[2] Sqrt[1 + 10 x + 50 x^2 - 10 y + 50 y^2])/(-1 + 10 x + 
     Sqrt[2] Sqrt[1 - 10 x + 50 x^2 - 10 y + 50 y^2])] - 
     x Log[-1 + 10 y + Sqrt[2] Sqrt[1 + 10 x + 50 x^2 - 10 y + 50 y^2]] -
    x Log[1 + 10 y + Sqrt[2] Sqrt[1 - 10 x + 50 x^2 + 10 y + 50 y^2]] +
    y Log[(1 + 10 x + 
    Sqrt[2] Sqrt[1 + 10 x + 50 x^2 + 10 y + 50 y^2])/(-1 + 10 x + 
    Sqrt[2] Sqrt[1 - 10 x + 50 x^2 + 10 y + 50 y^2])] + 
    x Log[1 + 10 y + Sqrt[2] Sqrt[1 + 10 x + 50 x^2 + 10 y + 50 y^2]] + 
    1/10 Log[((1 + 10 x + 
     Sqrt[2] Sqrt[1 + 10 x + 50 x^2 - 10 y + 50 y^2]) (1 + 10 y + 
     Sqrt[2] Sqrt[1 - 10 x + 50 x^2 + 10 y + 50 y^2]) (1 + 10 x + 
     Sqrt[2] Sqrt[1 + 10 x + 50 x^2 + 10 y + 50 y^2]) (1 + 10 y + 
     Sqrt[2] Sqrt[1 + 10 x + 50 x^2 + 10 y + 50 y^2]))/((-1 + 
     10 x + Sqrt[2] Sqrt[
      1 - 10 x + 50 x^2 - 10 y + 50 y^2]) (-1 + 10 y + 
     Sqrt[2] Sqrt[1 - 10 x + 50 x^2 - 10 y + 50 y^2]) (-1 + 10 y +
      Sqrt[2] Sqrt[1 + 10 x + 50 x^2 - 10 y + 50 y^2]) (-1 + 
     10 x + Sqrt[2] Sqrt[1 - 10 x + 50 x^2 + 10 y + 50 y^2]))]}    *)

NIntegrate is then very fast

  Helem[current2_, current1_, opts___] := 
  NIntegrate[current2*current1*int1[x, y], {y, -Ly,Ly},{x,-Lx,Lx},opts]

Helem[curr[x, y, 1, 1], curr[x, y, 1, 1]] // Timing

(*   {0.234, 1.45999- 2.52248*10^-18 I}   *)

Helem[curr[x, y, 1, 1], curr[x, y, 1, 2]] // Timing

This yields error messages, because this integral is zero, since the integrand is symetric in y.

int1[x, y]*curr[x, y, 1, 1]*curr[x, y, 1, 2] + 
int1[x, -y]*curr[x, -y, 1, 1]*curr[x, -y, 1, 2] // ExpToTrig // 
FullSimplify[#, Assumptions -> -Lx < x < Lx && 0 <= y < Ly] &

(*   0   *)

All not zero integrals evaluate very fast

Helem[curr[x, y, 1, 4], curr[x, y, 1, 4]] // Timing

(*   {0.235, 24.6493- 4.03598*10^-17 I}   *)
$\endgroup$
  • $\begingroup$ There are 2 additional complications: - I actually have 4 different kinds of curr-like functions; - The arguments of Helem should really be [n,m,i,j]: it should be like curr[x,y,i,j]*aenm[x,y,n,m] (or equivalent). Should the performance remain similar with these complications? Right now it takes a lot of time, and I'm not sure if I messed up somewhere or if these are simply long calculations. Btw, my computer is slower than yours - takes around 220seconds to run that command you wrote for aenm, but it shouldn't be that much slower. $\endgroup$ – miguel Jun 18 '18 at 9:10
  • $\begingroup$ I made an edit in reaction to your comment. $\endgroup$ – Akku14 Jun 18 '18 at 13:14
  • $\begingroup$ First of all, thank you very much for all the feedback! I seem to be having problems running the first integration command. The kernel simply quits after 5-10 minutes and I am not getting anywhere. Have you encountered any such problems? $\endgroup$ – miguel Jun 19 '18 at 13:00
  • $\begingroup$ Using Mathematica Version 8.0, I had no problems. Try split the int1 into two integrals (int0[x_, y_, y1_] = Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2], {x1, -Lx, Lx}, Assumptions -> -Lx <= x <= Lx && -Ly <= y <= Ly && -Ly <= y1 <= Ly]) // AbsoluteTiming and (int1[x_, y_] = Integrate[int0[x, y, y1], {y1, -Ly, Ly}, Assumptions -> -Lx < x < Lx && -Ly < y < Ly]) // AbsoluteTiming . This yields the same int1, for me even faster. $\endgroup$ – Akku14 Jun 19 '18 at 14:08
  • $\begingroup$ I rechecked my calculations and actually the current in the first integral depends on (x1,y1), so it should also be integrated in the first step. I edited the original post accordingly. Sorry for the imprecision $\endgroup$ – miguel Jun 24 '18 at 10:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.