Math Is Fun Forum

  Discussion about math, puzzles, games and fun.   Useful symbols: ÷ × ½ √ ∞ ≠ ≤ ≥ ≈ ⇒ ± ∈ Δ θ ∴ ∑ ∫ • π ƒ -¹ ² ³ °

You are not logged in.

#1 2025-10-09 11:19:02

lanxiyu
Member
Registered: 2022-05-10
Posts: 65

Mathematica test for Jacobi elliptic functions

Module[{Δ, m, k, kp, u, v, d, ϕu, ϕv, m1, am, sn, cn, dn, tn, Ε},
	If[Context[JacobiEpsilon] =!= "System`",
		JacobiEpsilon = EllipticE[JacobiAmplitude[#1, #2], #2] &;
	];

	Δ = Sqrt[1 - m*Sin[#]^2] &;
	am = JacobiAmplitude[#, m] &;
	sn = JacobiSN[#, m] &;
	cn = JacobiCN[#, m] &;
	dn = JacobiDN[#, m] &;
	tn = JacobiSC[#, m] &;
	Ε = JacobiEpsilon[#, m] &;

	m = RandomReal[WorkingPrecision -> 50];
	k = Sqrt[m];
	kp = Sqrt[1 - m];
	{u, v} = RandomReal[{0, EllipticK[m]}, 2, WorkingPrecision -> 50];
	{ϕu, ϕv} = JacobiAmplitude[{u, v}, m];

	(* Addition formulas *)
	d = 1 - m*Sin[ϕu]^2*Sin[ϕv]^2;

	Assert[sn[u + v] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕv] + Sin[ϕv]*Cos[ϕu]*Δ[ϕu])/d];
	Assert[sn[u - v] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕv] - Sin[ϕv]*Cos[ϕu]*Δ[ϕu])/d];
	Assert[cn[u + v] == (Cos[ϕu]*Cos[ϕv] - Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv])/d];
	Assert[cn[u - v] == (Cos[ϕu]*Cos[ϕv] + Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv])/d];
	Assert[dn[u + v] == (Δ[ϕu]*Δ[ϕv] - m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv])/d];
	Assert[dn[u - v] == (Δ[ϕu]*Δ[ϕv] + m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv])/d];
	Assert[sn[u + v]*sn[u - v] == (Sin[ϕu]^2 - Sin[ϕv]^2)/d];
	Assert[sn[u + v]*cn[u - v] == (Sin[ϕu]*Cos[ϕu]*Δ[ϕv] + Sin[ϕv]*Cos[ϕv]*Δ[ϕu])/d];
	Assert[sn[u + v]*dn[u - v] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕu] + Sin[ϕv]*Cos[ϕu]*Δ[ϕv])/d];
	Assert[cn[u + v]*sn[u - v] == (Sin[ϕu]*Cos[ϕu]*Δ[ϕv] - Sin[ϕv]*Cos[ϕv]*Δ[ϕu])/d];
	Assert[cn[u + v]*cn[u - v] == (1 - Sin[ϕu]^2 - Sin[ϕv]^2 + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
	Assert[cn[u + v]*dn[u - v] == (Cos[ϕu]*Cos[ϕv]*Δ[ϕu]*Δ[ϕv] - (1 - m)*Sin[ϕu]*Sin[ϕv])/d];
	Assert[dn[u + v]*sn[u - v] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕu] - Sin[ϕv]*Cos[ϕu]*Δ[ϕv])/d];
	Assert[dn[u + v]*cn[u - v] == (Cos[ϕu]*Cos[ϕv]*Δ[ϕu]*Δ[ϕv] + (1 - m)*Sin[ϕu]*Sin[ϕv])/d];
	Assert[dn[u + v]*dn[u - v] == (1 - m*Sin[ϕu]^2 - m*Sin[ϕv]^2 + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];

	Assert[sn[u + v] + sn[u - v] ==  2*Sin[ϕu]*Cos[ϕv]*Δ[ϕv]/d];
	Assert[sn[u + v] - sn[u - v] ==  2*Sin[ϕv]*Cos[ϕu]*Δ[ϕu]/d];
	Assert[cn[u + v] + cn[u - v] ==  2*Cos[ϕu]*Cos[ϕv]/d];
	Assert[cn[u + v] - cn[u - v] == -2*Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv]/d];
	Assert[dn[u + v] + dn[u - v] ==  2*Δ[ϕu]*Δ[ϕv]/d];
	Assert[dn[u + v] - dn[u - v] == -2*m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv]/d];
	Assert[sn[u + v] + cn[u - v] ==  (Cos[ϕu] + Sin[ϕu]*Δ[ϕv])*(Cos[ϕv] + Sin[ϕv]*Δ[ϕu])/d];
	Assert[sn[u + v] - cn[u - v] == -(Cos[ϕu] - Sin[ϕu]*Δ[ϕv])*(Cos[ϕv] - Sin[ϕv]*Δ[ϕu])/d];
	Assert[k*sn[u + v] + dn[u - v] ==  (Δ[ϕu] + k*Sin[ϕu]*Cos[ϕv])*(Δ[ϕv] + k*Sin[ϕv]*Cos[ϕu])/d];
	Assert[k*sn[u + v] - dn[u - v] == -(Δ[ϕu] - k*Sin[ϕu]*Cos[ϕv])*(Δ[ϕv] - k*Sin[ϕv]*Cos[ϕu])/d];
	Assert[Sin[am[u + v] + am[u - v]] == Sin[2*ϕu]*Δ[ϕv]/d];
	Assert[Sin[am[u + v] - am[u - v]] == Sin[2*ϕv]*Δ[ϕu]/d];
	Assert[Cos[am[u + v] + am[u - v]] == (Cos[2*ϕu] + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
	Assert[Cos[am[u + v] - am[u - v]] == (Cos[2*ϕv] + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];

	Assert[(1 + sn[u + v])*(1 + sn[u - v]) == (Cos[ϕv] + Sin[ϕu]*Δ[ϕv])^2/d];
	Assert[(1 - sn[u + v])*(1 - sn[u - v]) == (Cos[ϕv] - Sin[ϕu]*Δ[ϕv])^2/d];
	Assert[(1 + sn[u + v])*(1 - sn[u - v]) == (Cos[ϕu] + Sin[ϕv]*Δ[ϕu])^2/d];
	Assert[(1 - sn[u + v])*(1 + sn[u - v]) == (Cos[ϕu] - Sin[ϕv]*Δ[ϕu])^2/d];
	Assert[(1 + k*sn[u + v])*(1 + k*sn[u - v]) == (Δ[ϕv] + k*Sin[ϕu]*Cos[ϕv])^2/d];
	Assert[(1 - k*sn[u + v])*(1 - k*sn[u - v]) == (Δ[ϕv] - k*Sin[ϕu]*Cos[ϕv])^2/d];
	Assert[(1 + k*sn[u + v])*(1 - k*sn[u - v]) == (Δ[ϕu] + k*Sin[ϕv]*Cos[ϕu])^2/d];
	Assert[(1 - k*sn[u + v])*(1 + k*sn[u - v]) == (Δ[ϕu] - k*Sin[ϕv]*Cos[ϕu])^2/d];
	Assert[(1 + cn[u + v])*(1 + cn[u - v]) == (Cos[ϕu] + Cos[ϕv])^2/d];
	Assert[(1 - cn[u + v])*(1 - cn[u - v]) == (Cos[ϕu] - Cos[ϕv])^2/d];
	Assert[(1 + cn[u + v])*(1 - cn[u - v]) == (Sin[ϕu]*Δ[ϕv] - Sin[ϕv]*Δ[ϕu])^2/d];
	Assert[(1 - cn[u + v])*(1 + cn[u - v]) == (Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu])^2/d];
	Assert[(1 + dn[u + v])*(1 + dn[u - v]) == (Δ[ϕu] + Δ[ϕv])^2/d];
	Assert[(1 - dn[u + v])*(1 - dn[u - v]) == (Δ[ϕu] - Δ[ϕv])^2/d];
	Assert[(1 + dn[u + v])*(1 - dn[u - v]) == m*Sin[ϕu - ϕv]^2/d];
	Assert[(1 - dn[u + v])*(1 + dn[u - v]) == m*Sin[ϕu + ϕv]^2/d];

	Assert[dn[u + v] + k*cn[u + v] == (Δ[ϕu]*Δ[ϕv] + k*Cos[ϕu]*Cos[ϕv])/(1 + k*Sin[ϕu]*Sin[ϕv])];
	Assert[dn[u + v] - k*cn[u + v] == (Δ[ϕu]*Δ[ϕv] - k*Cos[ϕu]*Cos[ϕv])/(1 - k*Sin[ϕu]*Sin[ϕv])];
	Assert[dn[u - v] + k*cn[u - v] == (Δ[ϕu]*Δ[ϕv] + k*Cos[ϕu]*Cos[ϕv])/(1 - k*Sin[ϕu]*Sin[ϕv])];
	Assert[dn[u - v] - k*cn[u - v] == (Δ[ϕu]*Δ[ϕv] - k*Cos[ϕu]*Cos[ϕv])/(1 + k*Sin[ϕu]*Sin[ϕv])];

	Assert[Sin[(am[u + v] + am[u - v])/2] == Sin[ϕu]*Δ[ϕv]/Sqrt[d]];
	Assert[Sin[(am[u + v] - am[u - v])/2] == Sin[ϕv]*Δ[ϕu]/Sqrt[d]];
	Assert[Cos[(am[u + v] + am[u - v])/2] == Cos[ϕu]/Sqrt[d]];
	Assert[Cos[(am[u + v] - am[u - v])/2] == Cos[ϕv]/Sqrt[d]];
	Assert[Tan[(am[u + v] + am[u - v])/2] == Tan[ϕu]*Δ[ϕv]];
	Assert[Tan[(am[u + v] - am[u - v])/2] == Tan[ϕv]*Δ[ϕu]];

	Assert[Tan[am[u + v]/2] == (Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu])/(Cos[ϕu] + Cos[ϕv])];
	Assert[Tan[am[u - v]/2] == (Sin[ϕu]*Δ[ϕv] - Sin[ϕv]*Δ[ϕu])/(Cos[ϕu] + Cos[ϕv])];
	Assert[sn[u + v]/(1 + dn[u + v]) == Sin[ϕu + ϕv]/(Δ[ϕu] + Δ[ϕv])];
	Assert[sn[u - v]/(1 + dn[u - v]) == Sin[ϕu - ϕv]/(Δ[ϕu] + Δ[ϕv])];

	Assert[tn[u + v] == (Tan[ϕu]*Δ[ϕv] + Tan[ϕv]*Δ[ϕu])/(1 - Tan[ϕu]*Tan[ϕv]*Δ[ϕu]*Δ[ϕv])];
	Assert[tn[u - v] == (Tan[ϕu]*Δ[ϕv] - Tan[ϕv]*Δ[ϕu])/(1 + Tan[ϕu]*Tan[ϕv]*Δ[ϕu]*Δ[ϕv])];

	Assert[Ε[u + v] == EllipticE[ϕu, m] + EllipticE[ϕv, m] - m*Sin[ϕu]*Sin[ϕv]*sn[u + v]];
	Assert[Ε[u - v] == EllipticE[ϕu, m] - EllipticE[ϕv, m] + m*Sin[ϕu]*Sin[ϕv]*sn[u - v]];

	Clear[d];

	(* Double angle formulas *)
	d = 1 - m*Sin[ϕu]^4;
	Assert[sn[2*u] == Sin[2*ϕu]*Δ[ϕu]/d];
	Assert[cn[2*u] == (Cos[2*ϕu] + m*Sin[ϕu]^4)/d];
	Assert[dn[2*u] == (1 - 2*m*Sin[ϕu]^2 + m*Sin[ϕu]^4)/d];
	Assert[Ε[2*u] == 2*(EllipticE[ϕu, m] - m*Sin[ϕu]^3*Cos[ϕu]*Δ[ϕu]/d)];
	Clear[d];

	(* Triple angle formulas *)
	d = m*(m*(1 - m)^3*Sin[ϕu]^8 - m*Cos[ϕu]^8 + Δ[ϕu]^8);
	Assert[sn[3*u]/sn[u] == (-m^3*Cos[ϕu]^8 + Δ[ϕu]^8 - (1 - m)^3)/d];
	Assert[cn[3*u]/cn[u] == (1 - m)*(m^3*(1 - m)*Sin[ϕu]^8 + Δ[ϕu]^8 - (1 - m))/d];
	Assert[dn[3*u]/dn[u] == m*(1 - m)*(-m*(1 - m)*Sin[ϕu]^8 + m*Cos[ϕu]^8 + (1 - m))/d];
	Assert[Ε[3*u] == 3*EllipticE[ϕu, m] - m^2*(1 - m)*(Sin[2*ϕu]*Δ[ϕu])^3/d];
	Clear[d];

	(* Half angle formulas *)
	Assert[sn[u/2] == Sin[ϕu/2]*Sqrt[2/(1 + Δ[ϕu])]];
	Assert[cn[u/2] == Cos[ϕu/2]*Sqrt[2*(1 - m)/(1 - m - m*Cos[ϕu] + Δ[ϕu])]];
	Assert[dn[u/2] == kp*Sqrt[(1 + Δ[ϕu])/(1 - m - m*Cos[ϕu] + Δ[ϕu])]];
	Assert[Ε[u/2] == (1/2)*(EllipticE[ϕu, m] + m*Sin[ϕu]*(1 - Cos[ϕu])/(1 + Δ[ϕu]))];

	Assert[sn[(u + v)/2] == Sin[(ϕu + ϕv)/2]*Sqrt[2/(1 + m*Sin[ϕu]*Sin[ϕv] + Δ[ϕu]*Δ[ϕv])]];
	Assert[cn[(u + v)/2] == Cos[(ϕu + ϕv)/2]*Sqrt[2*(1 - m)/(1 - m - m*Cos[ϕu]*Cos[ϕv] + Δ[ϕu]*Δ[ϕv])]];

	(* Quadratic transformations *)
	m1 = 4*k/(1 + k)^2;
	v = (1 + k)*u;
	d = 1 + k*Sin[ϕu]^2;
	Assert[JacobiSN[v, m1] == (1 + k)*Sin[ϕu]/d];
	Assert[JacobiCN[v, m1] == Cos[ϕu]*Δ[ϕu]/d];
	Assert[JacobiDN[v, m1] == (1 - k*Sin[ϕu]^2)/d];
	Clear[d];

	Assert[JacobiSN[v/2, m1] == Sqrt[1 + k]*Sin[ϕu/2]*Sqrt[(1 - m - m*Cos[ϕu] + Δ[ϕu])/((1 + Δ[ϕu])*(Δ[ϕu] - k*Cos[ϕu]))]];
	Assert[JacobiCN[v/2, m1] == Sqrt[1 - k]*Cos[ϕu/2]*Sqrt[(1 + Δ[ϕu])*(Δ[ϕu] + k*Cos[ϕu])/(1 - m - m*Cos[ϕu] + Δ[ϕu])]];
	Assert[JacobiDN[v/2, m1] == (Δ[ϕu] + k*Cos[ϕu])/(1 + k)];

	m1 = ((1 - kp)/(1 + kp))^2;
	v = (1 + kp)*u;
	Assert[JacobiSN[v, m1] == (1 + kp)*Sin[ϕu]*Cos[ϕu]/Δ[ϕu]];
	Assert[JacobiCN[v, m1] == (1 - (1 + kp)*Sin[ϕu]^2)/Δ[ϕu]];
	Assert[JacobiDN[v, m1] == (1 - (1 - kp)*Sin[ϕu]^2)/Δ[ϕu]];

	Assert[JacobiSN[v/2, m1] == (1 + kp)*Sin[ϕu]/(1 + Δ[ϕu])];
	Assert[JacobiCN[v/2, m1] == Cos[ϕu]*Sqrt[2]*Sqrt[1 + kp]/(Sqrt[1 + Δ[ϕu]]*Sqrt[kp + Δ[ϕu]])];
	Assert[JacobiDN[v/2, m1] == Sqrt[2]*Sqrt[kp + Δ[ϕu]]/(Sqrt[1 + kp]*Sqrt[1 + Δ[ϕu]])];

] & // Block[{$AssertFunction = Automatic}, #[]] &;

Integration of elliptic function:

Module[{m, m1, u, n, f},
	m = RandomReal[WorkingPrecision -> 50];
	m1 = 1 - m;
	u = EllipticK[m]*RandomReal[WorkingPrecision -> 50];
	n = RandomInteger[{3, 10}];

	If[Context[JacobiEpsilon] =!= "System`",
		JacobiEpsilon = (EllipticE[JacobiAmplitude[#1, #2], #2]*JacobiDN[#1, #2]/Sqrt[1 - #2*JacobiSN[#1, #2]^2] &);
	];

	Assert[JacobiSN[u, m] == (-Log[JacobiDN[#, m] + Sqrt[m]*JacobiCN[#, m]]/Sqrt[m] &)'[u]];
	Assert[JacobiCN[u, m] == (ArcTan[Sqrt[m]*JacobiSN[#, m]/JacobiDN[#, m]]/Sqrt[m] &)'[u]];
	Assert[JacobiDN[u, m] == (ArcTan[JacobiSN[#, m]/JacobiCN[#, m]] &)'[u]];
	Assert[JacobiCD[u, m] == (Log[JacobiND[#, m] + Sqrt[m]*JacobiSD[#, m]]/Sqrt[m] &)'[u]];
	Assert[JacobiSD[u, m] == (-ArcTan[Sqrt[m]*JacobiCD[#, m]/(Sqrt[m1]*JacobiND[#, m])]/(Sqrt[m]*Sqrt[m1]) &)'[u]];
	Assert[JacobiND[u, m] == (ArcTan[Sqrt[m1]*JacobiSD[#, m]/JacobiCD[#, m]]/(Sqrt[m1]) &)'[u]];
	Assert[JacobiDC[u, m] == (Log[JacobiNC[#, m] + JacobiSC[#, m]] &)'[u]];
	Assert[JacobiNC[u, m] == (Log[JacobiDC[#, m] + Sqrt[m1]*JacobiSC[#, m]]/Sqrt[m1] &)'[u]];
	Assert[JacobiSC[u, m] == (Log[JacobiDC[#, m] + Sqrt[m1]*JacobiNC[#, m]]/Sqrt[m1] &)'[u]];
	Assert[JacobiNS[u, m] == (-Log[JacobiDS[#, m] + JacobiCS[#, m]] &)'[u]];
	Assert[JacobiDS[u, m] == (-Log[JacobiNS[#, m] + JacobiCS[#, m]] &)'[u]];
	Assert[JacobiCS[u, m] == (-Log[JacobiNS[#, m] + JacobiDS[#, m]] &)'[u]];

	Assert[JacobiSN[u, m]^2 == ((-JacobiEpsilon[#, m] + #)/m &)'[u]];
	Assert[JacobiCN[u, m]^2 == ((JacobiEpsilon[#, m] - m1*#)/m &)'[u]];
	Assert[JacobiDN[u, m]^2 == (JacobiEpsilon[#, m] &)'[u]];
	Assert[JacobiCD[u, m]^2 == ((-JacobiEpsilon[#, m] + # + m*JacobiSN[#, m]*JacobiCD[#, m])/m &)'[u]];
	Assert[JacobiSD[u, m]^2 == ((JacobiEpsilon[#, m] - m1*# - m*JacobiSN[#, m]*JacobiCD[#, m])/(m*m1) &)'[u]];
	Assert[JacobiND[u, m]^2 == ((JacobiEpsilon[#, m] - m*JacobiSN[#, m]*JacobiCD[#, m])/m1 &)'[u]];
	Assert[JacobiDC[u, m]^2 == (-JacobiEpsilon[#, m] + # + JacobiSN[#, m]*JacobiDC[#, m] &)'[u]];
	Assert[JacobiNC[u, m]^2 == ((-JacobiEpsilon[#, m] + m1*# + JacobiSN[#, m]*JacobiDC[#, m])/m1 &)'[u]];
	Assert[JacobiSC[u, m]^2 == ((-JacobiEpsilon[#, m] + JacobiSN[#, m]*JacobiDC[#, m])/m1 &)'[u]];
	Assert[JacobiNS[u, m]^2 == (-JacobiEpsilon[#, m] + # - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];
	Assert[JacobiDS[u, m]^2 == (-JacobiEpsilon[#, m] + m1*# - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];
	Assert[JacobiCS[u, m]^2 == (-JacobiEpsilon[#, m] - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];

	Assert[JacobiSN[2*u, m] == (ArcTanh[Sqrt[m]*JacobiSN[#, m]^2]/Sqrt[m] &)'[u]];
	Assert[JacobiCN[2*u, m] == (ArcTan[Sqrt[m]*JacobiCN[#, m]*JacobiSD[#, m]]/Sqrt[m] &)'[u]];
	Assert[JacobiDN[2*u, m] == (ArcTan[JacobiDN[#, m]*JacobiSC[#, m]] &)'[u]];
	Assert[JacobiCD[2*u, m] == (ArcTanh[Sqrt[m]*JacobiCD[#, m]*JacobiSN[#, m]]/Sqrt[m] &)'[u]];
	Assert[JacobiSD[2*u, m] == (ArcTan[Sqrt[m]*Sqrt[m1]*JacobiSD[#, m]^2]/(Sqrt[m]*Sqrt[m1]) &)'[u]];
	Assert[JacobiND[2*u, m] == (ArcTan[Sqrt[m1]*JacobiND[#, m]*JacobiSC[#, m]]/Sqrt[m1] &)'[u]];
	Assert[JacobiDC[2*u, m] == (ArcTanh[JacobiDC[#, m]*JacobiSN[#, m]] &)'[u]];
	Assert[JacobiNC[2*u, m] == (ArcTanh[Sqrt[m1]*JacobiNC[#, m]*JacobiSD[#, m]]/Sqrt[m1] &)'[u]];
	Assert[JacobiSC[2*u, m] == (ArcTanh[Sqrt[m1]*JacobiSC[#, m]^2]/Sqrt[m1] &)'[u]];
	Assert[JacobiNS[2*u, m] == (Log[JacobiNC[#, m]*JacobiSD[#, m]]/2 &)'[u]];
	Assert[JacobiDS[2*u, m] == (Log[JacobiSN[#, m]*JacobiDC[#, m]]/2 &)'[u]];
	Assert[JacobiCS[2*u, m] == (Log[JacobiSN[#, m]*JacobiCD[#, m]]/2 &)'[u]];

	f[u_] := JacobiSN[u, m]^(n - 3)*JacobiCN[u, m]*JacobiDN[u, m];
	f[u_, n_] := JacobiSN[u, m]^n;
	Assert[f[u, n] == 1/(m*(n - 1))*f'[u] + (1 + m)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(m*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiCN[u, m]^(n - 3)*JacobiSN[u, m]*JacobiDN[u, m];
	f[u_, n_] := JacobiCN[u, m]^n;
	Assert[f[u, n] == 1/(m*(n - 1))*f'[u] + (m - m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] + m1*(n - 3)/(m*(n - 1))*f[u, n - 4]];
	Clear[f];
	
	f[u_] := JacobiDN[u, m]^(n - 3)*JacobiSN[u, m]*JacobiCN[u, m];
	f[u_, n_] := JacobiDN[u, m]^n;
	Assert[f[u, n] == m/(n - 1)*f'[u] + (1 + m1)*(n - 2)/(n - 1)*f[u, n - 2] - m1*(n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiCD[u, m]^(n - 3)*JacobiSD[u, m]*JacobiND[u, m];
	f[u_, n_] := JacobiCD[u, m]^n;
	Assert[f[u, n] == -m1/(m*(n - 1))*f'[u] + (1 + m)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(m*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiSD[u, m]^(n - 3)*JacobiCD[u, m]*JacobiND[u, m];
	f[u_, n_] := JacobiSD[u, m]^n;
	Assert[f[u, n] == -1/(m*m1*(n - 1))*f'[u] + (m - m1)*(n-2)/(m*m1*(n-1))*f[u, n - 2] + (n - 3)/(m*m1*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiND[u, m]^(n - 3)*JacobiCD[u, m]*JacobiSD[u, m];
	f[u_, n_] := JacobiND[u, m]^n;
	Assert[f[u, n] == -m/(m1*(n - 1))*f'[u] + (1 + m1)*(n - 2)/(m1*(n - 1))*f[u, n - 2] - (n - 3)/(m1*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiDC[u, m]^(n - 3)*JacobiNC[u, m]*JacobiSC[u, m];
	f[u_, n_] := JacobiDC[u, m]^n;
	Assert[f[u, n] == m1/(n - 1)*f'[u] + (1 + m)*(n - 2)/(n - 1)*f[u, n - 2] - m*(n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiNC[u, m]^(n - 3)*JacobiDC[u, m]*JacobiSC[u, m];
	f[u_, n_] := JacobiNC[u, m]^n;
	Assert[f[u, n] == 1/(m1*(n - 1))*f'[u] + (m1 - m)*(n - 2)/(m1*(n - 1))*f[u, n - 2] + m*(n - 3)/(m1*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiSC[u, m]^(n - 3)*JacobiDC[u, m]*JacobiNC[u, m];
	f[u_, n_] := JacobiSC[u, m]^n;
	Assert[f[u, n] == 1/(m1*(n - 1))*f'[u] - (1 + m1)*(n - 2)/(m1*(n - 1))*f[u, n - 2] - (n - 3)/(m1*(n - 1))*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiNS[u, m]^(n - 3)*JacobiDS[u, m]*JacobiCS[u, m];
	f[u_, n_] := JacobiNS[u, m]^n;
	Assert[f[u, n] == -1/(n - 1)*f'[u] + (1 + m)*(n - 2)/(n - 1)*f[u, n - 2] - m*(n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiDS[u, m]^(n - 3)*JacobiNS[u, m]*JacobiCS[u, m];
	f[u_, n_] := JacobiDS[u, m]^n;
	Assert[f[u, n] == -1/(n - 1)*f'[u] + (m1 - m)*(n - 2)/(n - 1)*f[u, n - 2] + m*m1*(n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := JacobiCS[u, m]^(n - 3)*JacobiNS[u, m]*JacobiDS[u, m];
	f[u_, n_] := JacobiCS[u, m]^n;
	Assert[f[u, n] == -1/(n - 1)*f'[u] - (1 + m1)*(n - 2)/(n - 1)*f[u, n - 2] - m1*(n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := Cos[JacobiAmplitude[u, m]]*JacobiDN[u, m];
	f[u_, n_] := Sin[JacobiAmplitude[u, m]]
	Assert[f[u, n] == -4/(m*(n - 1))*f'[u] - 2*(1 + m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_] := Sin[JacobiAmplitude[u, m]]*JacobiDN[u, m];
	f[u_, n_] := Cos[JacobiAmplitude[u, m]]
	Assert[f[u, n] == 4/(m*(n - 1))*f'[u] - 2*(1 + m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(n - 1)*f[u, n - 4]];
	Clear[f];

	f[u_?InexactNumberQ] := -With[{k = Sqrt[m]}, Sqrt[2/k]*InverseJacobiSD[Sqrt[2*k/m1]*JacobiCN[u/2, m], (1 + k)/2]];
	Assert[f'[u] == Sin[(1/2)*JacobiAmplitude[u, m]]];
	Clear[f];

	f[u_?InexactNumberQ] := With[{k = Sqrt[m]}, Sqrt[2/k]*InverseJacobiSD[Sqrt[2*k]*JacobiSD[u/2, m], (1 + k)/2]];
	Assert[f'[u] == Cos[(1/2)*JacobiAmplitude[u, m]]];
	Clear[f];
] & // Block[{$AssertFunction = Automatic}, #[]] &;

Integration of zeta squared:

If[Context[JacobiZN] =!= "System`",
	JacobiZN[u_, m_] := JacobiZeta[JacobiAmplitude[u, m], m];
];

f[u_, m_, n0_] :=
Module[{pr, q, z, m1, t, n, res = $Failed},
	If[
		VectorQ[{u, m}, NumericQ] && (pr = Accuracy[{u, m}]) < Infinity,
		pr += 5;
		m1 = SetAccuracy[m, pr];
		q = EllipticNomeQ[m1];
		z = Exp[I*Pi*SetAccuracy[u, pr]/EllipticK[m1]];
		res = 1;
		n = n0;
		While[
			t = res;
			res *= ((1 - q^n*z)/(1 - q^n/z))^n;
			t != res,
			n += 2;
		];
		res = I*Log[res];
		res = SetAccuracy[res, Accuracy[res] - 5];
	];
	res /; NumberQ[res]
];

On[Assert];
m = RandomReal[{1/3, 2/3}, WorkingPrecision -> 50];
u = RandomReal[{0, EllipticK[m]}, WorkingPrecision -> 50];
Assert[
	Derivative[1, 0, 0][f][u, m, 1] ==
	(EllipticK[m]/Pi*(JacobiZN[u, m]^2 - m*JacobiSN[u, m]^2 + (1/3)*(1 + m)) - Pi/(12*EllipticK[m]))
];
Assert[
	Derivative[1, 0, 0][f][u, m, 2] ==
	(EllipticK[m]/Pi*(JacobiZN[u, m]^2 + 2*JacobiZN[u, m]*JacobiCN[u, m]*JacobiDS[u, m] + m*JacobiSN[u, m]^2 - (2/3)*(1 + m)) + Pi/(6*EllipticK[m]))
];

Repeated integration:

On[Assert];

f[z_] :=
Module[{n = 1, a = 1, b = 1, c, t, s = 0, r},
	c = 1 - SetPrecision[z, Precision[z] + 5];
	t = c/4;
	While[
		r = s;
		s += a*t;
		r != s,
		b *= n*(n+1);
		n += 2;
		a = a*n^2 + b;
		t *= c/(n + 1)^2;
	];
	SetPrecision[s, Precision[s] - 5]
];

m = RandomReal[{1/3, 2/3}, WorkingPrecision -> 50];
Assert[f[m] == Sqrt[m]*HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, m] - (Pi/2)*EllipticK[m] + (4*Catalan/Pi)*EllipticK[1 - m]];

With[{z = SetPrecision[EllipticK[m], Infinity], m1 = SetPrecision[m, Infinity], kp = Sqrt[1 - m]},
	(* ∫_0^K d(u') ∫_0^u' d(u) sn(u) *)
	Assert[NIntegrate[(z - t)*JacobiSN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == m^(-1/2)*ArcTanh[m^(1/2)]*z - (1 - m)^(-1)*HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, m/(m - 1)]];
	(* ∫_0^K d(u') ∫_0^u' d(u) cn(u) *)
	Assert[NIntegrate[(z - t)*JacobiCN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, m]];
	(* ∫_0^K d(u') ∫_0^u' d(u) dn(u) *)
	Assert[NIntegrate[(z - t)*JacobiDN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == (1/4)*Pi*z + m/(1 + kp)^3*HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, ((1 - kp)/(1 + kp))^2]];
];

References:
Transformations of the Jacobian amplitude function

Last edited by lanxiyu (Yesterday 13:47:05)

Offline

Board footer

Powered by FluxBB