Numerical scheme to solve ODE with singularity

Question

I've been struggling with the following. I'd like to find a numerical scheme to solve the following ODE:

$$ x'' + \frac{1}{x^2} = 0 $$

Most schemes work quite well when $x$ isn't too close to zero, but I need one that is robust to traversing the origin. E.g., I expect that if my initial conditions are $(x, x') = (1, 0)$, the solution would be a kind of oscillator between $1$ and $-1$.

I'm not sure this is even possible to find? I've done some research and tried to find a solution myself, but I'm hitting a wall. Does anyone have some insights on this problem?

Thanks for your help!

Answer 1

Before you look for a numerical scheme, ask whether there is a solution at all. The differential equation is undefined at $x=0$, so there is really no such thing as "traversing the origin" for its solutions. In fact, since $x'' < 0$ everywhere, it is impossible to continue a solution after it hits the origin with velocity $-\infty$, unless you make the velocity jump to $+\infty$. You could, for example, bounce off the origin (so $x(t_0+t) = x(t_0-t)$ where $t_0$ is a time you hit the origin). But if you somehow made it to some $x < 0$ with $x' < 0$, $x'$ will stay negative so you can never go back.

EDIT: In fact, by "conservation of energy", $E(x,v) = \frac{v^2}{2} - \frac{1}{x}$ is conserved, where $v = dx/dt$ is the velocity.
Trajectories in the $(x,v)$ plane look like this.

For positive $x$, solutions can come from $x \to 0, v \to +\infty$ at some finite time in the past and go back with $x \to 0, v \to -\infty$ at some finite time in the future, or come from $(0, +\infty)$ and go to $x \to +\infty$ as $t \to +\infty$, or come from $x \to +\infty$ as $t \to -\infty$ and go to $(0,-\infty)$. For negative $x$, solutions must come from $x \to -\infty, v \to const > 0$ as $t \to -\infty$ and go back to $x \to -\infty$, $v \to -const < 0$ as $t \to +\infty$ without ever approaching $x=0$.

Answer 2

Before comment for numerical solution it is of interest to see to what leads the analytic solution : $$\frac{d^2x}{dt^2}+\frac{1}{x^2}=0$$ $$2\frac{d^2x}{dt^2}\frac{dx}{dt}+\frac{2}{x^2}\frac{dx}{dt}=0$$ $$\left(\frac{dx}{dt}\right)^2-\frac{2}{x}=c$$ With initial condition $(x,x')=(1,0)$ at $t=0$ , equivalently $\begin{cases} x(0)=1 \\ x'(0)=0\end{cases}$ $$0^2-\frac21=c\quad\implies\quad c=-2$$ $$\left(\frac{dx}{dt}\right)^2=\frac{2}{x}-2\quad\implies\quad 0\leq x< 1$$ $$\frac{dx}{dt}=\pm\sqrt{\frac{2}{x}-2}$$ $$dt=\pm\sqrt{\frac{x}{2(1-x)}}dx\qquad 0\leq x< 1$$

Let $x(t)=\cos^2(\alpha(t))$ $$dt=\pm \sqrt{2}\cos^2(\alpha)d\alpha$$ $$t=\pm \sqrt{2}\int\cos^2(\alpha)d\alpha$$ $$t=\pm\frac{1}{\sqrt{2}}\big(\alpha+\sin(\alpha)\cos(\alpha) \big)+\text{constant}$$ The constant is eliminated with condition $x(0)=1$ then $\alpha(0)=0$.

The analytic solution expressed on parametric form is : $$\boxed{\begin{cases} t=\pm\frac{1}{\sqrt{2}}\big(\alpha+\sin(\alpha)\cos(\alpha) \big)\\ x=\cos^2(\alpha) \end{cases}}$$ One cannot express $x(t)$ on closed form because the equation $\quad t(x)=\pm\frac{1}{\sqrt{2}}\left(\cos^{-1}(\sqrt{x})+\sqrt{x(1-x)} \right)\quad$ cannot be inverted on the form of a finite number of elementary functions.

The above parametric form is the simplest to study the solution, which clearly is periodic. (Period$=\pi\sqrt{2}$ ).

COMMENT about numerical solution.

I suppose that the starting point is for $t=0$ with $x=1$ and $x'=0$ then $x''=-\frac{1}{1^2}=-1$.

There is no difficulty with successive small increments of $t$ to compute the successive values of $x''$ then $x'$ and next $x$.

The difficulty arrises when $t\simeq\frac{\pi}{\sqrt{2}}$ corresponding to $x\simeq 0$ because $x''$ becomes big. A possible way to overcome the difficulty is to change the origin of $t$ and restart the process backwards from another point, for example $t=\pi\sqrt{2}$, then forwards up to $t\simeq 3\frac{\pi}{\sqrt{2}}$, and so on.

But definitively the simplest method is to use the parametric equation and draw $\big(t(\alpha)\:,\:x(\alpha)\big)$ for $\alpha$ from $0$ to any large value of $\alpha$.