Re: Performance improvement needed - Help.
- To: mathgroup at smc.vnet.net
- Subject: [mg63824] Re: Performance improvement needed - Help.
- From: "antononcube" <antononcube at gmail.com>
- Date: Mon, 16 Jan 2006 02:16:51 -0500 (EST)
- References: <dpiotv$l11$1@smc.vnet.net><dq56cj$b0k$1@smc.vnet.net> <dq7tpa$3nq$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In my previous message I suggested the definition:
Clear[I1];
I1[n_, m_, X_, c_] :=
Block[{func},
func[a_?NumericQ] :=
Sqrt[Cos[a]] Exp[-Sin[a]^2 c^2] (LegendreP[n, m,
Cos[a]] (m^2 BesselJ[m, X Sin[a]](1 - Cos[a])/(X
Sin[a]) +
m Cos[a] BesselJ[m + 1,
X Sin[a]]) + (m BesselJ[m, X Sin[a]]/(X Sin[a]) -
BesselJ[m + 1, X Sin[a]])(n + m)(n - m + 1)Sin[a]
LegendreP[
n, m - 1, Cos[a]]);
NIntegrate[func[a], {a, 0, 1.1}]]
(Here I have removed the options given to NIntegrate since they are not
necessary for this problem.)
In this note I would like to clarify why it works.
Let's define func[] in the global context without the _?NumericQ
argument restriction.
Clear[func];
func[a_]:=Sqrt[
Cos[a]] Exp[-Sin[a]^2 c^2] (LegendreP[n,m,
Cos[a]] (m^2 BesselJ[m,X Sin[a]](1-Cos[a])/(X Sin[a])+
m Cos[a] BesselJ[m+1,
X Sin[a]])+(m BesselJ[m,X Sin[a]]/(X Sin[a])-
BesselJ[m+1,X Sin[a]])(n+m)(n-m+1)Sin[a]
LegendreP[n,m-1,
Cos[a]]);
If we plot func[x] with
Plot[func[a],{a,0,1.1}]
we can see that the integrand oscillates but not much, so NIntegrate
shouldn't have much problem integrating it. Nevertheless NIntegrate did
have problems (i.e. it was slow) as it was indicated in the openning
post.
This is how func[a] looks:
In[9]:= func[a]
Out[9]=
(Sqrt[Cos[a]]*((m*BesselJ[1 + m, X*Sin[a]]*Cos[a] + (m^2*BesselJ[m,
X*Sin[a]]*(1 - Cos[a])*Csc[a])/X)*
LegendreP[n, m, Cos[a]] + (1 - m + n)*(m + n)*(-BesselJ[1 + m,
X*Sin[a]] +
(m*BesselJ[m, X*Sin[a]]*Csc[a])/X)*LegendreP[n, -1 + m,
Cos[a]]*Sin[a]))/E^(c^2*Sin[a]^2)
But if we give values to n, m, X, and c,
In[13]:=
n = 40; m = 1;
X = 1; c = 1;
func[a] will look quite different:
In[17]:= func[a]
Out[17]=
(Sqrt[Cos[a]]*(-((1/34359738368)*(205*Sqrt[1 -
Cos[a]^2]*(-34461632205*Cos[a] + 9385051170495*Cos[a]^3 -
760189144810095*Cos[a]^5 + 28923386985869805*Cos[a]^7 -
629887094358942420*Cos[a]^9 +
8761156857901653660*Cos[a]^11 - 83343312673884961740*Cos[a]^13
+ 567528272017407120420*Cos[a]^15 -
2854333368087547576230*Cos[a]^17 +
10833113192332271210370*Cos[a]^19 -
31467614511060406849170*Cos[a]^21 +
70522282323206524440630*Cos[a]^23 -
122238622693557975697092*Cos[a]^25 +
163333088442389431914348*Cos[a]^27 -
166551474421549814809212*Cos[a]^29 +
127152200902473514531764*Cos[a]^31 -
70319020196064898188021*Cos[a]^33 +
26591226124562356457655*Cos[a]^35 -
6148721956730634976695*Cos[a]^37 +
655531760569123027205*Cos[a]^39)*
(BesselJ[2, Sin[a]]*Cos[a] + BesselJ[1, Sin[a]]*(1 -
Cos[a])*Csc[a]))) +
(1/34359738368)*(205*(34461632205 - 28258538408100*Cos[a]^2 +
3847870979902950*Cos[a]^4 -
207785032914759300*Cos[a]^6 + 5929294332103310025*Cos[a]^8 -
103301483474866556880*Cos[a]^10 +
1197358103913226000200*Cos[a]^12 -
9763073770369381232400*Cos[a]^14 +
58171647881784229843050*Cos[a]^16 -
260061484647976556945400*Cos[a]^18 +
888315281771246239250340*Cos[a]^20 -
2345767627188139419665400*Cos[a]^22 +
4819022625419112503443050*Cos[a]^24 -
7710436200670580005508880*Cos[a]^26 +
9566652323054238154983240*Cos[a]^28 -
9104813935044723209570256*Cos[a]^30 +
6516550296251767619752905*Cos[a]^32 -
3391858621221953912598660*Cos[a]^34 +
1211378079007840683070950*Cos[a]^36 -
265365894974690562152100*Cos[a]^38 +
26876802183334044115405*Cos[a]^40)*(-BesselJ[2, Sin[a]] +
BesselJ[1, Sin[a]]*Csc[a])*Sin[a])))/
E^Sin[a]^2
This last expression although mathematically equivalent to the one we
got first by evaluating func[a] is prone to numerical instabilities
when evaluated for numerical a. This can be seen by comparing these two
plots:
Plot[func[a], {a, 0.03, 0.05}]
Plot[func[a] // Evaluate, {a, 0.03, 0.05}]
In the function I1 NIntegrate will use the first expression of func[a]
since its symbolic evaluation is prevented with the pattern. _?NumericQ