<<RiemannHilbert`;

This notebook discusses the computation of the elliptic integral

g(z;a) = (∫_0^z)(1/((Sqrt[1-(s^2)])(Sqrt[(a^2)-(s^2)])))ds

throughout the complex plane, as outlined in Sections 1.2 and 5.5.

Compute g

The computation of g(z;a) is accomplished by reducing the problem to calculating two Cauchy transforms of a function with square root singularities on [-1,1].  This is possible by noting that g(z;a) has two branch cuts: One on [1,a] and the other on [-a,-1].  Thus we conformally map these two intervals, independently, to [-1,1] using M[x]:

a=1.5;
  (* M[x] maps [-1,1] to [a,1] *)
M[x_]:=(1+a+x-a x)/2;  
(*Mi is the inverse of M, hence maps [a,1] to [-1,1] *)
Mi[z_]:=((1+a-2 z)/(-1+a));

We can factor out the square root singularities that arise, and expand the remaining smooth part of the integrand into a 20-term Chebyshev expansion:

f1=Fun[-2I/(Sqrt[1+M[#]]Sqrt[a+M[#]])&,Line[{-1,1}],20];
(*Chebyshev coefficients for non-singular part*)
coefs=f1//DCT;
(*Convert to vanishing basis*)

To obtain simple, closed formulas we convert the expansion from the Chebyshev basis (Subscript[T, k])(x) to the vanishing basis (Subscript[T, 0]^z)(x)  = 1, (Subscript[T, 1]^z)(x)  = x, (Subscript[T, k]^z)(x) = (Subscript[T, k])(x)-(Subscript[T, k-2])(x) which vanishes at ±1 for k≥2, as described in Section 4.4:

ToVanishingCoefsMat[n_]:=Table[Table[If[Mod[i+k,2]==0&& i≥k,1,0],{i,1,n}],{k,1,n}];
(* Vanishing basis coefficients *)
coefs=ToVanishingCoefsMat[Length[coefs]].coefs;

We can now use Lemma 5.10 which gives an explicit formula for the Cauchy transform of (Subscript[T, k]^z) in terms of one of the inverses of the Joukowsky map ((z+(z^(-1)))/2), which we denote Jpi:

(* Define the inverse of (Subscript[J, +]) *)
Jpi[x_]:=x -(Sqrt[x-1]) (Sqrt[1+x]);
(* Map back to the slit plane on [1,a] *)
fp[z_]:=Jpi[Mi[z]];
(* Map back to the slit plane on [-a,-1] *)
fm[z_]:=Jpi[Mi[-z]];

The following commands use the explicit formula for the Cauchy transform, along with the indefinite integration formulas 5.4.4:

h[z_]:=coefs[[1]](Log[fp[z]]-Log[fm[z]])+coefs[[2]](fp[z]-fm[z])+coefs[[3]](fp[z]^2/2-Log[fp[z]]-fm[z]^2/2+Log[fm[z]])+Sum[coefs[[i]](((fp[z]^(i-1)-fm[z]^(i-1))/((i-1)))-((fp[z]^(i-3)-fm[z]^(i-3))/((i-3)))),{i,4,Length[coefs]}];
(* gg[z] Gives g(z; a)*)
gg[z_]:=I/2 h[z];

We can compare the results with NIntegrate:

gN[a_,x_]:=NIntegrate[(1/(Sqrt[1-(z^2)]Sqrt[(a^2)-(z^2)])),{z,0,x}]

(* Check at a given value *)
gg[5I+5]
gN[a,5I+5]

The benefit of gg over gN is that it remains accurate within machine precision up to the branch cut:

gN[a,1.2+0.0001 I]

gg[1.2+0.0001 I]

Mathematica requires high precision arithmetic to achieve the same results:

NIntegrate[(1/(Sqrt[1-(z^2)]Sqrt[(a^2)-(z^2)])),{z,0,1.2+0.0001 I},WorkingPrecision->20,MaxRecursion->20 ]

Note that gg is 300x faster:

NIntegrate[(1/(Sqrt[1-(z^2)]Sqrt[(a^2)-(z^2)])),{z,0,1.2+0.0001 I},WorkingPrecision->20,MaxRecursion->20 ]//Timing

gg[1.2+0.0001 I]//Timing

Plotting

Here we use gg to plot the elliptic integral, which takes advantage of the accuracy up to the branch cut:

mStyle[x_]:=Style[x,14,Bold];
p1=TablePlot3D[gg[x + I y]//Re,{x,-2,2,.1/2},{y,-2,2,.1/4},AxesLabel->mStyle/@{Re[z],Im[z],Re[g[z]]}]
p2=TablePlot3D[gg[x + I y]//Im,{x,-2,2,.1/2},{y,-2,2,.1/4},AxesLabel->mStyle/@{Re[z],Im[z],Im[g[z]]}]