Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2003
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

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>
  • Prev by Date: RE: Unevaluated
  • Next by Date: ideas on drawing a vector from origin to point
  • Previous by thread: Re: PlotLabel
  • Next by thread: Derivative of a funtion evaluated at a point in 3D