Re: Efficient BesselJ and LegendreP Evaluation
- To: mathgroup at smc.vnet.net
- Subject: [mg74788] Re: Efficient BesselJ and LegendreP Evaluation
- From: "Antonio" <aneves at gmail.com>
- Date: Thu, 5 Apr 2007 04:11:19 -0400 (EDT)
- References: <eunq6b$7eu$1@smc.vnet.net><euqnj9$88g$1@smc.vnet.net>
On 2 abr, 11:57, "dimitris" <dimmec... at yahoo.com> wrote:
> 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 hasBesselJ, 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}];]- Ocultar texto entre aspas -
>
> - Mostrar texto entre aspas -
Thank you Dimitris, any help is appreciated. Glad also if I could
receive some feedback from the readers of Mathematica.Group.
Anyone ever encountered a similar problem?
Antonio