Re: Numerical Integration
- To: mathgroup at smc.vnet.net
- Subject: [mg67384] Re: Numerical Integration
- From: "antononcube" <antononcube at gmail.com>
- Date: Wed, 21 Jun 2006 02:12:57 -0400 (EDT)
- References: <e784ps$fr0$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
The second integration fails because FirstIntegration[0.] fails.
FirstIntegration[0.] fails because
the function int[x_,y_] is Indeterminate for int[0,0]. If you specify a
singular point 0 for the second integration it works.
In[8]:= FirstIntegration[0.]
NIntegrate::inum:
4 2
-(- - -------------------------)
3 2
Sqrt[0.01 + (3.2 - kkx) ]
Integrand -------------------------------- + <<1>> is not numerical
at
2
4 (9.6 - 3 kkx) Pi
{kkx} = {0.}.
2 2
Out[8]= NIntegrate[int[kkx, 0.], {kkx, -Sqrt[1 - 0. ], Sqrt[1 - 0. ]},
> MaxRecursion -> 20]
In[9]:= int[0, 0]
1
Power::infy: Infinite expression - encountered.
0
Infinity::indet: Indeterminate expression 0 ComplexInfinity
encountered.
1
Power::infy: Infinite expression - encountered.
0
Infinity::indet: Indeterminate expression 0 ComplexInfinity
encountered.
1
Power::infy: Infinite expression - encountered.
0
General::stop: Further output of Power::infy
will be suppressed during this calculation.
Infinity::indet: Indeterminate expression 0 ComplexInfinity
encountered.
General::stop: Further output of Infinity::indet
will be suppressed during this calculation.
Out[9]= Indeterminate
In[10]:= NIntegrate[FirstIntegration[kky], {kky, -1, 0, 1}]
Out[10]= 1.88673 10^-6
Anton Antonov,
Wolfram Research, Inc.
Stefano Chesi wrote:
> Hi,
> I have to calculate a double integral of a function (the definition is
> long and I copy it
> at the end of the mail). I can plot the result of the first integration:
>
> INPUT:
> FirstIntegration[kky_?NumberQ]:=
>
> NIntegrate[int[kkx,kky],{kkx,-Sqrt[1-kky^2],Sqrt[1-kky^2]},MaxRecursion->20]
> Plot[FirstIntegration[kky],{kky,-1,1}]
> OUTPUT: Graphics
>
> But the second integration doesn't work:
>
> INPUT:
> NIntegrate[FirstIntegration[kky],{kky,-1,1}]
> OUTPUT: NIntegrate::inum: .... is not numerical at {kky} = {0.}
> NIntegrate[FirstIntegration[kky],{kky,-1,1}]
>
> The actual code is slightly more complicated, and I need to do the
> integration
> in two different steps (sometimes in the first integration I have poles
> and I need
> to take the principal value), in the way it is coded above.
>
> Thank you very much,
> Stefano
>
>
> Here's the function definition. If I make it a simple function, e.g. by
> adding " ; kky^2 "
> before the last parenthesis "]", the second integration works. So, I
> don't understand
> what's the problem.
>
> aa=0.1;
>
> kx=0.2;
> ky=0.1;
> q=3;
>
>
> int[kkx_,kky_]:=
>
> Module[
>
> {nn,n1,n2,n3,n4,c1,c2,c3,c4,s1,s2,s3,s4},
>
> nn=Sqrt[(q+kx-kkx)^2+(ky-kky)^2];
>
> n1=-Sqrt[kx^2+ky^2];
> n2=Sqrt[(kx+q)^2+ky^2];
> n3=-Sqrt[kkx^2+kky^2];
> n4=Sqrt[(kkx-q)^2+kky^2];
>
> c1=((kx+q)kx+ky ky)/n1/n2;
> c2=(kkx(kx+q)+kky ky)/n2/n3;
> c3=((kkx-q)kkx+kky kky)/n3/n4;
> c4=(kx(kkx-q)+ky kky)/n4/n1;
>
> s1=(ky kx-(kx+q) ky)/n1/n2;
> s2=(kky(kx+q)-kkx ky)/n2/n3;
> s3=(kky kkx-(kkx-q) kky)/n3/n4;
> s4=(ky (kkx-q)-kx kky)/n4/n1;
>
> (
> (c1 c1 c3 c3/q- c1 c2 c3
> c4/nn)/(q^2+kx*q-kkx*q-aa(n1+n2+n3+n4))+
>
> (s1 s1 c3 c3/q+ s1 c2 c3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1+n2+n3+n4))+
> (s1 s1 c3 c3/q+s1 s2 c3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1-n2+n3+n4))+
> (c1 c1 s3 s3/q+ c1 s2 s3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1+n2-n3+n4))+
> (c1 c1 s3 s3/q+ c1 c2 s3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1+n2+n3-n4))+
>
> (c1 c1 c3 c3/q+ c1 s2 c3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1-n2+n3+n4))+
> (s1 s1 s3 s3/q- s1 s2 s3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1+n2-n3+n4))+
> (s1 s1 s3 s3/q+ s1 c2 s3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1+n2+n3-n4))+
> (s1 s1 s3 s3/q+ s1 c2 s3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1-n2-n3+n4))+
> (s1 s1 s3 s3/q- s1 s2 s3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1-n2+n3-n4))+
> (c1 c1 c3 c3/q+ c1 s2 c3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1+n2-n3-n4))+
>
> (s1 s1 c3 c3/q+ s1 c2 c3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(+n1-n2-n3-n4))+
> (s1 s1 c3 c3/q+ s1 s2 c3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1+n2-n3-n4))+
> (c1 c1 s3 s3/q+ c1 s2 s3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1-n2+n3-n4))+
> (c1 c1 s3 s3/q+ c1 c2 s3 s4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1-n2-n3+n4))+
>
> (c1 c1 c3 c3/q- c1 c2 c3 c4/nn)/(q^2+kx*q-kkx*q-
> aa(-n1-n2-n3-n4))
>
> )/4/Pi^2-
>
> (4/q - 2/nn)/(q^2 + kx*q - kkx*q)/4/Pi^2
>
> ]