2
$\begingroup$

I am trying to compute an integral of the form

$\int_{\frac{1}{10000}}^{\frac{1}{1000}}g(x)x e^{3ix}$

in Mathematica both numerically and analytically. $g(x)$ is equal to

g[x_] := (1.647918829487756` + 0.0004377541512921337` I)-(5.950382035839695`*^-21 +1.0244561383672523`*^-17 I)/x^3 + ( 4.917597072242726`*^-16 + 6.250640697230879`*^-13 I)/x^2-( 1.7623301652915415`*^-11 + 1.8809797793209515`*^-8 I)/x -(5.184867552096762-20.947979274232033` I) x - (107.18849192181267` + 17075.00566589623` I) x^2 + (142847.93815916454` +  3.083551530852119`*^7 I) x^3 - (2.2921426862077513`*^8 + 
 5.329598975685562`*^10 I) x^4 + (3.2507619599780566`*^11 + 
 7.42469714410255`*^13 I) x^5 - (3.5772293756902906`*^14 + 
 7.82535940524534`*^16 I) x^6 + (2.853115060240398`*^17 + 
 5.938191548262084`*^19 I) x^7 - (1.5363326543684728`*^20 + 
 3.0415782075611755`*^22 I) x^8 + (4.969151746509189`*^22 + 
 9.378105083129474`*^24 I) x^9 - (7.267395382752606`*^24 + 
 1.3113949262912268`*^27 I) x^10

I find $8.13965*10^-7 + 6.95987*10^-9 I$ for the following integral

NIntegrate[g[x]*x*Exp[I*x*3], {x, 1/10000, 1/1000}]

but when I use Integrate instead of NIntegrate, I get 0. I am wondering if someone can explain this to me. Thank you!

I have one more question. The reason I asked this question was that I was trying to compute a triple integral including the integral I mentioned before. The triple integral is

PVtest[ts_] := NIntegrate[(g[z]*(I*z)*Exp[I*(-z)*4]*Exp[I*z*ts])*
   NIntegrate[
g[y]*(I*y)*Exp[I*(-y)*4]*(Exp[I*y*ts]/(y - x)), {y, 1/10000, x, 
 1/1000},PrincipalValue -> True], {z, 1/10000, 1/1000}, {x, 1/10000, 1/1000}, Method -> {Automatic, "SymbolicProcessing" -> 0}]

I need to evaluate this principal value integral for different values of $ts$, For example, 3, 10, 15, ... Mathematica is very slow in computing the integral. So I wanted to integrate over variable $z$ first and then numerically compute the rest of the integrals. I asked this question before. Somebody told me you could take the integral over $x$ first analytically and use the principal value method, I did it and the answer for $ts=3$ was close to what I had found using NIntegrate of the triple integral, but when I increase $ts$ a little, the answer was different. This makes sense because I think the principal value integral is for variable $y$ and we cannot change the order of integration because of the singularity. I am wondering if you know how to accelerate the numerical computation of the aforementioned triple integral. Thank you!

$\endgroup$
1
  • 2
    $\begingroup$ Neither is necessarily correct although I might believe NIntegrate did something sensible. The integrand is ripe for cancellation error. And possibly it is ill-conditioned, given the scales of numbers and polynomial degree involved $\endgroup$ Jun 12, 2023 at 3:27

1 Answer 1

1
$\begingroup$

I think Integrate is wrong. It seems due to use floating point with Integrate which is not good idea any way.

g[x_] := (1.647918829487756` + 
    0.0004377541512921337` I) - (5.950382035839695`*^-21 + 
     1.0244561383672523`*^-17 I)/
   x^3 + (4.917597072242726`*^-16 + 6.250640697230879`*^-13 I)/
   x^2 - (1.7623301652915415`*^-11 + 1.8809797793209515`*^-8 I)/
   x - (5.184867552096762 - 
     20.947979274232033` I) x - (107.18849192181267` + 
     17075.00566589623` I) x^2 + (142847.93815916454` + 
     3.083551530852119`*^7 I) x^3 - (2.2921426862077513`*^8 + 
     5.329598975685562`*^10 I) x^4 + (3.2507619599780566`*^11 + 
     7.42469714410255`*^13 I) x^5 - (3.5772293756902906`*^14 + 
     7.82535940524534`*^16 I) x^6 + (2.853115060240398`*^17 + 
     5.938191548262084`*^19 I) x^7 - (1.5363326543684728`*^20 + 
     3.0415782075611755`*^22 I) x^8 + (4.969151746509189`*^22 + 
     9.378105083129474`*^24 I) x^9 - (7.267395382752606`*^24 + 
     1.3113949262912268`*^27 I) x^10

anti=Integrate[SetAccuracy[g[x],Infinity]*x*Exp[I*x*3],x];
(*OK to do this, as integrand has no poles in between and continous*)
N[Limit[anti,x->1/1000]-Limit[anti,x->1/10000],16]

Mathematica graphics

You also do not need to take the limit manually. This works

anti = Integrate[SetAccuracy[g[x], Infinity]*x*Exp[I*x*3], {x, 1/10000, 1/1000}];
N[anti, 16]

Mathematica graphics

The zero result shows up only when using real numbers:

anti = Integrate[g[x]*x*Exp[I*x*3], {x, 1/10000, 1/1000}]

Mathematica graphics

Compare to when using NIntegrate

NIntegrate[g[x]*x*Exp[I*x*3],{x,1/10000,1/1000}]

Mathematica graphics

But you do not need all of this. One can see that the total areas under the curve of the integrand can not be zero (of the real and imaginary parts).

Here is plot of of the integrand, the real part and its imaginary part

Grid[{{Plot[Re[g[x]*x*Exp[I*x*3]],{x,1/10000,1/1000},
  AxesOrigin->{0,0}],Plot[Im[g[x]*x*Exp[I*x*3]], 
  {x,1/10000,1/1000},AxesOrigin->{0,0}]}}]

Mathematica graphics

They are both positive and do not cross the x-axis. so total are can not be zero. By ballpark estimate the above total area is

1/2*(1/1000-1/10000)*(0.0015-0.0001)+ I*(1/2*(1/1000 - 1/10000)*(0.00002 - 3*10^(-6)))

Mathematica graphics

Which is close to what NIntegrate gave. And this is just by looking at the Plot. So zero result is wrong.

The bottom line: Avoid using non-exact numbers with exact solvers. It can lead to less accurate results in some cases.

$\endgroup$
3
  • $\begingroup$ Thank you so much for your detailed answer! @Nasser $\endgroup$
    – HadamardN2
    Jun 12, 2023 at 1:00
  • $\begingroup$ I have one more question. I couldn't type it here so I added that to the main question. Could you please help me with that? Thanks! For some reason, I cannot mention you at the beginning of my comment. Sorry about that. @Nasser $\endgroup$
    – HadamardN2
    Jun 12, 2023 at 13:46
  • $\begingroup$ @HadamardN2 it looks like your second question you added is not really related to the first one. It will be better to make that separate question, it will get more attention. $\endgroup$
    – Nasser
    Jun 12, 2023 at 19:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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