# Using NDSolve on the Painlevé equations

In an earlier question of mine, I was looking for a way to handle certain kinds of singularities (poles) when using NDSolve. Michael E2's answer, which relied on projective geometry was the most upvoted answer.

I have tried generalizing his method to systems of second order \begin{align*} u'&=f(u,v,t), \\ v'&=g(u,v,t), \\ u(t_0)&=u_0, \\ v(t_0)&=v_0,\end{align*} by considering the pair of dependent variables $$(u,v)$$ as an element of $$\mathbb{P}^1 \times \mathbb{P}^1$$, that is \begin{align*}u=\frac{u_1}{u_2}, \\ v=\frac{v_1}{v_2}.\end{align*} This transforms the system to \begin{align*}u_1' u_2-u_1 u_2'&=u_2^2 f \left( \frac{u_1}{u_2}, \frac{v_1}{v_2},t \right), \\ v_1' v_2-v_1 v_2'&=v_2^2 g \left( \frac{u_1}{u_2},\frac{v_1}{v_2},t \right). \end{align*} In order to handle the under-determined nature of these equations, we add the constraints \begin{align*}\frac{\mathrm{d}}{\mathrm{d} t} \left( u_1^2+u_2^2 \right)&\equiv 2u_1 u_1'+2u_2 u_2'=0, \\ \frac{\mathrm{d}}{\mathrm{d} t} \left( v_1^2+v_2^2 \right)&\equiv 2v_1 v_1'+2v_2 v_2'=0. \end{align*}

As a test, here is the code to solve the $$\sec$$-$$\tan$$ system \begin{align*}u'&=uv \\ v'&=u^2, \\ u(0)&=1, \\ v(0)&=0. \end{align*}

sol = NDSolve[{u1'[t] u2[t] - u1[t] u2'[t] ==
u2[t]^2 *u1[t]/u2[t]*v1[t]/v2[t],
v1'[t] v2[t] - v1[t] v2'[t] == v2[t]^2* u1[t]^2/u2[t]^2,
u1[t] u1'[t] + u2[t] u2'[t] == 0, v1[t] v1'[t] + v2[t] v2'[t] == 0,
u1[0] == 1, u2[0] == 1, v1[0] == 0, v2[0] == 1}, {u1[t], u2[t],
v1[t], v2[t]}, {t, -4, 4}, WorkingPrecision -> 30,
MaxSteps -> Infinity]


This works like a charm, in the sense that the quotients $$u_1/u_2$$ and $$v_1/v_2$$ behave like $$\sec$$ and $$\tan$$. However, when I tried handling more complex problems, such as the first Painlevé equation \begin{align*}u'&=v, \\ v'&=6u^2+t,\end{align*} NDSolve couldn't integrate beyond the poles (with the error: NDSolve::ndsz: At t == -1.25774722466261350439119698118, step size is effectively zero; singularity or stiff system suspected.). Here is my (failed) attempt at doing this:

sol = NDSolve[{u1'[t] u2[t] - u1[t] u2'[t] == u2[t]^2*v1[t]/v2[t],
v1'[t] v2[t] - v1[t] v2'[t] == v2[t]^2*(6 u1[t]^2/u2[t]^2 + t),
u1[t] u1'[t] + u2[t] u2'[t] == 0, v1[t] v1'[t] + v2[t] v2'[t] == 0,
u1[0] == 1, u2[0] == 1, v1[0] == 0, v2[0] == 1}, {u1[t], u2[t],
v1[t], v2[t]}, {t, -4, 4}, WorkingPrecision -> 30,
MaxSteps -> Infinity]


My question:

Is there any way to adapt Michael E2's method in order to produce a numerical integrator for the first Painlevé equation that can integrate beyond poles?

Thank you!

• The solution is to integrate around the poles. I wish there were a suboption to NDSolve's RungeKutta such as IntegrationContour to manually specify how to reach the end point... Commented May 30, 2018 at 1:22
• I think trying a different constraint may help. Commented Jun 7, 2018 at 3:10
• One approach I have used when I was investigating the Painlevé equations was to construct the (nonlinear!) ODE for the associated tau functions, which are analytic. As I recall, the problem was that they also had considerable exponential behavior that made the approach usable only for moderately-sized arguments. Commented Jul 25, 2019 at 13:52

By the method in this article A method for calculating the Painlevé transcendents, poles of second-order ODE can be passed by auxiliary functions u[x] and v[x]. Different from Michael E2's answer, for the $$k$$th order pole of y[x], the expressions of u[x] and v[x] are $$y(x)/y^{\prime}(x)$$ and $$(y^{\prime}(x))^k/(y(x))^{k+1}$$. Since all of the poles of the first type of Painlevé equations are of second order, the auxiliary functions are chosen to be $$u(x)=\frac{2y(x)}{y^{\prime}(x)},\quad\nu(x)=\frac{y^{\prime2}(x)}{4y^3(x)},$$ and correspondingly, $$y(x)= \frac{1}{u(x)^2 v(x)},\quad y'(x)=\frac{2}{u(x)^3 v(x)}.$$ The relations between $$u(x)$$, $$v(x)$$, and $$y(x)$$ suggest that auxiliary functions u[x] and v[x] cannot pass by the roots of y[x] and y'[x], and that means we need to solve the equations in two methods region by region.

In auxiliary functions u[x] and v[x], the original ODE is changed to

{6/v[x] + x  u[x]^4  v[x] + 2  Derivative[1][u][x] == 4,
6 + x  u[x]^4  v[x]^2 == 6  v[x] + u[x]  Derivative[1][v][x]}


First, we define two major functions that solve the ODE in different ways:

soly[x0_, y0_, yp0_, xmin_ : -10, xmax_ : 10] :=
Reap@First@
NDSolve[{y''[x] == 6 y[x]^2 + x, y[x0] == y0, y'[x0] == yp0,
WhenEvent[y[x]^2 > 10^9 || x == xmax || x == xmin,
Sow[{x, y[x], y'[x]}]; "StopIntegration"]}, y, {x, xmin, xmax}]

soluv[x0_, u0_, v0_, xmin_ : -10, xmax_ : 10] :=
Reap@First@
NDSolve[{6/v[x] + x  u[x]^4  v[x] + 2  Derivative[1][u][x] == 4,
6 + x  u[x]^4  v[x]^2 == 6  v[x] + u[x]  Derivative[1][v][x],
u[x0] == u0, v[x0] == v0,
WhenEvent[u[x]^2 + v[x]^2 > 10^9 || x == xmax || x == xmin,
Sow[{x, u[x], v[x]}]; "StopIntegration"]}, {u, v}, {x, xmin,
xmax}]


and the transformation of the boundary conditions:

ytouv[x0_, y0_, yp0_] := {x0, 2 y0/yp0, yp0^2/(4 y0^3)};
ytouv[x_List] := ytouv @@ x;
uvtoy[x0_, u0_, v0_] := {x0, 1/(u0^2 v0), 2/(u0^3 v0)};
uvtoy[x_List] := uvtoy @@ x;


Now we can solve the equations as

sy[1] = soly[0, 1, 0]
suv[1] = soluv[Sequence @@ ytouv@sy[1][[2, 1, 2]]]
sy[2] = soly[Sequence @@ uvtoy@suv[1][[2, 1, 2]]]
suv[2] = soluv[Sequence @@ ytouv@sy[2][[2, 1, 2]]]

Show[Plot[
y[x] /. sy[1][[1]], {x, sy[1][[2, 1, 1, 1]], sy[1][[2, 1, 2, 1]]},
PlotRange -> {0, 100}],
Plot[y[x] /. sy[2][[1]], {x, sy[2][[2, 1, 1, 1]],
sy[2][[2, 1, 2, 1]]}, PlotRange -> {0, 100}],
PlotRange -> Automatic, Frame -> True]



To make the above process automatic, we can use the following code

Clear[sy, suv]

sy[1] = soly[0, 1, 0, -10, 10];
Last@Reap@
Do[If[sy[i][[2, 1, 2, 1]] > 10, Sow[i]; Break[]];
suv[i] = soluv[Sequence @@ ytouv@sy[i][[2, 1, 2]], -10, 20];
sy[i + 1] =
soly[Sequence @@ uvtoy@suv[i][[2, 1, 2]], -10, 20];, {i, 1, 100}]

Last@Reap@
Do[If[sy[-i + 1][[2, 1, 1, 1]] < -10, Sow[-i + 1]; Break[]];
suv[-i] = soluv[Sequence @@ ytouv@sy[-i + 1][[2, 1, 1]], -20, 20];
sy[-i] = soly[Sequence @@ uvtoy@suv[-i][[2, 1, 1]], -20, 20];, {i,
0, 100}]

(* {{19}}*)

(* {{-9}}*)



and plot

Show[Table[
Plot[y[x] /. sy[i][[1]], {x, sy[i][[2, 1, 1, 1]],
sy[i][[2, 1, 2, 1]]}, PlotRange -> {-40, 40}], {i, -9, 19}],
PlotRange -> {Automatic, Automatic}, Frame -> True]


• Interesting. One small question: why do you add the coefficients 2 and 4 in the transformation? Is there a deep reason? BTW, you can use DSolveChangeVariables for the transformation :) . Commented Apr 24 at 9:22
• The coefficients 2 and 4 is from the original paper to simplify transformations and the ultimate formulas. @xzczd Commented Apr 24 at 10:54