Re: NIntegrate singularity problem v4.0 vs v5.0
- To: mathgroup at smc.vnet.net
- Subject: [mg44471] Re: NIntegrate singularity problem v4.0 vs v5.0
- From: Anton Antonov <antonov at wolfram.com>
- Date: Wed, 12 Nov 2003 08:01:22 -0500 (EST)
- References: <boieu0$o9d$1@smc.vnet.net> <200311110055.TAA25147@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
The evaluation semantics of NIntegrate in V5.0 is changed. (As in NSum, FindRoot,
and FindMinimum.) The integrand is evaluated before proceeding with the numerical
computations. Let's consider an example based on the Mathematica code in the question
(it exhibits the same problem; I have replaced \[Xi] with g):
In[41]:=
Clear[K,Terme,PsiC,Cr]
K[r_,g_]:=(Print["r =!= g"];1/r)/;r=!=g;
K[r_,g_]:=(Print["r === g"];r)/;r===g;
Terme[x_]:=NIntegrate[K[r,0.],{r,-1/2,1/2}]
Terme[0.3]//InputForm
r =!= g
\!\(\*
RowBox[{\(NIntegrate::"inum"\), \(\(:\)\(\ \)\), "\<\"Integrand \\!\\(K[\\(\
\\(r, 0.`\\)\\)]\\) is not numerical at \\!\\({r}\\) = \\!\\({0.`}\\). \\!\\(\
\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", \
ButtonFrame->None, ButtonData:>\\\"NIntegrate::inum\\\"]\\)\"\>"}]\)
Out[45]//InputForm=
NIntegrate[K[r, 0.], {r, -2^(-1), 1/2}]
Note that the definition of the function K with the condition /;r===g has not been
picked up.
The following evaluation is done inside NIntegrate
In[46]:=
Evaluate[K[r,0.]]//InputForm
r =!= g
Out[46]//InputForm=
r^(-1)
It returns as a result
Out[197]//InputForm=
r^(-1)
So if we prevent the evaluation with the the inclusion of _?NumberQ in the function
definitions everything works fine:
In[47]:=
Clear[K,Terme,PsiC,Cr]
K[r_?NumericQ,g_?NumericQ]:=(Print["r =!= g"];1/r)/;r=!=g;
K[r_?NumericQ,g_?NumericQ]:=(Print["r === g"];r)/;r===g;
Terme[x_?NumericQ]:=NIntegrate[K[r,0.],{r,-1/2,1/2}]
Terme[0.3]
r === g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
r =!= g
\!\(\*
RowBox[{\(NIntegrate::"ploss"\), \(\(:\)\(\ \)\), "\<\"Numerical \
integration stopping due to loss of precision. Achieved neither the requested \
PrecisionGoal nor AccuracyGoal; suspect one of the following: highly \
oscillatory integrand or the true value of the integral is 0. If your \
integrand is oscillatory try using the option Method->Oscillatory in \
NIntegrate. \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", \
ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, \
ButtonData:>\\\"NIntegrate::ploss\\\"]\\)\"\>"}]\)
Out[51]=
0.
In[52]:=
Evaluate[K[r,0.]]//InputForm
Out[52]//InputForm=
K[r, 0.]
The Mathematica code in the question had three errors, I guess, from pasting.
On three places I had to put missing multiplication sign, like r*\[Alpha] instead
of r\[Alpha]. I also rewrote the nested NIntegration in the definition
of Cr, NIntegrate[NIntegrate[...]...], into a single one. If we put _?NumberQ in
the function definitions everything works fine. Again, without _?NumberQ specifications,
the case /;r===\[Xi] for the function K is not picked up.
In[53]:=
$Version
Clear[K, Terme, PsiC, Cr]
K[(r_)?NumberQ, (\[Xi]_)?NumberQ, \[Alpha]_, \[Lambda]_] :=
(1/(2*Pi*\[Lambda]*(r^2 - \[Xi]^2)))*(r*\[Alpha]*BesselJ[0, (2*Pi*\[Alpha]*\[Xi])/\[Lambda]]*
BesselJ[1, (2*Pi*r*\[Alpha])/\[Lambda]] - \[Alpha]*\[Xi]*BesselJ[0, (2*Pi*r*\[Alpha])/\[Lambda]]*
BesselJ[1, (2*Pi*\[Alpha]*\[Xi])/\[Lambda]]) /; r =!= \[Xi];
K[(r_)?NumberQ, (\[Xi]_)?NumberQ, \[Alpha]_, \[Lambda]_] :=
(\[Alpha]^2*(BesselJ[0, (2*Pi*\[Alpha]*r)/\[Lambda]]^2 + BesselJ[1, (2*Pi*\[Alpha]*r)/\[Lambda]]^2))/(2*\[Lambda]^2) /;
r === \[Xi];
Terme[(x_)?NumberQ, (\[Alpha]_)?NumberQ, (\[Lambda]_)?NumberQ] :=
(2*Pi)*2*NIntegrate[r*K[r, x, \[Alpha], \[Lambda]], {r, 0, 1/2}];
PsiC[(x_)?NumberQ, (r1_)?NumberQ, (r2_)?NumberQ, (z1_)?NumberQ, (z2_)?NumberQ,
(\[Lambda]_)?NumberQ] := 1 + Terme[x, r2, \[Lambda]]*(E^(I*2*Pi*2*(z2/\[Lambda])) - 1) +
Terme[x, r1, \[Lambda]]*(E^(I*2*Pi*2*(z1/\[Lambda])) - E^(I*2*Pi*2*(z2/\[Lambda])));
Cr[r1_, r2_, z1_, z2_] := NIntegrate[2*Pi*x*Abs[PsiC[x, r1, r2, z1, z2, \[Lambda]]]^2, {x, 0, 0.5},
{\[Lambda], 0.8, 1.2}];
Cr[0.3, 0.5, 2.1, 1.5]
Out[53]=
5.0 for Linux (June 9, 2003)
Out[60]=
0.230245
One can also use the option:
Developer`SetSystemOptions["EvaluateNumericalFunctionArgument"->False]
Related discussion is in
http://support.wolfram.com/mathematica/mathematics/numerics/nsumerror.html
--
Anton Antonov
Wolfram Research, Inc.
> Hi,
> the following code runs on v4.0.1 but not on v5.0, because of the
> changes in NIntegrate function. It seems that v5.0 cannot handle the
> conditional definition of the function K[].
>
> There is a singularity for r=xi and I tried to use the new option in
> v5.0, removing the conditional definition of K[] and giving the
> singularity point r=x directly to NIntegrate :
> \!\(Terme[x_, \[Alpha]_, \[Lambda]_] := \(\((2 \[Pi])\)\^2\)
> NIntegrate[r\ \ K[r, x, \[Alpha], \[Lambda]], {r, 0, x, 1\/2}]\)
> it works fine for this function, but produces an other error in the
> following integrals.
> If someone has ideas about this problem or could help me, I will
> greatly appreciate !
> Thanks,
> Remi
>
> here is the code:
>
> \!\(K[r_, \[Xi]_, \[Alpha]_, \[Lambda]_] := \(1\/\(2\ \[Pi]\
> \[Lambda]\ \
> \((r\^2 - \[Xi]\^2)\)\)\) \((r\ \[Alpha]\ BesselJ[
> 0, \(2\ \[Pi]\ \[Alpha]\ \[Xi]\)\/\[Lambda]]\ BesselJ[
> 1, \(2\ \[Pi]\ r\ \[Alpha]\)\/\[Lambda]] - \[Alpha]\
> \[Xi]\ \
> BesselJ[0, \(2\ \[Pi]\ r\ \[Alpha]\)\/\[Lambda]]\ BesselJ[
> 1, \(2\ \[Pi]\ \[Alpha]\ \[Xi]\)\/\[Lambda]])\) /;
> UnsameQ[r, \[Xi]]\n
> K[r_, \[Xi]_, \[Alpha]_, \[Lambda]_] := \(\[Alpha]\^2\ \((BesselJ[0,
> \(2\ \
> \[Pi]\ \[Alpha]\ r\)\/\[Lambda]]\^2 + BesselJ[1, \(2\ \[Pi]\ \[Alpha]\
> r\)\/\
> \[Lambda]]\^2)\)\)\/\(2\ \[Lambda]\^2\) /; SameQ[r, \[Xi]]\n
> Terme[x_, \[Alpha]_, \[Lambda]_] := \(\((2 \[Pi])\)\^2\)
> NIntegrate[r\ \ K[r, x, \[Alpha], \[Lambda]], {r, 0, 1\/2}]\n
> \(PsiC[x_, r1_, r2_, z1_, z2_, \[Lambda]_] :=
> 1 + Terme[x,
> r2, \[Lambda]] \((\[ExponentialE]\^\(\[ImaginaryI]\ \ 2
> \[Pi]\ 2\
> \ z2/\[Lambda]\) - 1)\) +
> Terme[x,
> r1, \[Lambda]] \((\[ExponentialE]\^\(\[ImaginaryI]\ \ 2
> \[Pi]\ 2\
> \ z1/\[Lambda]\) - \[ExponentialE]\^\(\[ImaginaryI]\ \ 2 \[Pi]\ 2\
> z2/\
> \[Lambda]\))\);\)\n
> \(Cr[r1_, r2_, z1_, z2_] :=
> NIntegrate[
> NIntegrate[
> 2 \[Pi]\ x\ Abs[PsiC[x, r1, r2, z1, z2, \[Lambda]]]\^2, {x,
> 0, .5}], {\[Lambda], .8, 1.2}];\)\n
> Cr[ .3, .5, 2.1, 1.5]\)
- References:
- Re: PlotLabel
- From: Jens-Peer Kuska <kuska@informatik.uni-leipzig.de>
- Re: PlotLabel