MathGroup Archive 2006

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

Search the Archive

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


  • Prev by Date: Re: Question about UpSet vs SetDelayed
  • Next by Date: JLink / VTK problem
  • Previous by thread: Re: Performance improvement needed - Help.
  • Next by thread: Re: How to hide a cell?