Re: Elliptic integral problems in mma
- To: mathgroup at yoda.physics.unc.edu
- Subject: Re: Elliptic integral problems in mma
- From: keiper
- Date: Wed, 28 Oct 92 10:43:48 CST
> Note that the parameter m > 1 (periodic solution).
> The case m < 1 corresponds to a strictly increasing
> y=y[t] (initial push z0 large enough), and in this
> case there is only one curve
There are indeed two cases and the code defining JacobiAmplitude[ ]
(in StartUp`Elliptic`) is wrong. It tries to deal with the more
complicated case of m < 1 (which it gets right), but messes up on
the much easier case of m > 1. Below is the fix which should be put
in the file StartUp`Elliptic` near line 1008. (Delete the old code.)
I take full responsibility for this error and appologize for any
inconvenience it may have caused.
Jerry B. Keiper
keiper at wri.com
(* ----------------------- JacobiAmplitude ------------------------ *)
JacobiAmplitude[u_, 0] := u;
JacobiAmplitude[u_, 1] := 2 ArcTan[Exp[u]] - Pi/2;
JacobiAmplitude[u_?NumberQ, m_?NumberQ] :=
Module[{am = ArcSin[JacobiSN[u,m]], pi, s, k2, kp, km},
k2 = 2 $EllipK[N[m, Precision[am]]];
s = EllipticF[am, m];
pi = N[Pi, Accuracy[am]+2];
km = Re[(u - s)/k2];
kp = Re[(u + s)/k2];
If[Abs[km - Round[km]] > Abs[kp - Round[kp]],
Round[kp] pi - am,
(* else *)
Round[km] pi + am]
] /; (Precision[{u, m}] < Infinity && Head[m] =!= Complex &&
Abs[m] < 1)
JacobiAmplitude[u_?NumberQ, m_?NumberQ] :=
ArcSin[JacobiSN[u,m]] /; Precision[{u, m}] < Infinity
JacobiAmplitude/: Derivative[1,0][JacobiAmplitude] = JacobiDN
InverseFunction[JacobiAmplitude, 1, 2] = EllipticF
InverseFunction[EllipticF, 1, 2] = JacobiAmplitude