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 theorem:

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, 2*(y + Sqrt[x]*Sqrt[y])/4]];
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/β]];

Rhombic lattice -> Rectangular lattice:

Module[{a0, a1, a2, ap, 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];

	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]];
	{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[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[ϕ]*Δ[ϕ] == Sqrt[Abs[a0 - a1]/a0]];
	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[1 - n == (bp - b2)/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]
	];

	(* 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])
		== Abs[a0 - a1]^(1/2)*(
			+ EllipticE[ϕ, m]
			+ (1/2)*(Cot[ϕ]^2 - m1/Δ[ϕ]^2)*EllipticF[ϕ, m]
			+ (1/2)*Cot[ϕ]*(1 - 2*m*Sin[ϕ]^2)/Δ[ϕ]
		)
	];

	(* 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[ϕ]*Δ[ϕ]
		)
	];

	(* 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]
	];

	(* Inequalities *)
	Assert@implies[Re[a1] > 0,
		rc[a0, Abs[a1]] <= Re@rf[a1, a2, a0] <= rc[a0, Re[a1]]];

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

Last edited by lanxiyu (Today 11:50:03)

Offline

Board footer

Powered by FluxBB