<<RiemannHilbert`;

Direct and inverse monodromy problem

This notebook discusses the solution to the monodromy problem as discussed in Sections 1.4 and 6.4.  Here we consider the 2×2 linear ODE:

Y'(z) =(((Subscript[A, 1])/(z-(Subscript[z, 1])))+((Subscript[A, 2])/(z-(Subscript[z, 2])))+((Subscript[A, 3])/(z-(Subscript[z, 3]))))Y(z),  Y(0) = I

which is defined so that it has a branch cut along [1,∞).    Along this branch cut we have the following relationships between the limit from above and below:

(SubPlus[Y])(x) = (SubMinus[Y])(x)(Subscript[M, 1   ])for 1 ≤ x ≤ 2
(SubPlus[Y])(x) = (SubMinus[Y])(x)(Subscript[M, 2   ])for 2 ≤ x ≤ 3
(SubPlus[Y])(x) = (SubMinus[Y])(x)(Subscript[M, 3   ])for 3 ≤ x

We focus un the case where (Subscript[A, 1])+(Subscript[A, 2 ])+(Subscript[A, 3]) = 0, which implies (Subscript[M, 3])=I, and hence the branch cut only exists along [1,3].

The direct monodromy problem is to relate the input matrices {(Subscript[A, 1]),(Subscript[A, 2])} to the monodromy matrices {(Subscript[M, 1]),(Subscript[M, 2])}.  The inverse monodromy problem asks the opposite: recover {(Subscript[A, 1]),(Subscript[A, 2])} from the monodromy matrices {(Subscript[M, 1]),(Subscript[M, 2])}.  We will see that the direct problem can be solved using simple time integration methods, while the inverse problem can be solved using Riemann–Hilbert problems.

Solve the direct monodromy problem

Consider the following matrices:

(*Setup RHS*)
A1 = {{0.,.1},{-.1,0.}};
A2={{0.,.2},{-.6,0.}};
A3=-(A1+A2);
F[z_,y_]:=(A1/(z-1)+A2/(z-2)+A3/(z-3)).y;

We wish to time integrate the ODE from the origin to each of the two branch cuts.  We will use a simple Runge–Kutta 4 scheme, where the time stepping is conducted in the complex plane to avoid the poles:

(* Use RK4 to integrate the ODE *)
RKStep[F_,t_,dt_][y_]:=Module[{k1,k2,k3,k4},
k1=dt F[t,y];
k2=dt F[t+1/2dt,y+1/2k1];
k3=dt F[t+1/2dt,y+1/2k2];
k4 =dt F[t+dt,y+k3];
y+1/6 (k1+2k2+2k3+k4)
];
(** Time step in the complex plane **)
FuchsianIntegrate[{a_,y0_},b_,dt_]:=Module[{θ,fdt,n,ang,Y,Z,i},
ang=Exp[I Arg[b-a]];
n =Abs[b-a]/dt//Floor;
fdt=Abs[b-a]- n dt;
Y={y0};
Z={a};
For[i=1,i≤n,i++,
Y=Join[Y,{RKStep[F,Z[[-1]],ang dt][Y[[-1]]]}];
Z=Join[Z,{Z[[-1]]+ ang dt}];
];
If[fdt > 0,
Y=Join[Y,{RKStep[F,Z[[-1]],ang fdt][Y[[-1]]]}];
Z=Join[Z,{Z[[-1]]+ ang fdt}];
];
{Z[[-1]],Y[[-1]]}
];
FuchsianIntegrate[{a_,y0_},b_List,dt_]:=Module[{yp,i},
yp={a,y0};
For[i=1,i≤Length[b],i++,
yp=FuchsianIntegrate[yp,b[[i]],dt];
];
yp
];

The following commands calculate (SubPlus[Y]) and (SubMinus[Y]) on the branch cut:

(*Integrate the ODE to compute the monodromy data*)
dt=.001;up=1.;
Y1p=FuchsianIntegrate[{0.,IdentityMatrix[2]},{up I,up I+1.5,1.5},dt]//Last;
Y1m=FuchsianIntegrate[{0.,IdentityMatrix[2]},{-up I,-up I+1.5,1.5},dt]//Last;
Y2p=FuchsianIntegrate[{0.,IdentityMatrix[2]},{up I,up I+2.5,2.5},dt]//Last;
Y2m=FuchsianIntegrate[{0.,IdentityMatrix[2]},{-up I,-up I+2.5,2.5},dt]//Last;

We use this limits to recover (Subscript[M, k])=(Subscript[Y, -]^(-1))(x)(SubPlus[Y])(x) on the corresponding branch cut:

M1=Inverse[M1m].M1p;
M2=Inverse[M2m].M2p;

As a measure of the error in the approximation let’s check that (Subscript[M, 3])=I:

Y3p=FuchsianIntegrate[{0.,IdentityMatrix[2]},{up I,up I+3.5,3.5},dt]//Last;
Y3m=FuchsianIntegrate[{0.,IdentityMatrix[2]},{-up I,-up I+3.5,3.5},dt]//Last;
M3=Inverse[M3m].M3p;
M3-IdentityMatrix[2]//Norm

Set up local parametrix

We wish to set up a Riemann–Hilbert problem for the inverse monodromy problem.  As discussed in Section 6.4 and depicted Figure 6.6, we need three local parametrices (Subscript[P, 1]), (Subscript[P, mid]) and (Subscript[P, 2]) that factor out the jump near these singularities.  These are constructed as follows:

b//Clear;
b[a_,c_,λ_,_][z_]:=Exp[Cauchy[Fun[Log[Abs[λ]]&,Line[{a,c}],2],z]];
LocalParametrix[M_,aa_,bb_,{s1_,s2_}]:=Module[{U,V,Ui},
{V,U}=Eigensystem[M];
U=Transpose[U];
Ui = Inverse[U];
U.DiagonalMatrix[{b[aa,bb,V[[1]],s1][#],b[aa,bb,V[[2]],s2][#]}].Ui&
];

(* Setup local parametrices to deal with non-smooth points where the jump matrix changes abruply*)
eps=$MachineEpsilon;
iM2M1=Inverse[M2].M1;
p1=LocalParametrix[M1,1,2,{-1,-1}];
p2=LocalParametrix[M2,2,3,{1,1}];
pmid=LocalParametrix[iM2M1,1,2,{1,1}];

We confirm numerically that the parametrices have factored out the appropriate jumps:

p1[1.5+I eps]-p1[1.5-I eps].M1//Norm
p2[2.5+I eps]-p2[2.5-I eps].M2//Norm
pmid[1.5+I eps]-pmid[1.5-I eps].iM2M1//Norm

Setup and solve the RH problem

With the singularities removed, we can solve a RH problem that has a smooth solution.  The following code sets up the jump matrix and depicts the contour used, as in Figure 6.6(d):

nnn=40;
rad=.15;
FunDet[a_]:=a[[1,1]]a[[2,2]]-a[[1,2]]a[[2,1]];
RHIndex[a_]:=DomainIntegrate[Table[Log[FunDet[a[[i]]]]',{i,1,Length[a]}]]/(2Pi I);
RHP={Fun[Inverse[p1[#+I eps]]&,Arc[1,rad,{0,Pi}],nnn],
Fun[Inverse[p1[#-I eps]]&,Arc[1,rad,{Pi,2Pi}],nnn],
Fun[M1&,Line[{1+rad,2-rad}],nnn],
Fun[Inverse[pmid[#+I eps]]&,Arc[2,rad,{0,Pi}],nnn],
Fun[M2.Inverse[pmid[#-I eps]]&,Arc[2,rad,{Pi,2Pi}],nnn],
Fun[M2&,Line[{2+rad,3-rad}],nnn],
Fun[Inverse[p2[#+I eps]]&,Arc[3,rad,{0,Pi}],nnn],
Fun[Inverse[p2[#-I eps]]&,Arc[3,rad,{Pi,2Pi}],nnn]
};
DomainPlot[RHP]

We can now use the RHSolve command, verifying convergence by plotting the calculated coefficients on each contour:

sol=RHP//RHSolve;
{DCTPlot/@RHP//Show,DCTPlot/@sol//Show}

We thus can evaluate Y and Y' as a Cauchy transform of sol:

y[z_]:=IdentityMatrix[2]+Cauchy[1,sol,z];
Dy[z_]:=CauchyD[sol,z];
(** normalize at zero to be I **)
AA=Inverse[y[0]];
Y[z_]:=AA.y[z];
DY[z_]:=AA.Dy[z];

To get back the matrices (Subscript[A, 1]), (Subscript[A, 2]) and (Subscript[A, 3]), we evaluate numerically a contour integral associated with Y'(z) ((Y(z))^(-1)):

A[z_]:=DY[z].Inverse[Y[z]];
(*estimate the residues by contour integrals around the singular points*)
rA1=1/(2Pi I)LFun[A[#+I eps]&,Circle[1,2rad],50]//DomainIntegrate;
rA2=1/(2Pi I)LFun[A[#+I eps]&,Circle[2,2rad],50]//DomainIntegrate;
rA3=1/(2Pi I)LFun[A[#+I eps]&,Circle[3,2rad],50]//DomainIntegrate;

We now verify we have correctly recovered the input matrices:

(*Check the errors*)
A1-rA1//Norm
A2-rA2//Norm
A3-rA3//Norm

Plot solution

Here we use Y[z] to plot the solution in the complex plane:

erad=(1-10000eps);X//Clear;
X[z_]:=X[z]=({{\[Piecewise], ({{Y[(z-1)erad+1].p1[z], Abs[z-1]<rad}, {Y[(z-2)erad+2].pmid[z], Abs[z-2]<rad && Im[z]≥0}, {Y[(z-2)erad+2].pmid[z].Inverse[M2], Abs[z-2]<rad && Im[z]<0}, {Y[(z-3)erad+3].p2[z], Abs[z-3]<rad}, {Y[z], True}})}})

ReImTablePlot[X[1.9+I yy][[1,1]],{yy,-1,1,.01}]
ReImTablePlot[pmid[1.9+I yy][[1,1]],{yy,-1,1,.01}]

mStyle[x_]:=Style[x,14,Bold];
P1=TablePlot3D[X[x + I yy][[1,1]]//Re,{x,0-Pi/100,4,.1},{yy,-.5+10eps,.5,.1/8},AxesLabel->mStyle/@{Re[z],Im[z],Re[(Subscript[Y, 11])[z]]}]
P2=TablePlot3D[X[x + I yy][[1,1]]//Im,{x,0-Pi/100,4,.1},{yy,-.5+10eps,.5,.1/8},AxesLabel->mStyle/@{Re[z],Im[z],Im[(Subscript[Y, 11])[z]]}]