Solving ODE with multiple equilibriums

Question

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.

Answer 1

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.

Answer 2

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

Answer 3

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.