# Analytic continuation of an ODE solution with ParametricNDSolveValue

I am working on the analytic continuation of a solution of the $$2\times 2$$ system $$dT=TP,\quad T(1)=Id$$ where $$P$$ is of the form $$P=\begin{pmatrix} 0 & t^{-1}\\ \frac{t}{4}+(t-1)^2\left(\frac{r z^2}{4}+\frac{r}{4 z^2}+\frac{s z}{2}+\frac{s}{2 z}+u\right) & 0\end{pmatrix}\frac{dz}{z},$$ where the variable $$z$$ is complex, $$t$$ is a complex parameter of modulus $$1$$ and $$r,s,u$$ are real numbers of our choice (more or less).

The idea of analytic continuation is that one can extend analytically the solution along any closed path enclosing the irregular singularity at $$z=0$$, but after that the solution in general won't be the same. In particular, in this setup, when you go around $$z=0$$ once then $$z\mapsto z+2\pi i$$, so this would be the kind of point at which I'd like to evaluate the solution.

I don't know if it's possible to find this 'analytically continued' solution numerically in Mathematica. I've been trying with ParametricNDSolveValue as follows:

r = 1/25; s = 0; u = 1/16;
P[z_] = (1/z){{0, 1/t},
{t/4 + (t - 1)^2 (u + r/(4 z^2) + s/(2 z) + (s z)/2 + (r z^2)/4), 0}}];
pf = ParametricNDSolveValue[{T'[z] == T[z].P[z],
T[1] == IdentityMatrix[2]}, T[1 + 2 Pi I], {z, 0.5, 2 \[Pi]}, t];


Mathematica deals with solving this so then I ask for the value at $$\;t=1$$ getting

pf[1]
(* InterpolatingFunction[{{0.5, 8.28319}}, <>][1 + 2 I \[Pi]] *)


with which I don't know how to deal with.

If I change the value $$1+2\pi i$$ for a real number, then Mathematica gives a solution that I can treat. So my guess is that for ParametricNDSolveValue evaluating at complex points is an issue.

Then I thought that since we can pick any closed curve around $$0$$, I could change my equation to polar (with radius $$1$$) and evaluate at $$2\pi$$ in ParametricNDSolveValue.

Then I run

Q[\[Theta]_] = FullSimplify[I Exp[I \[Theta]] P[Exp[I \[Theta]]];
pfpolar = ParametricNDSolveValue[{S'[\[Theta]] == S[\[Theta]].Q,
S[0] == IdentityMatrix[2]}, S[2 \[Pi]], {\[Theta], 0, 2 \[Pi]}, t];


which works nicely...

pfpolar[1]
(* {{-1. + 0. I, 0. + 7.06964*10^-8 I}, {0. + 1.76741*10^-8 I, -1. + 0. I}} *)


...but there is something wrong: I have proved that this matrix has a maximum trace of $$2$$ at $$t=1$$ and then lower values never less than $$-2$$, while instead in Mathematica I get

Tr[pfpolar[1]]
Tr[pfpolar[I]]
Tr[pfpolar[-1]]
Tr[pfpolar[-I]]
(* -2. + 0. I *)
(* -1.21463 + 0. I *)
(* 1.96848 + 0. I *)
(* -1.21463 + 0. I *)


It is extremely weird as it looks like I'm getting some short of singed changed solution... I don't know if I'm using wrongly ParametricNDSolveValue or making any mistake you can spot. Any help is very welcome.

EDIT This question has been tagged as a possible duplicate. I will try to explain in what follows why I still can't solve my problem.

I have tried to use the interpolation methods in the related question, first by computing the solution without evaluating, that is,

pf = ParametricNDSolveValue[{T'[z] == T[z].P[z],
T[1] == IdentityMatrix[2]}, T, {z, 0.5, 1.5}, t];


Then this pf depends on $$z$$ and on the parameter $$t$$. In order to use the intepolation method shown in the related question, I need to define a function depending on the one variable for which the interpolation will be done. The problem is that I don't know how to do this without 'losing' the dependence on the parameter $$t$$. If I do something like

func[z_] := pfz[t][z];
func[1]
ParametricFunction[ <> ][t][1]


I get a parametric function depending on $$t$$, but this is not what I want since I would like to construct a function depending on $$t$$, once I can evaluate at $$z=1+2\pi i$$. I have tried also

func[z_, t_] := pfz[t][z];


But then surprisingly (at least for me)

func[1, 0.9]
pfz[1][0.9]
{{1., 0.}, {0., 1.}}
{{1.00139, -0.105409}, {-0.0263523, 1.00139}}


func and pfz don't give the same output when evaluated at the same points.

I still wanted to try the StatisticsLibraryBarycentricInterpolation suggested in the related question so I've forced my parameter $$t$$ to depend on $$z$$ (this dependence is not true, so even if this worked I wouldn't use it for my problem). I define then

func[z_] := pfz[Exp[I z]][z][[1]][[1]]; (* take for the example only
the first entry in the matrix solution *)


and run the interpolation method

(*Fourier interpolation on a complex disk*)
nn = 64; (*number of interpolation points*)
z0 = 1; (*center of circle*)
rr = 1/2; (*radius of circle*)
wp = MachinePrecision;(*working precision*)
tj = 2 Pi*Range[0, nn - 1]/nn;
zj = N[z0 + rr Exp[I tj], wp]; (*interpolation nodes*)
fj := func /@ zj; (*function values on nodes*)if :=
StatisticsLibraryBarycentricInterpolation[zj, fj,
Weights -> Exp[2 Pi I Range[0., nn - 1]/nn]];


If look inside if the output is

if

Part::partd: Part specification (1.49759 +0.0490086 I)[[1]]
is longer than depth of object. >>

Part::partd: Part specification (1.49039 +0.0975452 I)[[1]]
is longer than depth of object. >>

Part::partd: Part specification (1.47847 +0.145142 I)[[1]]
is longer than depth of object. >>

General::stop: Further output of Part::partd will be suppressed
during this calculation. >>


I guess this is because the interpolation process is evaluating the result of my initial ParametricNDSolveValue at complex values, which was my problem originally. ParametricNDSolveValue seems to behave badly with this evaluations, at least for this function. Actually, if I evaluate myself at a complex value

func[I]
Part::partd: Part specification I[[1]]
is longer than depth of object. >>


same happens.

Just for the sake of completeness, looking inside if gave also a long output of the form

StatisticsLibraryBarycentricInterpolation[{1.5, 1.49759 + 0.0490086 I,
1.49039 + 0.0975452 I, 1.47847 + 0.145142 I,
(* keeps going *)
{1.00743 + 7.53749*10^-19 I, (1.49759 + 0.0490086 I)[[1]], (1.49039 + 0.0975452 I)[[1]],
(* keeps going *)
Weights -> {1. + 0. I, 0.995185 + 0.0980171 I, 0.980785 + 0.19509 I,
(* keeps going *)
0.92388 - 0.382683 I, 0.95694 - 0.290285 I, 0.980785 - 0.19509 I, 0.995185 - 0.0980171 I}]


which I can't really interpret.

Any help with this or with the polar version above is much appreciated.

• Possible duplicate of Complex continuation of an interpolated function – Carl Woll May 7 at 23:16
• @CarlWoll Could you help me to understand how I can apply the interpolation methods in those answers? What I'm expecting as output of `ParametricNDSolveValue' is a matrix solution of the ODE already evaluated at a point and that depends on the parameter $t$... – Edu May 8 at 20:58