The equation $\dot{x}=x^3$.¶

If we try to solve $\dot{x}=x^3$ we find that the solution goes to infinity at a finite time. In this notebook we investigate the details.

In [2]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

We first solve $\dot{x}=x^3$ for $0\leq t\leq 0.1$ with $x=1$ at $t=0$. The result is well-behaved.

In [3]:
t = np.linspace(0, 0.1, 1000)
sol = solve_ivp(lambda t, x: x**3, [0, t[-1]], [1], t_eval=t)
plt.plot(sol.t, sol.y[0])
Out[3]:
[<matplotlib.lines.Line2D at 0x29c99ae5ed0>]
No description has been provided for this image

We now try to go to $t=1$. We see that the solution escapes to $\infty$ at $t=0.5$.

In [4]:
t = np.linspace(0, 1, 1000)
sol = solve_ivp(lambda t, x: x**3, [0, t[-1]], [1], t_eval=t)
plt.plot(sol.t, sol.y[0])
Out[4]:
[<matplotlib.lines.Line2D at 0x29cab303d50>]
No description has been provided for this image

We define a function cutoff(a) which estimates the time at which the solution escapes to $\infty$ if we start with $x=a$ at $t=0$.

In [5]:
def cutoff(a):
    t = np.linspace(0, 1, 10000)
    sol = solve_ivp(lambda t, x: x**3, [0, 1], [a], t_eval=t)
    return sol.t[-1]

for a in np.arange(1,10):
    print(f"Cutoff time for a={a} is t={cutoff(a)}")

Cutoff time for a=1 is t=0.4999499949995
Cutoff time for a=2 is t=0.12491249124912492
Cutoff time for a=3 is t=0.05550555055505551
Cutoff time for a=4 is t=0.031203120312031204
Cutoff time for a=5 is t=0.019901990199019903
Cutoff time for a=6 is t=0.013801380138013802
Cutoff time for a=7 is t=0.0102010201020102
Cutoff time for a=8 is t=0.007800780078007801
Cutoff time for a=9 is t=0.006100610061006101
In [9]:
aa = np.linspace(1, 10, 100)
tt = np.array([cutoff(a) for a in aa])
plt.plot(aa, tt, label='cutoff time')
# plt.plot(aa, 1/(2 * aa ** 2),'r')
plt.xlabel('a')
plt.ylabel('t')
plt.legend()
Out[9]:
<matplotlib.legend.Legend at 0x29cab3b5090>
No description has been provided for this image

After a bit of experimentation we see that the graph of 1/np.sqrt(2*cutoff(a)) is a straight line of slope $1$, so 1/np.sqrt(2*cutoff(a))=a, which rearranges to give cutoff(a)=1/(2a**2).

In [11]:
plt.plot(aa,1/np.sqrt(2*tt))
Out[11]:
[<matplotlib.lines.Line2D at 0x29cab62f310>]
No description has been provided for this image

We now define a function to calculate the cutoff formula and check that it agrees with our earlier approach within an error of about $10^{-4}$.

In [14]:
def cutoff_formula(a):
    return 1/(2*a**2)

np.max(np.abs([cutoff(a) - cutoff_formula(a) for a in np.arange(1,10)]))
Out[14]:
9.800980098009782e-05

We can obtain the above formula analytically as follows. We can separate variables in the equation $\dot{x}=x^3$ to get $dt=x^{-3}dx$, and then integrate to get $t=c-x^{-2}/2$ for some constant $c$. We want $x=a$ when $t=0$ which gives $c=a^{-2}/2$. We now have $t=(a^{-2}-x^{-2})/2$ which gives $x^{-2}=a^{-2}-2t$ and then $$ x = (a^{-2}-2t)^{-1/2} = \frac{a}{\sqrt{1-2a^2t}}. $$ This goes to infinity when $1-2a^2t=0$ i.e. $t=1/(2a^2)$.

In [ ]: