# Solving ODE with multiple equilibriums

Consider an ODE of the form: $$u'(t)=-\frac{1}{\varepsilon}u(u-\frac{1}{2})(u-1)$$ with the initial value $$u(0)=u_0.$$ Here $\varepsilon>0$ is a constant.

It is easy to verify that $u\equiv0$ and $u\equiv1$ are two stable equilibriums for the ODE. And the ODE has the property:

• if $u_0<\frac{1}{2}$, then $u(t)\rightarrow0$ as $t\rightarrow\infty$;
• if $u_0=\frac{1}{2}$, then $u(t)\equiv\frac{1}{2}$;
• if $u_0>\frac{1}{2}$, then $u(t)\rightarrow1$ as $t\rightarrow\infty$;

In the problem, $\varepsilon$ is very small. I want to design a numerical scheme that has two property:

• the time step size $\Delta t$ is independent of $\varepsilon$.
• the numerical solutions enjoy has the same asymptotic property as the exact solution.

I have tried the Euler backward scheme: $$u_{n+1} = u_n -\frac{\Delta t}{\varepsilon}u_{n+1}(u_{n+1}-\frac{1}{2})(u_{n+1}-1).$$ However, this scheme, which is a nonlinear algebraic equation for $u_{n+1}$, does not necessarily have a unique solution for any $\varepsilon$ and $\Delta t$.

I want to find a high-order scheme (at least third-order) for this ODE. It seems that the Runge-Kutta method does not fit my purpose.

• Are you sure this is possible? Jul 25, 2016 at 20:36
• @ChrisRackauckas I am not sure...However, it is not easy to prove the nonexistence. I have searched a lot and found no papers on this topic. :( Jul 25, 2016 at 22:41
• @DavidKetcheson I edited the question just now. I want to find a high-order scheme (at least third-order) for this ODE. It seems that the Runge-Kutta method does not fit my purpose. Jul 26, 2016 at 13:16

There are a couple of questions implicit in your post:

How does one deal with non-uniqueness of the algebraic equations generated by any implicit numerical method?

Typically you have a very good initial guess -- the previous time step solution. This will usually ensure that e.g. Newton's method converges to the correct root. For extremely large values of $\Delta t$ this is not true; I don't believe that there is a completely generic way to get around this other than by not taking $\Delta t$ too large. In your particular case you know that the correct solution is real and lies in the interval $[0,1]$, which may help.

How does one ensure high-order convergence for arbitrarily stiff problems?

Use a Runge-Kutta method with high stage order, or use a high order linear multistep method.

• Yes, you are right. There seems no hope to find a high-order method which has nice asymptotic property unconditionally. Jul 26, 2016 at 20:48

If you relax the independence between $\epsilon$ and $\Delta t$ to be "approximate", i.e. use a method which has a really large stable region even if it's not A-stable, then there are plenty of methods which can get good performance. A standard set of methods for this are the Backward Differentiation Formulas (BDF) or the Numerical Differentiation Formulas (NDF). Usually these implementations get quite complicated because the lower order BDF/NDF methods have better stability, so they are normally implemented as variable-timestep and variable-order. For lower tolerance but faster solving one usually uses Rosenbrock methods of orders 2/3 (which is also easier to implement). [Note that although the BDF methods are not A-stable, the stability region always includes the full negative real axis, meaning that for a lone real-valued ODE like the one you have here, these methods are stable for any $\epsilon$ and $\Delta t$. But that fact doesn't generalize to systems of equations or complex-valued ODEs].

So these give implicit methods. However, as you noticed, the implicit methods do not necessarily have a unique solution for what is the next step. How do you choose one? Normally, you assume that your previous timestep is close enough that whatever root finder you use (Newton's method, or others) is going to converge to "the right one" (here you have a positive, a zero, and a negative solution so it's easy to tell when it fails). This has a timestep assumption in there that the derivation of the implicit method doesn't account for, which could lead to some unexpected issues. If this ends up giving you issues, one way to deal with this is to make some kind of a Predictor-Corrector method where you use a Runge-Kutta 4/5 step to produce an estimate for the step (or split the timestep and do a few low-order but stable explicit steps to get an estimate), and use that estimate as the initial guess for Newton's method. This estimate should be closer to the true step, at least enough that Newton's method will now converge to the right place.

Dealing with stiff problems has a lot of research behind it. Not only that, but it has years of software development. So if you're really looking for solvers which make a good compromise between speed and accuracy on hard problems, you should try to stick to tried-and-tested software. You can directly call some of the classics, or you can use wrappers in a higher level language (MATLAB,Python, Julia (wrappers to Hairer and other algorithms coming really soon!)).

• Your advice is great! I think I can try some BDF methods. By the way, the Second Dahlquist barrier is only for the multistep method. For other kinds of method such as Runge-Kutta method, there is no such barrier. Jul 26, 2016 at 20:47
• This answer contains multiple misleading statements. The first sentence of this answer is false. The Dahlquist barrier only refers to linear multistep methods. There are plenty of A-stable methods with order (and stage order) higher than 3. Furthermore, you do not necessarily need A-stability in order to have a stable timestep that is independent of $\epsilon$; you just need a wedge of stability that includes the region where the eigenvalues of your problem are. Jul 27, 2016 at 5:22
• Yes I realized I was being silly. Edited to fix that. Jul 27, 2016 at 5:25
• Your answer is still false. The stability regions of the BDF methods are infinite. Jul 27, 2016 at 5:26
• No, that's false. Only the orders 1 and 2 BDF methods are A-stable. But I should probably mention that it's usually done in a variable order fashion. Jul 27, 2016 at 5:28

I don't quite understand what the problem is with the direct Euler method. Let $f(u) =u(u-1/2)(u-1)$ so the equation is $u'=-(1/ \varepsilon)f(u)$ and use the scheme $u_{n+1}=u_n - (\Delta t / \varepsilon)f(u_n) = F(u_n)$. To study the behavior of the numerical solutions, we study the dynamics of the map $u \mapsto F(u)$. Check that $F(0)=0, F(1/2)=1/2, F(1)=1$. Then $F(u)=u(1-(\Delta t/\varepsilon)(u-1/2)(u-1))=ug(u)$. For $u \in (0,1/2)$, $g(u)<1, g'(u)<0$, so it's strictly decreasing, and thus. $F(u)<u$. Hence for any $u_n \in (0,1/2), u_{n+1}=F(u_n)<u_n,$ so $g(u_{n+1})<g(u_n)<g(u_{n-1})<...g(u_0):=q <1$ as long as $u_0 \in (0,1/2)$. Now $u_{n+1}=F(u_n)=u_ng(u_n)<q u_n < q^2 u_{n-1}<...<q^{n+1} u_0$. Now you see that when $n$ goes to infinity $u_n$ converges to $0$. The rest is handled the same way.

• Yes, if $u_0\in(0,1/2)$, then $u_1<u_0$. However, $u_1$ is not necessarily in the range of $(0,1/2)$. Your proof does not work. Jul 26, 2016 at 14:58
• Isn't $\Delta t$ as small as you want it compared to $\varepsilon$? As long as $\Delta t < 2\varepsilon$ the interval $(0,1/2)$ is invariant under the map $F$. Jul 26, 2016 at 15:59
• @Futurologist The OP is hoping for a method where $\Delta t$ can be taken arbitrarily large. Jul 26, 2016 at 16:28
• Well, independent doesn't necessarily mean unbounded. After all, most likely the actual solutions of the equation converge to step functions as epsilon goes to zero. How is one going to capture that abrupt jump without adjusting the time step? Without that you are going to miss it. Jul 26, 2016 at 17:05