Math Is Fun Forum

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

You are not logged in.

#1 2025-11-02 17:56:32

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

Mathematica test for Elliptic integrals

Duplication theorems:

c[x_, y_] := CarlsonRC[x, y];
d[x_, y_, z_] := CarlsonRD[x, y, z];
f[x_, y_, z_] := CarlsonRF[x, y, z];
g[x_, y_, z_] := CarlsonRG[x, y, z] - (1/6)*(x + y + z)*CarlsonRF[x, y, z];
j[x_, y_, z_, p_] := CarlsonRJ[x, y, z, p];
i[x_, y_, z_] := With[{a = (x + y + z)/3},
	WeierstrassSigma[CarlsonRF[x, y, z], {(2/3)*((x - y)^2 + (x - z)^2 + (y - z)^2), 4*(a - x)*(a - y)*(a - z)}]
];

{x, y, z, p} = RandomReal[{1, 2}, 4, WorkingPrecision -> 50];
λ = Sqrt[x]*Sqrt[y] + Sqrt[x]*Sqrt[z] + Sqrt[y]*Sqrt[z];

On[Assert];
Assert[c[x, y] == c[(Sqrt[x] + Sqrt[y])^2/4, (y + Sqrt[x]*Sqrt[y])/2]];
Assert[d[x, y, z] == (1/4)*d[(x + λ)/4, (y + λ)/4, (z + λ)/4] + 3/(Sqrt[z]*(z + λ))];
Assert[f[x, y, z] == f[(x + λ)/4, (y + λ)/4, (z + λ)/4]];
Assert[g[x, y, z] == 4*g[(x + λ)/4, (y + λ)/4, (z + λ)/4] - (1/2)*(Sqrt[x] + Sqrt[y] + Sqrt[z])];
Assert[i[x, y, z] == (1/8)*i[(x + λ)/4, (y + λ)/4, (z + λ)/4]^4*(Sqrt[x] + Sqrt[y])*(Sqrt[x] + Sqrt[z])*(Sqrt[y] + Sqrt[z])];

α = p*(Sqrt[x] + Sqrt[y] + Sqrt[z]) + Sqrt[x]*Sqrt[y]*Sqrt[z];
β = Sqrt[p]*(p + λ);
Assert[j[x, y, z, p] == (1/4)*j[(x + λ)/4, (y + λ)/4, (z + λ)/4, (p + λ)/4] + 3*CarlsonRC[α^2, β^2]];

δ = (p - x)*(p - y)*(p - z); (* β^2 - α^2 *)
a = Sqrt@Abs[δ];  (* sqrt(|β^2 - α^2|) *)
b = (Sqrt[p] + Sqrt[x])*(Sqrt[p] + Sqrt[y])*(Sqrt[p] + Sqrt[z]); (* β + α *)

(* Circular case *)
If[δ > 0, Assert[CarlsonRC[α^2, β^2] == 2*ArcTan[b, a]/a]];
(* Hyperbolic case *)
If[δ < 0, Assert[CarlsonRC[α^2, β^2] == Log[1 + a*(1 + a/b)/β]/a]];
(* Degenerate case *)
If[δ == 0, Assert[CarlsonRC[α^2, β^2] == 1/α == 1/β]];

Addition theorems:

x = RandomReal[{0, 1}, WorkingPrecision -> 50];
{y, z, p} = x + RandomReal[{0, 1}, 3, WorkingPrecision -> 50];

a = y - x;
b = z - x;
c = p - x;
u = a*b/x;
v = a*z/x;
w = b*y/x;
q = u + c;

On[Assert];

Assert[CarlsonRC[x, y] + CarlsonRC[a^2/x, a*y/x] == Pi/(2*Sqrt[a])];

Assert[CarlsonRF[x, y, z] + CarlsonRF[u, v, w] == CarlsonRF[0, a, b]];

Assert[CarlsonRG[x, y, z] + CarlsonRG[u, v, w] == CarlsonRG[0, a, b]
	+ (1/2)*((x - u)*CarlsonRF[x, y, z] + u*CarlsonRF[0, a, b] + Sqrt[y*z/x])];

Assert[CarlsonRJ[x, y, z, p] + CarlsonRJ[u, v, w, q] == CarlsonRJ[0, a, b, c]
	- 3*CarlsonRC[c^2*y*z/x, c*p*q]];

Assert[CarlsonRD[x, y, z] + CarlsonRD[u, v, w] == CarlsonRD[0, a, b]
	- 3/(b*Sqrt[y*z/x])];

Assert[CarlsonRD[x, z, y] + CarlsonRD[u, w, v] == CarlsonRD[0, b, a]
	- 3/(a*Sqrt[y*z/x])];

Assert[a*b*(CarlsonRD[y, z, x] + CarlsonRD[v, w, u]) == -6*CarlsonRG[0, a, b]
	+ 3*Sqrt[y*z/x]];

Addition theorems (General case):

Module[{x, y, z, p, λ, μ, ν, γ, δ, rfλ, rfμ, rfν, rdλ, rdμ, rdν, rgλ, rgμ, rgν, rjλ, rjμ, rjν},

	{x, y, z} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50];

	{λ, μ} = x + y + z + RandomReal[{0, 1}, 2, WorkingPrecision -> 50];
	ν = ((Sqrt[(x + λ)*(y + λ)*(z + λ)] - Sqrt[(x + μ)*(y + μ)*(z + μ)])/(λ - μ))^2 - λ - μ - x - y - z;

	Assert[Sqrt[z + ν]
		== (Sqrt[(x + λ)*(y + λ)*(z + μ)] - Sqrt[(x + μ)*(y + μ)*(z + λ)])/(λ - μ)
		== ((z + λ)*(z + μ) - (z - x)*(z - y))/(Sqrt[(x + λ)*(y + λ)*(z + μ)] + Sqrt[(x + μ)*(y + μ)*(z + λ)])
	];

	p = Min[x, y, z] + RandomReal[{0, 1}, WorkingPrecision -> 50];

	rfλ = CarlsonRF[x + λ, y + λ, z + λ];
	rfμ = CarlsonRF[x + μ, y + μ, z + μ];
	rfν = CarlsonRF[x + ν, y + ν, z + ν];

	rdλ = CarlsonRD[x + λ, y + λ, z + λ];
	rdμ = CarlsonRD[x + μ, y + μ, z + μ];
	rdν = CarlsonRD[x + ν, y + ν, z + ν];

	rgλ = CarlsonRG[x + λ, y + λ, z + λ];
	rgμ = CarlsonRG[x + μ, y + μ, z + μ];
	rgν = CarlsonRG[x + ν, y + ν, z + ν];

	rjλ = CarlsonRJ[x + λ, y + λ, z + λ, p + λ];
	rjμ = CarlsonRJ[x + μ, y + μ, z + μ, p + μ];
	rjν = CarlsonRJ[x + ν, y + ν, z + ν, p + ν];

	Assert[rfλ + rfμ == rfν];

	Assert[rdλ + rdμ == rdν - 3/Sqrt[(z + λ)*(z + μ)*(z + ν)]];

	Assert[rgλ + rgμ == rgν + (1/2)*(λ*rfλ + μ*rfμ - ν*rfν + Sqrt[λ + μ + ν + x + y + z])];

	γ = (p + λ)*(p + μ)*(p + ν);
	δ = (p - x)*(p - y)*(p - z);
	Assert[rjλ + rjμ == rjν - 3*CarlsonRC[γ - δ, γ]];

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

Carlson form ↔ Legendre form:

Module[{x, y, z, p, Δ, ϕ, m, m1, n, pr},
	Δ = Sqrt[1 - m*Sin[#]^2] &;

	pr = 50;
	ϕ = RandomReal[{0, Pi/2}, WorkingPrecision -> pr];
	m = RandomReal[{0, 1}, WorkingPrecision -> pr];
	m1 = 1 - m;
	n = 1 - 1/RandomReal[{0, 1}, WorkingPrecision -> pr];

	Assert[EllipticF[ϕ, m] ==
		+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
	];
	Assert[EllipticE[ϕ, m] ==
		+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
		- (1/3)*m*Sin[ϕ]^3*CarlsonRD[Cos[ϕ]^2, Δ[ϕ]^2, 1]
	];
	Assert[EllipticPi[n, ϕ, m] ==
		+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
		+ (1/3)*n*Sin[ϕ]^3*CarlsonRJ[Cos[ϕ]^2, Δ[ϕ]^2, 1, 1 - n*Sin[ϕ]^2]
	];
	Assert[JacobiZeta[ϕ, m] ==
		+ (1/3)*m*Sin[ϕ]*Cos[ϕ]*Δ[ϕ]*CarlsonRJ[0, m1, 1, Δ[ϕ]^2]/EllipticK[m]
	];

	{x, y, z, p} = {Cos[ϕ]^2, Δ[ϕ]^2, 1, 1 - n*Sin[ϕ]^2}*
		RandomReal[{0, 1}, WorkingPrecision -> pr];

	Assert[ϕ == ArcCos@Sqrt[x/z]];
	Assert[m == (z - y)/(z - x)];
	Assert[m1 == (y - x)/(z - x)];
	Assert[n == (z - p)/(z - x)];
	Assert[m - n == (p - y)/(z - x)];
	Assert[1 - n == (p - x)/(z - x)];

	Assert[CarlsonRC[x, z] ==
		(z - x)^(-1/2)*ϕ];
	Assert[CarlsonRF[x, y, z] ==
		(z - x)^(-1/2)*EllipticF[ϕ, m]];
	Assert[CarlsonRD[x, y, z] ==
		(3/m)*(z - x)^(-3/2)*(EllipticF[ϕ, m] - EllipticE[ϕ, m])];
	Assert[CarlsonRJ[x, y, z, p] ==
		(3/n)*(z - x)^(-3/2)*(EllipticPi[n, ϕ, m] - EllipticF[ϕ, m])];
	Assert[2*CarlsonRG[x, y, z] ==
		(z - x)^(1/2)*(EllipticE[ϕ, m] + Cot[ϕ]^2*EllipticF[ϕ, m] + Cot[ϕ]*Δ[ϕ])];
] & // Block[{$AssertFunction = Automatic}, #[]] &;

Rhombic lattice → Rectangular lattice:

Module[{a0, a1, a2, ap, aq, b0, b1, b2, bp, bq, Δ, ϕ, ψ, m, m1, n, pr, rc, rd, rf, rg, rj, implies},
	{rc, rd, rf, rg, rj} = Function[{f}, Module[{args, res},
		args = {##};
		res = f @@ SetAccuracy[args, Accuracy@args + 5];
		SetAccuracy[res, Accuracy@res - 5]
	] &] /@ {CarlsonRC, CarlsonRD, CarlsonRF, CarlsonRG, CarlsonRJ};

	implies = Function[Null, !TrueQ[#1] || #2, HoldAll];

	Δ = Sqrt[1 - m*Sin[#]^2] &;

	pr = 50;
	a0 = RandomReal[{0, 1}, WorkingPrecision -> pr];
	a1 = RandomComplex[{0, 1 + I}, WorkingPrecision -> pr];
	a2 = Conjugate[a1];
	ap = RandomReal[{0, 1}, WorkingPrecision -> pr];
	aq = a0 + Abs[a0 - a1]^2/(ap - a0);

	b0 = (1/2)*(Abs[a1] + Re[a1]);
	b1 = (1/2)*(Abs[a1] + Abs[a0 - a1] + a0);
	b2 = (1/2)*(Abs[a1] - Abs[a0 - a1] + a0);
	bp = (1/2)*(Abs[a1] + Abs[ap - a1] + ap);
	bq = (1/2)*(Abs[a1] - Abs[ap - a1] + ap);

	ϕ = 2*ArcTan[Abs[Sqrt[a0] + Sqrt[a1]], Sqrt@Abs[a0 - a1]];
	ψ = 2*ArcTan[Sqrt[a0], Sqrt@Abs[a0 - a1]];
	{m, m1} = ReIm@Sqrt[a0 - a1]^2/Abs[a0 - a1];
	n = (b1 - bp)/Abs[a0 - a1];

	Assert[0 <= b2 <= b0 <= b1];
	Assert[0 <= bq <= b0 <= bp];

	Assert[Sqrt[b0] == Re@Sqrt[a1]];
	Assert[Sqrt[b1] == (1/2)*(Abs[Sqrt[a0] + Sqrt[a1]] + Abs[Sqrt[a0] - Sqrt[a1]])];
	Assert[Sqrt[b2] == (1/2)*(Abs[Sqrt[a0] + Sqrt[a1]] - Abs[Sqrt[a0] - Sqrt[a1]])];
	Assert[Sqrt[bp] == (1/2)*(Abs[Sqrt[ap] + Sqrt[a1]] + Abs[Sqrt[ap] - Sqrt[a1]])];
	Assert[Sqrt[bq] == (1/2)*(Abs[Sqrt[ap] + Sqrt[a1]] - Abs[Sqrt[ap] - Sqrt[a1]])];

	Assert[Sqrt[b1] + Sqrt[b2] == Abs[Sqrt[a0] + Sqrt[a1]]];
	Assert[Sqrt[b1] - Sqrt[b2] == Abs[Sqrt[a0] - Sqrt[a1]]];

	Assert[(bp - b0)*(b0 - bq) == (b1 - b0)*(b0 - b2)];

	Assert[a0 == b1*b2/b0];
	Assert[ap == bp*bq/b0];
	Assert[Re[a1] == 2*b0 - Abs[a1]];
	Assert[Abs@Im[a1] == 2*Sqrt[b1 - b0]*Sqrt[b0 - b2]];
	Assert[Abs[a1] == b2 + b1 - a0];

	Assert[Abs[a0 - a1] == b1 - b2];
	Assert[Re@Sqrt[a0 - a1]^2 == b1 - b0];
	Assert[Im@Sqrt[a0 - a1]^2 == b0 - b2];

	Assert[Abs[ap - a1] == bp - bq];
	Assert[Re@Sqrt[ap - a1]^2 == bp - b0];
	Assert[Im@Sqrt[ap - a1]^2 == b0 - bq];

	Assert[Sin[ϕ] == Sqrt[Abs[a0 - a1]/b1]];
	Assert[Cos[ϕ] == Sqrt[b2/b1]];
	Assert[Δ[ϕ] == Sqrt[b0/b1]];

	Assert[Cos[ϕ]/Δ[ϕ] == Sqrt[a0/b1]];
	Assert[Tan[ϕ]*Δ[ϕ] == Tan[ψ/2] == Sqrt[Abs[a0 - a1]/a0]];
	Assert[Sin[ψ] == 2*Sqrt[a0*Abs[a0 - a1]]/(a0 + Abs[a0 - a1])];
	Assert[Cos[ψ] == (a0 - Abs[a0 - a1])/(a0 + Abs[a0 - a1])];
	Assert[Δ[ψ] == Abs[a1]/(a0 + Abs[a0 - a1])];

	Assert[1 - m*(Sin[ϕ]*Cos[ϕ]/Δ[ϕ])^2 == Abs[a1]/b1];
	Assert[Csc[ϕ]^2 - m*(Cos[ϕ]/Δ[ϕ])^2 == Abs[a1]/Abs[a0 - a1]];

	Assert[Max[0, n] <= m <= 1];
	Assert[m - n == (bp - b0)/Abs[a0 - a1]];
	Assert[1 - n == (bp - b2)/Abs[a0 - a1]];
	Assert[m*m1/(m - n) == (b0 - bq)/Abs[a0 - a1]];
	Assert[m*(1 - n)/(m - n) == (b1 - bq)/Abs[a0 - a1]];
	Assert[m1*n/(m - n) == (b2 - bq)/Abs[a0 - a1]];

	(* CarlsonRF *)
	Assert[rf[a1, a2, a0]
		== rf[b1, b2, b0]
		== Abs[a0 - a1]^(-1/2)*EllipticF[ϕ, m]
		== (1/2)*Abs[a0 - a1]^(-1/2)*EllipticF[ψ, m]
	];

	(* CarlsonRD *)
	Assert[Abs[a0 - a1]*rd[a1, a2, a0]
		== 2*(b1 - b0)*rd[b0, b2, b1] - 3*rf[b1, b2, b0] + 3/Sqrt[a0]
		== 2*(b0 - b2)*rd[b0, b1, b2] + 3*rf[b1, b2, b0] - 3/Sqrt[a0]
		== 3*Abs[a0 - a1]^(-1/2)*(
			- 2*EllipticE[ϕ, m]
			+ EllipticF[ϕ, m]
			+ Tan[ϕ]*Δ[ϕ]
		)
		== 3*Abs[a0 - a1]^(-1/2)*(
			- EllipticE[ψ, m]
			+ (1/2)*EllipticF[ψ, m]
			+ Δ[ψ]*Tan[ψ/2]
		)
	];

	(* CarlsonRG *)
	Assert[rg[a1, a2, a0]
		== 2*rg[b1, b2, b0] - (1/2)*(Abs[a1]*rf[b1, b2, b0] + Sqrt[a0])
		== (1/12)*Im[a1]^2*rd[b1, b2, b0] + (1/2)*(Re[a1]*rf[b1, b2, b0] + Sqrt[a0])
		== (1/2)*Abs[a0 - a1]^(1/2)*(
			+ 2*EllipticE[ϕ, m]
			+ (Cot[ϕ]^2 - m1/Δ[ϕ]^2)*EllipticF[ϕ, m]
			+ Cot[ϕ]*(1 - 2*m*Sin[ϕ]^2)/Δ[ϕ]
		)
		== (1/2)*Abs[a0 - a1]^(1/2)*(
			+ EllipticE[ψ, m]
			+ Cot[ψ/2]*(Cot[ψ]*EllipticF[ψ, m] + Δ[ψ])
		)
	];

	(* CarlsonRJ *)
	Assert[Abs[ap - a1]*rj[a1, a2, a0, ap]
		== 2*(bp - b0)*rj[b1, b2, b0, bp] - 3*rf[b1, b2, b0] + 3*rc[a0, ap]
		== 2*(b0 - bq)*rj[b1, b2, b0, bq] + 3*rf[b1, b2, b0] - 3*rc[a0, ap]
		== 3*Abs[a0 - a1]^(-1/2)*(
			+ 2*(m/n - 1)*EllipticPi[n, ϕ, m]
			+ (1 - 2*m/n)*EllipticF[ϕ, m]
		) + 3*rc[a0, ap]
	];

	Assert@Implies[0 < aq, With[{γ = n*(1 - n), δ = (m - n^2)},
		(ap - a0)*rj[a0, a1, a2, ap]
		== (3/2)*Abs[a0 - a1]^(-1/2)*(
			+ (1/2 - γ/δ)*EllipticPi[δ^2/(4*γ*(m - n)), ψ, m]
			+ (γ/δ)*EllipticF[ψ, m]
		) - (3/2)*rc[Abs[a1]^2/a0, ap*aq/a0]
	]];

	(* Branch cut of CarlsonRF *)
	Assert@implies[a0 <= 0,
		Re@rf[a1, a2, a0] == rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]] &&
		Abs@Im@rf[a1, a2, a0] == Sqrt[-a0]*rf[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
	];

	(* Branch cut of CarlsonRD *)
	Assert@implies[a0 < 0,
		Re@rd[a1, a2, a0] == -Abs[a1 - a0]^(-1)*(
			+ 3*rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
			- 2*Im@Sqrt[a1 - a0]^2*rd[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
		) &&
		Abs@Im@rd[a1, a2, a0] == Abs[
			+ 3*Abs[a1]/(Sqrt[-a0]*Abs[a1 - a0]^2)
			- (-a0)^(3/2)*rd[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
		]
	];

	(* Branch cut of CarlsonRG *)
	Assert@implies[a0 <= 0,
		Re@rg[a1, a2, a0] == (
			+ (1/2)*(a0 + Abs[a1 - a0])*rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
			- (1/3)*Abs[a1 - a0]*Im@Sqrt[a1 - a0]^2*rd[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
		) &&
		Abs@Im@rg[a1, a2, a0] == (1/2)*(-a0)^(3/2)*Abs[
			+ rf[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
			- (1/3)*Abs[a1 - a0]^2*rd[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
		]
	];

	(* Inequalities *)
	Assert@implies[a0 >= 0 && Re[a1] > 0,
		rc[a0, Abs[a1]] <= Re@rf[a1, a2, a0] <= rc[a0, Re[a1]] &&
		implies[0 < a0 < Re[a1],
			rc[a0, Re[a1]] == Abs[a0 - a1]^(-1/2)*(m1 - m)^(-1/2)*ArcTan[Sqrt[m1 - m]*Tan[ϕ]*Δ[ϕ]]] &&
		implies[a0 > Re[a1],
			rc[a0, Re[a1]] == Abs[a0 - a1]^(-1/2)*(m - m1)^(-1/2)*ArcTanh[Sqrt[m - m1]*Tan[ϕ]*Δ[ϕ]]]
	];

	Assert@implies[0 < a0 < Re[a1],
		Re@rf[a1, a2, a0] >= (2/Pi)*Abs[a0 - a1]^(-1/2)*EllipticK[m]*ArcCos@Sqrt[a0/Re[a1]]
	];

	Assert@implies[a0 >= Re[a1],
		Re@rf[a1, a2, a0] >= a0^(-1/2)*(1 - (Im[a1]/a0)^2/10)
	];

	Assert@implies[0 < a0 <= Re[a1],
		Re@rf[a1, a2, a0] <= a0^(-1/2)*(1 - Arg[a1]^2/10)
	];

	Assert[Re@rf[a1, a2, a0] < (Log[Pi/Abs@Arg[-a1]] + (1/2)*Log[E^Pi + 4*a0/Abs[a1]])/Sqrt@Abs[a0 - a1]];

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

AGM method:

Module[{pr, ϕ, m, f},

	pr = 50;
	ϕ = RandomReal[{0, Pi/2}, WorkingPrecision -> pr];
	m = RandomReal[{0, 1}, WorkingPrecision -> pr];

	f = Module[{ϕ, m, m1, n, a, b, a1, b1, r, t},
		ϕ = #1;
		m = #2;
		m1 = 1 - m;
		n = Round[Re[ϕ]/Pi];
		ϕ -= n*Pi;
		t = Tan[ϕ/2];
		r = t^2;
		m1 = 1 - m;
		a = 1 + (m1 - m)*r;
		b = Sqrt[a^2 + 4*m*m1*r^2];
		If[ReIm[a].ReIm[b] < 0, b = -b];
		r = Sqrt[1 + 2*m*r/(a + b)];
		a = Sqrt[m1];
		b = 1;
		While[
			a != b || r != 1,
			t *= r;
			a1 = a;
			b1 = b;
			a = (a1 + b1)/2;
			b = Sqrt[a1*b1];
			r = Sqrt[(b1 + a1*r)/(a*(1 + r))];
		];
		Return[(2*ArcTan[a*t] + n*Pi)/a, Module];
	] &;

	Assert[EllipticF[ϕ, m] == f[ϕ, m]];

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

Last edited by lanxiyu (Yesterday 18:06:41)

Offline

Board footer

Powered by FluxBB