MathGroup Archive 2007

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

Search the Archive

Re: Efficient BesselJ and LegendreP Evaluation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg74721] Re: Efficient BesselJ and LegendreP Evaluation
  • From: "dimitris" <dimmechan at yahoo.com>
  • Date: Mon, 2 Apr 2007 06:55:27 -0400 (EDT)
  • References: <eunq6b$7eu$1@smc.vnet.net>

Hi Antonio.

Your post apperared in quite unreadable format (take a look at the
functions arguments).
Anyway from our private communication I was able to see your integrals
you are interested in.

I will what I can do.

Dimitris

=CF/=C7 Antonio =DD=E3=F1=E1=F8=E5:
> 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, =C=
E=B1Lim = =CE=B1Ma=
> x];
> cosA2[(=CE=B1_)?NumericQ] := Sqrt[1 - (nGlass*(Sin[=CE=B1]/nWat=
er))^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[=C=
E=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*D=
Legendre[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}];]



  • Prev by Date: Re: Simplification with Integers assumption
  • Next by Date: Finding unknown parameters using Mathematica
  • Previous by thread: Re: Efficient BesselJ and LegendreP Evaluation
  • Next by thread: Re: Efficient BesselJ and LegendreP Evaluation