Using NDSolve on the Painlevé equations

Question

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!

Answer 1

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]