# Numerical approximation of non-unique ODE solution

Suppose you have an autonomous first order IVP $$x'(t)=f(x(t)), \quad x(t_0)=x_0$$ where $x(t)\in\mathbb{R}$ and $x_0$ is a fixed point of the equation. Clearly, $$x^0(t) \equiv x_0$$ is a solution to this ODE. However, suppose there exists some other non-trivial solution $x^1(t)$ and that $f'(x(t))$ does not exist when $t=t_0$.

Is there a numerical procedure which will approximately solve the IVP and give $x^1(t)$ as a result?

Euler's method, $x(t+h) = x(t)+hf(x(t))%$, clearly converges to the trivial solution $x^0$, which can be seen via a simple induction argument. Other Taylor series methods are out since $f$ is not differentiable at $t_0$. It looks like basically any numerical procedure which converges to a solution of the the equation will converge to $x^0$.

Example: $x'(t) = \sqrt{|x(t)|}, \quad x(0)=0$

$x^0(t)\equiv 0$ is one solution, and $x^1(t) = \frac{1}{4}t^2\text{sign}(t)$ is another (Grimshaw, "Nonlinear Ordinary Differential Equations", pg. 21 ex. 3). Any numerical procedure I throw at it converges to $x^0$, but I want one which approximates $x^1$.

Problems like this are to my knowledge often not discussed in numerical analysis. The problem you look at is clearly ill-posed.

However, there is the fairly new branch of mathematics called "probabilistic numerics" that tries to combine methods from statistics and numerical analysis. They in principal propose to reformulate the problem to a well-posed statistical one. This has been done for all kinds of numerical analysis, including the discretisation of ODEs. I suggest you look at Conrad et al. 2016 for probabilistic ODE solvers and maybe at http://www.probabilistic-numerics.org for the bigger image.

• Then why was it on my final exam :( May 3, 2018 at 22:41
• Was the general question on your final exam or the question with the square-root? May 4, 2018 at 6:59
• There was a different equation with the same issue. I believe it was y'=y^(1/3). May 6, 2018 at 23:42

It looks like are trying explicit numerical methods. By their nature they can only give a single value at the new time -- multiple solutions are not obtained when sing explicit methods $$y_{n+1} = g(y_n)$$. Implicit methods involve solving an implicit equation for $$y_{n+1}$$ at each time-step. That should allow for multiple solutions since rootfinding methods must be able to handle multiple solutions.

Try Backward Euler (at first step): $$y_1 - y_0 = dt\; y_1^{1/3}$$. Factoring on the left side yields $$(1-dt\;y_1^{-2/3})y_1=y_0$$. We take $$y_0=0$$ and get two solutions for $$y_1$$: $$y_1=0$$ and $$y_1=dt^{3/2}$$. The second solution approximates the solution you want. Depending on the rootfinder you use, you might get either of these solutions.

In this case a semi-implicit solution that weights the right side evaluation as 1/3 at old time and 2/3 at new time will give zero numerical error on the first time-step for any step size - but that is specific to this example.