Using backward and forward Euler method to solve a certain stiff ODE

Question

When using the backward and forward Euler methods to solve a certain stiff differential equation, what criteria does one look at before drawing the conclusion that one is more stable than the other?

Answer 1

General observations

Convergence

If $L$ is a Lipschitz constant for the system that is valid for the medium term, then with $Lh\le 1$ both methods give with good probability (but not certainty) sensible results, with $Lh\le 0.1$ or smaller you should also get quantitatively correct results.

Stability

The difference occurs when $Lh>1$. The implicit method will continue to give qualitatively correct results, while the explicit Euler method will give wild oscillations with increasing amplitude.

Solution uniqueness for the implicit step equation

Under good conditions the equation $z-hf(z)=x$, with $x=x_n$ of the implicit Euler method, $z=x_{n+1}$, should only have one solution $x+O(h)$ that lies on a continuous curve $h\mapsto z$. All the other solutions of the non-linear system have distances $\sim 1/h$, that is, are far away.

When $h$ becomes too large, so that there is no meaningful difference in distances $h$ and $1/h$, it can happen that the implicit equation has multiple solutions close-by and the implicit solver converges to a wrong one, that connects to the infinite sphere, not to $x$.

Exploration of the last points on the "flame dynamic" example

Statement of the IVP

$$ x'(t)=f(x(t))=x(t)^2-x(t)^3,~~ 0<x(0)=\delta\ll 1, t\in[0,2/δ] $$ with a stable equilibrium at $x=1$ that the solution jumps to at about $t=1/δ$

Explicit Euler

Reaching a point close and below the equilibrium, $x=1-\varepsilon$, the Euler step is in first order $$ x_{+1}=x+hf(x)=1-ε+hε+O(hε^2). $$ So for $h>1$ the $x_{+1}$ will end up on the other side of the equilibrium and for $h>1$ its distance from the equilibrium will be larger against the stable nature of the equilibrium.

Implicit Euler

The step equation to solve is $$ z-hz^2(1-z)=x\iff hz^3-hz^2+z-x=0 \\ \iff (hz^2+1)(z-1)+1-x=0 $$ For $h\approx 0$ the roots will approximately be $$ z=1+\frac{x-1}{1+hz^2}=1+\frac{x-1}{1+hx^2}+O(h^2) $$ for the close-by root and $$ \sqrt{h}z-i=\frac{x-1}{(z-1)(\sqrt{h}z+i)} =\sqrt{h}\frac{x-1}{2i(i-\sqrt{h})}+O(h) \\ \implies z=\frac{i}{\sqrt{h}}-\frac{x-1}2+O(\sqrt{h}) $$ and its complex conjugate for the far-away roots. For a purely real solver there is no danger of converging to the complex roots. Nevertheless it demonstrates the idea of the spurious roots coming into interference range

Step size control for forward Euler

The motivating blog post made a point about variable-step explicit methods. Close to an equilibrium these tend to extend the step size towards the boundary of the stability region, here $2$, and a little beyond, then bouncing back, etc. A complication is that at the stability boundary the error estimate becomes unreliable, luckily in most cases in direction of overestimating the error.

For the Euler method with the Heun method for the second order error estimator this can be implemented with a relatively short script:

problem data

def f(y): return y*y*(1-y)

y0 = 1e-3
t0, tf = 0, 2/y0
tol = 1e-4
h0 = 1e-2/y0

time loop

t,y,h = t0,y0,h0
sol_t = [t]
sol_y = [y]
sol_th = [t]
sol_hh = [h]
dy=dt=0
while t < tf:
    # compute step and error estimate
    k1 = f(y)            # Euler y+h*k1
    k2 = f(y+h*k1)       # Heun  y+0.5*h*(k1+k2)
    err = abs(k2-k1)*0.5 # unit step error, compares directly to tol
    # compute optimal step size factor
    fac = (tol*(abs(y)+h*abs(k1)))/(tol**3+err)
    fac = max(0.2,min(2,fac))
    # accept or reject step
    if fac >=0.8:
        y += h*k1; t += h
        sol_t.append(t)
        sol_y.append(y)
    # new step size
    h = h*fac
    sol_th.append(t)
    sol_hh.append(h)
    
sol_t = np.asarray(sol_t)
sol_y = np.asarray(sol_y)
sol_th = np.asarray(sol_th)
sol_hh = np.asarray(sol_hh)

plots

The timing in the solution is wrong, the jump should happen slightly before $1/δ=1000$. But the nature of the explicit Euler method makes all curves a little wider. In the zoom part one can observe the expected oscillations around the equilibrium with amplitude about the error tolerance. The corresponding step sizes also almost follow the prediction, oscillating between $1.5$ and $3$, larger steps getting rejected.

Answer 2

One of the many definitions of stiffness is that implicit (A-stable ones in particular) methods solve the system dynamics much more efficiently than explicit ones. In the case of stiff ODEs, an A-stable implicit method is typically always more stable than an explicit method. Such methods (and in particular methods which are also L-stable) dampen the solution modes associated with eigenvalues of large magnitudes. In the case of dissipative systems, these larger eigenvalues typically have a large negative real part, which means that, in the exact solution, they disappear almost instantly. L-stable methods replicate this behaviour for any value of the time step. Then the Time step only needs to be selected based on the slower solution modes. On the opposite, explicit methods must lower the time step to much lower values since all the eigenvalues must be contained in the method's stability domain (which grows linearly in $1/\Delta t$).

A certain exception I can think of is the unsteady heat equation, where certain explicit methods (ROCK4) can be very stable and efficient in the mildly stiff case.

Highly stable implicit methods may however give wrong results if the larger eigenvalues have a large imaginary part and a small negative real part, or a large real part. I. This case, provided the solution can be computed (its uniqueness may no longer be guaranteed as explained in the other answer), the implicit method will dampen the associated modes if the time step is too large, likely removing important dynamics from the numerical solution.

I am curious to see what example has led you to ask this question.