<<RiemannHilbert`;

Calculating Erfc using RHPackage

The error function can be expressed in terms of the Cauchy transform of the Gaussian:

 	erf(z)=({{{, ({{2-((e^(-(z^2)))/(2 π i)) (∫_(-i ∞)^(i ∞))((2(e^(s^2)))/(s-z))ds, Re[z] < 0}, {-((e^(-(z^2)))/(2 π i)) (∫_(-i ∞)^(i ∞))((2(e^(s^2)))/(s-z))ds, Re[z] > 0}})}})= ({{{, ({{2-((e^(-(z^2)))/(2 π i)) (∫_(-∞)^∞)((2(e^(-(x^2))))/(x+i z))dx, Re[z] < 0}, {-((e^(-(z^2)))/(2 π i)) (∫_(-∞)^∞)((2(e^(-(x^2))))/(x+i z))dx, Re[z] > 0}})}})
 
  We can represent the Gaussian on the real line in RHPackage as follows:

f=LFun[Exp[(((I #))^2)]&,RealLine]

We can then define our own error function using the Cauchy function:

erfc[z_]/;Re[z]>0:=-Exp[-(z^2)]Cauchy[2f,-I z];
erfc[z_]/;Re[z]<0:=2-Exp[-(z^2)]Cauchy[2f,-I z];

This is accurate uniformly in the complex plane, which we can check at two points:

erfc[1.+I]-Erfc[1.+I]

erfc[-1.+I]-Erfc[-1.+I]

Here is a plot of the real and imaginary parts of the approximation:

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

Code to compute Erfc

We now dig into the details of what is happening behind the scenes.  The command LFun constructs a rational approximation of the Gaussian.  This is accomplished by conformally mapping to the unit circle and expanding in Laurent series, which in turn is accomplished by mapping to the periodic interval and using Fourier series.  

The following code can be used to construct an approximation of a 2π-periodic function:

(* Evenly spaced points on the periodic interval *)
Uni[n_]:=2Pi Table[i,{i,0,n-1,1}]/(n)-Pi;
(* The correct normalization for the FFT *)
mFourier[x_]:=Module[{a},
a =Fourier[x];
Table[a[[i]](-1)^(i+1),{i,1,Length[a]}]/Sqrt[Length[a]]
];
(*Negative coefficients are at the end, positive coefficients at the frunt *)
IS[x_]:=If[x≥0,x+1,x];

Fapprox evaluates a finite Fourier series at a point θ:

Fapprox[θ_,data_]:=Sum[data[[k//IS]] Exp[-I k θ],{k,-Floor[1/2(Length[data]-1)],Floor[1/2(Length[data])]}]

We now construct the Mobius transformation Μ that maps (-i∞, i∞) to (-π,π] and its invere T:

(* Boolean to help with large values *)
LargeQ[z_]:=Abs[z]>40;
LargeQ[ComplexInfinity]:=True;
M[z_]:=-I (z-1)/(1+z); (* map points on (-i∞,i∞) to (-π,π] *)
iM[z_]:=(I-z)/(I+z);
T[z_]:=I M[Exp[I z]]; (* map points on (-π,π] to (-i∞,i∞)*)

We can use this map to construct a rational approximation of a Gaussian on the imaginary axis:

f[z_]:=2 Exp[z^2];
(* So f[∞] makes sense *)
f[z_?LargeQ]=0;
f[z_List]:=f/@z;

F[z_]:=f[T[z]];
fapprox[x_,data_]:=Sum[data[[k//IS]]iM[-I x]^(-k),{k,-Floor[1/2(Length[data]-1)],Floor[1/2(Length[data])]}]

n=2^8;
(*Compute coefficients via FFT*)
coefs=F[Uni[n]]//mFourier//Quiet;

The functions Cp and Cm give the Cauchy transform evaluated at x:

(*Define the Cauchy operators on the coefficients*)
Cp[x_,data_]:=Sum[data[[k//IS]](iM[-I x]^(-k)-(-1)^(-k)),{k,-Floor[1/2(Length[data]-1)],-1}]
Cm[x_,data_]:=-Sum[data[[k//IS]](iM[-I x]^(-k)-(-1)^(-k)),{k,1,Floor[1/2(Length[data])]}]

The following demonstrates that Plemelj’s lemma holds to machine accuracy:

(* Check Plemelj, stability and analyticity *)
{DiscretePlot[Cp[I x,coefs]-Cm[I x,coefs]-f[I x]//Abs,{x,-5,5,.01}],DiscretePlot[Cp[ x,coefs]//Re,{x,-2,0,.01}],DiscretePlot[Cp[ x Exp[3 Pi I/4],coefs]//Im,{x,0,2,.01}]}

We can thus create an approximation erfc to Erfc from the rational approximation to a Gaussian:

ψ[z_]:=If[Re[z]>0,Cm[z,coefs],Cp[z,coefs]];
(* Exact formula for ψ *)
ψexact[z_]:=If[Re[z]>0,-Exp[z^2]Erfc[z],2Exp[z^2]-Exp[z^2]Erfc[z]];
(* Asymptotic expansion for ψ *)
ψexactseries[z_]:=If[Re[z]>0,-(1/((Sqrt[π]) z))+(1/(2 (Sqrt[π]) (z^3)))-(3/(4 (Sqrt[π]) (z^5)))+(15/(8 (Sqrt[π]) (z^7))),-(1/((Sqrt[π]) z))+(1/(2 (Sqrt[π]) (z^3)))-(3/(4 (Sqrt[π]) (z^5)))+(15/(8 (Sqrt[π]) (z^7)))];
ψexact[s_,z_]:=N[ψexact[SetPrecision[z-s $MachineEpsilon,100]],100];
ψ[s_,z_]:=ψ[z-s $MachineEpsilon];
(* Code to compute Erfc anywhere in the comples plane!*)
erfc[z_]:=If[Re[z]>0,-Exp[-z^2]ψ[-1,z],-Exp[-z^2]ψ[1,z]+2]

Figures that appear in the book and on the webpage

asymerror[z_]:=Abs[ψ[z]-ψexactseries[z]];

TableLogLogPlot[Abs[z asymerror[-z]],{z,100,1000000,100},Axes->False,Frame->True,FrameStyle->Directive[FontSize->15],PlotStyle->Black]

TableLogLogPlot[error[z],{z,100,1000000,500},Axes->False,Frame->True,FrameStyle->Directive[14,Black],PlotStyle->Black,FrameTicks->{{None,{1/100,1/10000,1/1000000}},{{1,10,1000000},None}},FrameLabel->{Style[z,16],Style["Absolute  Error",16,Black]}]

TableLogLogPlot[error[z],{z,100,1000000,100},Axes->False,Frame->True,FrameStyle->Directive[FontSize->15],PlotStyle->Black]
TableLogLogPlot[error[z Exp[I Pi/4]],{z,100,1000000,100},Axes->False,Frame->True,FrameStyle->Directive[FontSize->15],PlotStyle->Black]