<<RiemannHilbert`;
eps=$MachineEpsilon;

This notebook discusses the compution of solutions to the Airy equation

	Y''(z) = z Y(z)
	
throughout the complex plane, as outlined in Sections 1.3 and 6.3.

Setup and solve the RH problem

Here we setup the Riemann–Hilbert Problem that encodes the Airy function in its solution.  While we aim to solve Problem 1.3.1, the singularity at the origin must be removed using the parametrix H(z) (6.6).  The following constructs this parametrix as H[z]. We also use H[z, s] where the parameter s denotes which limit is taken along the branch cut:

H[z_]:=1/2(({{1, (1-1/Sqrt[z])}, {1, (1+1/Sqrt[z])}}));
H[z_,s_]:=1/2(({{1, (1-1/Sqrt[z+s I eps])}, {1, (1+1/Sqrt[z+s I eps])}}));
Hi[z_,s_]:=H[z,s]//Inverse;

The following are the jump matrices corresponding to the local behaviour of H(z):

(Subscript[J, 1])=(({{1, 0 }, {-I, 1}}));
(Subscript[J, 2])=(({{1, -I }, {0, 1}}));
(Subscript[J, 0])=(({{0, -I }, {-I, 0}}));
(*Check product condition*)
Inverse[(Subscript[J, 0])].(Subscript[J, 1]).(Subscript[J, 2]).(Subscript[J, 1])

We now define R[z,s] corresponding to R(z):

R[z_,s_]:=2Sqrt[Pi](z+s I eps)^(1/4)(({{Exp[2/3 (z+s I eps)^(3/2)], 0 }, {0, Exp[-2/3 (z+s I eps)^(3/2)]}})).H[z,s];
R[z_]:=2Sqrt[Pi](z)^(1/4)(({{Exp[2/3 (z)^(3/2)], 0 }, {0, Exp[-2/3 (z)^(3/2)]}})).H[z];
Ri[z_,s_]:=1/(2Sqrt[Pi])(z+s I eps)^(-1/4)Hi[z,s].(({{Exp[-2/3 (z+s I eps)^(3/2)], 0 }, {0, Exp[2/3 (z+s I eps)^(3/2)]}}));
inner[z_]:=({{\[Piecewise], ({{Inverse[(Subscript[J, 1])], 0<Arg[z]≤2Pi/3}, {IdentityMatrix[2], 2Pi/3<Arg[z]≤Pi}, {Inverse[(Subscript[J, 1]).(Subscript[J, 2]).(Subscript[J, 1])], -Pi<Arg[z]<-2Pi/3}, {Inverse[(Subscript[J, 2]).(Subscript[J, 1])], -2Pi/3≤Arg[z]≤0}})}})

Thus we can construct the jumps for the RH Problem depicted in Figure 6.3:

(*Define Jump matrices*)
ry1[z_,s_]:=Ri[z,s].(Subscript[J, 1]);
ry2[z_,s_]:=Ri[z,s].(Subscript[J, 2]);
ryI[z_,s_]:=Ri[z,s];
ry21[z_,s_]:=Ri[z,s].(Subscript[J, 2]).(Subscript[J, 1]);
ry121[z_,s_]:=Ri[z,s].(Subscript[J, 1]).(Subscript[J, 2]).(Subscript[J, 1]);
ry[z_,s_]:=If[Im[z]>0,ryI[z,s],ry121[z,s]];
ry1r[z_,s_]:=Ri[z,s].(Subscript[J, 1]).R[z,s];
ry2r[z_,s_]:=Ri[z,s].(Subscript[J, 2]).R[z,s];

We now approximate the jump functions at n mapped Chebyshev points on each component of the domain:

(* Setup parameters: r is the size of the circle around the origin, L is the outer truncation radius and n is the number of collocation points per contour *)
r=1.5;
L=10;
n=42;
(* Define the RHP *)
RHP={Fun[ry2r[#,1]&,Line[{r,L}],n],
Fun[ry1[#,1]&,Arc[0,r,{0,2Pi/3}],n],
Fun[ry[#,1]&,Arc[0,r,{2Pi/3,4Pi/3}],n],
Fun[ry21[#,1]&,Arc[0,r,{4Pi/3,2Pi}],n],
Fun[ry1r[#,1]&,Line[{r,L}Exp[2Pi I/3]],n],
Fun[ry1r[#,1]&,Line[{r,L}Exp[4Pi I/3]],n]};
P[z_]:=If[Abs[z]≤r,inner[z],Ri[z,1]];

Here is a depiction of the domain that the jump function lives on.

DomainPlot[RHP]

We can now approximate the solution to the RH Problem by finding sol, whose Cauchy transform incodes the solution.  RHSolve sets up and solve the corresponding collocation system.

sol=RHP//RHSolve;

RHSolve gives a matrix-valued approximation to the RH Problem that behaves like the identity at ∞:

Φ[z_]:=IdentityMatrix[2]+Cauchy[sol,z];
Norm[Φ[3.+eps I] - Φ[3.-eps I].RHP[[1]][3.]]

The solution we want behaves like (1,1) at infinity, so we multiply by {1,1}:

ϕ[z_]:=ϕ[z]={1,1}.Φ[z];
{Norm[ϕ[3.+eps I] - ϕ[3.-eps I].RHP[[1]][3.]],ϕ[100.+100.I]}

From ϕ we can recover the solution to the Airy equation as y:

ω=Exp[-2 Pi I/3];
y[z_]:=ϕ[z].P[z];
yexact[z_]:={AiryAi[z],I ω AiryAi[ω z]};

y[2.+2. I]-yexact[2.+2.I]

Error analysis

We can compare our numerical approximation with the exact solution the RH Problem and the asymptotic solution:

ϕexact[z_]=yexact[z].Rsym[z,1]//Simplify;
(* Find a symbolic asymptotic expansion*)
ϕexactseries[z_]=Series[ϕexact[z],{z,∞,6}]//Normal//Simplify;
error[z_]:=Norm[ϕ[z]-ϕexactseries[z]]

TableLogPlot[error[x],{x,1000,5000000,1000},PlotRange->All,PlotStyle->Black,Axes->False,Frame->True,FrameStyle->Directive[FontSize->14]]

TableLogPlot[error[x],{x,300,5000,10},PlotRange->{{300,5000},{10^(-16),10^(-14)}},PlotStyle->Black,Axes->False,Frame->True,FrameStyle->Directive[FontSize->14],FrameLabel->{Style[x,17],Style["Absolute Error",17]}]

Other plotting

(*Airy in the complex plane, we use eps here a lot to make sure we do not accidentally evaluate a limit from the wrong side of the contour*)
mAi[z_]:=({{\[Piecewise], ({{y[z+I eps][[1]], -2Pi/3≤Arg[z]≤2Pi/3}, {y[z+(I-1) eps][[1]]+I y[z+(I-1) eps][[2]], 2Pi/3<Arg[z]≤Pi}, {y[z-(I+1) eps][[1]]-I y[z-(I+1) eps][[2]], -Pi<Arg[z]<-2Pi/3}})}})
AiryAiRH[z_]:=If[Abs[z]<1,mAi[(1-100eps)z],mAi[(1+100eps)z]];

mStyle[x_]:=Style[x,14,Bold];
p1=TablePlot3D[AiryAiRH[x + I Y]//Re,{x,-6,2,.2},{Y,-2,2,.1},AxesLabel->mStyle/@{Re[z],Im[z],Re[Ai[z]]}]
p2=TablePlot3D[AiryAiRH[x + I Y]//Im,{x,-6,2,.2},{Y,-2,2,.1},AxesLabel->mStyle/@{Re[z],Im[z],Im[Ai[z]]}]