Implicit methods and stiffness

Implicit methods and stiffness#

Up until this point, we have been discussing explicit solver schemes, in which \(y_{i+1}\) is determined using the information \(x_i\) and \(y_i\) with kowledge of \(f(x_i, y_i)\).

Implicit methods uses information at the new time step, in order to determine it; i.e.: using \(f(x_{i+1}, y_{i+1})\)! The new value of \(y_{i+1}\) is therefore calcualted implicitly, for which we will generally need to use a root finder. The root finder adds significant computational expense but is substantially less sensative to numerical instability.

The Backward Euler method#

The simplest implicit scheme is the implicit Euler method (backward Euler). Like the forward Euler method we assumes a constant slope over the timestep but this time it is the slope at the end of the timestep:

\[ y_{i+1} = y_i + f(x_{i+1}, y_{i+1}) h\]

which implicitly defines \(y_{i+1}\).

Numerical instability#

An algorithm is numerically unstable when small errors that occur during computation grow. Such errors can occur due to user-choices (e.g. - too large a step size) or simply round-off/truncation error.

NB: numerical instability is a property of the algorithm, not the equation.

Example: Initial value problem#

Integrate

\[\frac{\partial y}{\partial t} = - y\]

from 0 to 20 using the forward and backward Euler methods and compare to the exact solution,

\[ y(t) = e^{-t}\]

for varying step sizes.

Interactive Comparison: Explicit vs. Implicit Euler#

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider


def plot_with_slider(h):
  t0 =0
  t_end = 20
  y0 = 2
  # Time points
  t_points = np.arange(t0, t_end + h, h)

  # Analytical solution
  y_analytical = y0 * np.exp(-t_points)

  # Explicit Euler method
  y_explicit = np.zeros_like(t_points)
  y_explicit[0] = y0
  for i in range(1, len(t_points)):
      y_explicit[i] = y_explicit[i-1] * (1 - h)

  # Implicit Euler method
  y_implicit = np.zeros_like(t_points)
  y_implicit[0] = y0
  for i in range(1, len(t_points)):
      y_implicit[i] = y_implicit[i-1] / (1 + h)

  # Plotting the results
  plt.figure(figsize=(10, 6))
  plt.plot(t_points, y_analytical, label='Analytical Solution', color='black', linestyle='--')
  plt.plot(t_points, y_explicit, label='Explicit Euler', color='blue', marker='o')
  plt.plot(t_points, y_implicit, label='Implicit Euler', color='red', marker='x')
  plt.xlabel('Time')
  plt.ylabel('y(t)')
  plt.xlim(0, 20)
  plt.title('Comparison of Explicit and Implicit Euler Methods')
  plt.legend()
  plt.grid(True)
  plt.show()

interact(plot_with_slider, h=FloatSlider(min=.1, max=5, step=.1, value=0.8));

Gadzooks!

Let’s consider the Forward Euler scheme analytically. For generality, say \(\frac{dy}{dt} = -ay\) for some positive constant \(a\),

\[\begin{split}\begin{align} y_{i+1} &= y_i+h \frac{dy}{dt}(y_{i}) \\ &= y_i - y_i a h \\ &= y_i [1-a h] \end{align} \end{split}\]

In order for the solution to follow the smoothly decreasing behaviour, \(h\lt\frac{1}{a}\). If \(\frac{1}{a}\lt h \lt \frac{2}{a}\), the solution will overshoot but eventually correct itself. But if \(h \gt \frac{2}{a}\) the method will catasrophically overcorrect!

Compare this to the implicit case,

\[\begin{split}\begin{align} y_{i+1} &= y_i+h \frac{dy}{dt}(y_{i+1}) \\ &= y_i - y_{i+1} a h \\ y_{i+1}&= \frac{y_i}{1+ah} \end{align} \end{split}\]

which is nicely decreasing for all \(h>0\) and is therefore unconditionally stable!

Numerical stability has therefore placed a limit on the step size for which a method will converge on the correct answer. There are different types of stability, which is outside the scope of this course.

The stability of the explicit vs implicit Euler methods helps show some intuition on what’s going on: since implicit schemes focus on the end of the time step and how we got there, they are less sensative to errors introduced by large time steps or numerical noise.

Stiffness#

One might think that the adaptive time steppers may dynamically avoid numerical instability but you would be dissappointed…

Consider the equation

\[ \frac{dy}{dt} = -1000 y + 3000 - 2000 e^{-t} \]

with

\[ y(0)=0 \]

The analytical solution is

\[ y(t) = 3-0.998 e^{-1000 t} - 2.0-2 e^{-t} \]

which features a fasts transient due to \(e^{-1000 t}\) followed by a slow progression:

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider

def plot_function(limit):
  t = np.linspace(0, limit, 500)
  y = 3 - 0.998 * np.exp(-1000 * t) - 2 - 2 * np.exp(-t)
  plt.figure(figsize=(10, 6))
  plt.plot(t, y)
  plt.xlabel('t')
  plt.ylabel('y(t)')
  plt.title('Plot of y(t) = 3 - 0.998e^(-1000t) - 2 - 2e^(-t)')
  plt.grid(True)
  plt.show()

interact(plot_function, limit=FloatSlider(min=.1, max=5, step=0.1, value=3));

The stability of the fast term with forward Euler requires \(h\lt \frac{2}{1000}\) to not catastrophically overcorrect, but at least that is only for the first little bit right?

Wrong. Even though the transient only dominates the behaviour for the first little bit, it is still there for the rest, and will still tear shirt up!

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider


def plot_with_slider(h):
  t0 =0.01
  t_end = .05
  y0 = 0.
  # Time points
  t_points = np.arange(t0, t_end + h, h)

  # Analytical solution
  def y_true(t):
    return 3. - 0.998*np.exp(-1000*t)  - 2.002*np.exp(-t)
  y_analytical = y_true(t_points)

  y0 = y_true(t0)

  def f(y,t):
    return -1000.*y + 3000. - 2000.*np.exp(-t)

  # Explicit Euler method
  y_explicit = np.zeros_like(t_points)
  y_explicit[0] = y0
  for i in range(1, len(t_points)):
      y_explicit[i] = y_explicit[i-1] + f(y_explicit[i-1], t_points[i-1]) * h

  # Implicit Euler method
  y_implicit = np.zeros_like(t_points)
  y_implicit[0] = y0
  def imp(y, t, y0, f, h):
    return y - (y0 + h*f(y, t))

  from scipy.optimize import root
  for i in range(1, len(t_points)):
      sol = root(imp, y_implicit[i-1], args=(t_points[i], y_implicit[i-1], f, h))
      y_implicit[i] = sol.x[0]

  # Plotting the results
  plt.figure(figsize=(10, 6))
  plt.plot(t_points, y_analytical, label='Analytical Solution', color='black', linestyle='--')
  plt.plot(t_points, y_explicit, label='Explicit Euler', color='blue', marker='o')
  plt.plot(t_points, y_implicit, label='Implicit Euler', color='red', marker='x')
  plt.xlabel('Time')
  plt.ylabel('y(t)')
  plt.xlim(t0, t_end)
  plt.title('Comparison of Explicit and Implicit Euler Methods')
  plt.legend()
  plt.grid(True)
  plt.show()

interact(plot_with_slider, h=FloatSlider(min=.0001, max=.003, step=.0001, value=0.0005, readout_format='.4f'));

Definition of stiffness#

A precise mathematical definition of stiffness is difficult, but in general: A stiff equation is one in which accuracy of the solution is determined by numerical stability rather than the behaviour of the solution.

This is expecially problematic in physics where short transients as a system reaches a ‘well behaved state’ are common. If we include the phenomena which describe the transients, we are bound by their timescales weather or not the transient has died off!

One approach is to ignore the transients entirely! Examples include:

  • Quasistatic elasticity

  • diffusive heat / mass transport

  • etc.

Alternately, we can use implicit methods which are more stable, but come with added computational expense.

Contents