2
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 9, 2021 at 15:34
  • $\begingroup$ This is something covered in any introductory numerical analysis textbook. $\endgroup$ Dec 10, 2021 at 11:24

1 Answer 1

3
$\begingroup$

First, you may define a differentiable interpolant for the data produced by the ODE stepper, and then plot the residual

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.