Numerical scheme to solve ODE with singularity

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?

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$$.

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}$$ ).

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$$.

• Does this particular analytic/parametric approach generalize to, say, a fourth-order equation? Does it have a name so I can look it up? It looks promising for the problem I mentioned here. Commented May 23 at 17:06
• @Daniel Hatton . The above approach is simply a convenient change of variable followed by solving for the inverse function. So nothing new. The key is the change of vatiable which simplifies the ODE to an ODE of separable kind. The case of the fourth-order ODE is much more complicated. There is no evident convenient change of variable. Probably numerical solving is the only way. Commented May 25 at 8:44
• I was really mainly talking about the order-reduction trick in the first three steps. But on reflection, I think that relies on the equation being autonomous, and mine isn't. (It was, but I already reduced it from fifth-order to fourth-order at the price of sacrificing autonomous-ness (and homogeneous-ness).) Commented May 25 at 11:42
• @Daniel Hatton. I am pessimistic about analytic solving because even in the simplified case $P=0$ apparently one cannot find the solution in terms of a finite number of standard functions. A fortiori with $P\neq 0$. Commented May 25 at 17:12