# Algorithm to numerically determine whether my computed solution for a 1st order ODE is stable/unstable?

We were given an assignment where we had to determine the numerical solution of Dahlquist's equation $$\dot x$$ = $$\lambda x$$, ($$\lambda$$ = $$-7$$) for time steps $${0.5,0.25,0.125}$$ using explicit euler method. When I coded it and plotted the solution on MATLAB, I saw that the solution unphysically oscillates for time steps $$0.5$$ and $$0.25$$. I was thinking whether I could devise an algorithm which predicts whether a solution is stable for any general ODE $$\dot x$$ = $$f(x)$$.

Any ideas of a robust algorithm that could work for general ODEs of the type above?

P.S : I was thinking that I could maybe see whether the slope of numerical solution and the analytical derivative at the solution points are of the same sign (i.e if $$\frac{x^{n+1} - x^n}{\delta t}f(x^n)<0$$ at any point ($$t^n$$,$$x^n$$) implies there are unphysical oscillations) but it would not work when both numerical and analytical solutions monotonically increase/decrease (with the numerical method going to very large positive or negative values for increasing time $$t$$).

• You need that $λh$ is in the stability region for any eigenvalue of the Jacobian. Or simplified, that $Lh<2$ where $L$ is a (local) Lipschitz constant of $f$. Usually, you also get stability if the error estimator of your choice gives sensible values. Dec 9, 2021 at 15:34
• This is something covered in any introductory numerical analysis textbook. Dec 10, 2021 at 11:24

$$\Delta(t) := \tilde{y}'(t) - f(\tilde{y},t)$$
where $$\tilde{y}$$ is the interpolation of the solution skeleton. (In matlab, isn't the default interpolator for solution skeletons pchip? This should do; though cubic Hermite interpolation might be better as $$\tilde{y}'(t)$$ is so easily computable.) This will almost certainly show you your stepper didn't do what you wanted. If you are looking to programmatically determine whether to resolve with a different method, I might compute $$\left\| \Delta \right\|_{\infty}$$.
Using the residual is in IMO the easiest way for you to obtain the information you want. However, strictly speaking I have not answered your question as stated, since you want to know if there is an a posteriori method to determine if a stepper's interaction with a rhs $$f$$ leads to it being unstable. It turns out that if a stepper is unstable on the initial value problem, it is very often the case that it is stable on the terminal value problem. You solve the equation backwards in time from the final location of your IVP solution. If these two solutions are not roughly the same, then you know the stepper is unstable in one or the other direction. I would use the residual to determine which is which.