Services & ResourcesWolfram ForumsMathGroup Archive

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

Efficient BesselJ and LegendreP Evaluation

• To: mathgroup at smc.vnet.net
• Subject: [mg74706] Efficient BesselJ and LegendreP Evaluation
• From: "Antonio" <aneves at gmail.com>
• Date: Sun, 1 Apr 2007 04:18:12 -0400 (EDT)

Dear group members,

I have a function that I want to integrate for various variable
parameter in the most precise and fastest way using mathematica
instead of external compile C program. Below is a code to evaluate a
numerical integral, which has BesselJ, LegendreP and their derivative.
These needs to be evaluated for a set of {n,m}. Below in this example
I have {n,1,20} and {m,0,n}.

There is also a singularity issue to be solved:"NIntegrate::singd :
NIntegrate's singularity handling has failed at point
{=CE=B1}={6.007268658180993`*^-9} for the specified precision goal. Try
using larger values for any of $MaxExtraPrecision or the options
WorkingPrecision, or SingularityDepth and MaxRecursion."

Any comments are welcomed.

Antonio



Clear["Global`*"];
Off[General::spell];
Off[General::spell1];
=CE=BB = 0.8;
w = 2500;
NA = 1.25;
f = 1700;
d = 50.;
nGlass = 1.4883 + 1.9023/(10^2*=CE=BB^2) - 4.25016/(10^3*=CE=BB^4);
nWater = 1.3253 + 2.7553/(10^3*=CE=BB^2) + 3.779/(10^5*=CE=BB^4);
nAir = 1.0001;
=CE=B1Max = ArcSin[NA/nWater];
ko = (2*Pi*nAir)/=CE=BB;
Nmax = 20;
GridNM = Table[{n, m}, {n, 1, Nmax}, {m, 0, n}];
=CE=B1Crit = ArcSin[nWater/nGlass];
If[=CE=B1Crit < =CE=B1Max, =CE=B1Lim = =CE=B1Crit, =CE=B1Lim = =CE=B1Ma=
x];
cosA2[(=CE=B1_)?NumericQ] := Sqrt[1 - (nGlass*(Sin[=CE=B1]/nWater))^2];
ts[(=CE=B1_)?NumericQ] := 2*nGlass*(Cos[=CE=B1]/(nGlass*Cos[=CE=B1] +
nWater*cosA2[=CE=B1]));
tp[(=CE=B1_)?NumericQ] := 2*nGlass*(Cos[=CE=B1]/(nWater*Cos[=CE=B1] +
nGlass*cosA2[=CE=B1]));
DBessel[m_, x_] := (1/2)*(BesselJ[-1 + m, x] - BesselJ[1 + m, x]);
DLegendre[n_, m_, x_] := ((-1 - n)*x*LegendreP[n, m, x] + (1 - m +
n)*LegendreP[1 + n, m, x])/(-1 + x^2);
ITM[(n_Integer)?Positive, (m_Integer)?NonNegative, zo_Real,
=CF=81o_Real] :=
   Block[{func}, func[(=CE=B1_)?NumericQ] := Sqrt[Cos[=CE=B1]]*Exp[-(f*(S=
in[=CE=B1]/
w))^2]*Exp[(-I)*ko*nWater*zo*cosA2[=CE=B1]]*
       Exp[(-I)*ko*d*(nGlass*Cos[=CE=B1] -
nWater*cosA2[=CE=B1])]*(m^2*ts[=CE=B1]*LegendreP[n, m, cosA2[=CE=B1]]*(Bess=
elJ[m,
ko*nGlass*=CF=81o*Sin[=CE=B1]]/(ko*nGlass*=CF=81o*Sin[=CE=B1])) -
        tp[=CE=B1]*(nGlass/nWater)^2*Sin[=CE=B1]^2*DBessel[m,
ko*nGlass*=CF=81o*Sin[=CE=B1]]*DLegendre[n, m, cosA2[=CE=B1]] +
        I*m*(ts[=CE=B1]*DBessel[m, ko*nGlass*=CF=81o*Sin[=CE=B1]]*LegendreP=
[n, m,
cosA2[=CE=B1]] - tp[=CE=B1]*(nGlass/nWater)^2*Sin[=CE=B1]^2*DLegendre[n, m,
cosA2[=CE=B1]]*
           (BesselJ[m, ko*nGlass*=CF=81o*Sin[=CE=B1]]/(ko*nGlass*=CF=81o*Si=
n[=CE=B1]))));
NIntegrate[func[=CE=B1], {=CE=B1, 0, =CE=B1Lim}, Compiled -> True, Method ->
DoubleExponential,
      MaxRecursion -> 1000, SingularityDepth -> 1000]];
zo = Random[Real, {-10., 10.}];
=CF=81o = Random[Real, {0., 10.}];
Timing[MI1 = Apply[ITM[#1, #2, zo, =CF=81o] & , GridNM, {2}];]

• Follow-Ups:
  • Re: Efficient BesselJ and LegendreP Evaluation
    • From: "Antonio Neves" <aneves@gmail.com>