Lipschitz Condition and Convexity

Definition

A function \(f(t,y)\) is said to satisfy a Lipschitz condition in the variable \(y\) on a set \(D\subset \mathbb{R}^2\) if a constant \(L>0\) exists with \[ \begin{aligned} |f(t,y_1)-f(t,y_2)|\le L |y_1-y_2|, \end{aligned} \] whenever \((t,y_1),(t,y_2)\in D\). The constant \(L\) is called a Lipschitz constant for \(f\).

Definition

A set \(D\subset \mathbb{R}^2\) is said to be convex if whenever \((t_1,y_1)\) and \((t_2,y_2)\) belong to \(D\) and \(\lambda\) is in \([0,1]\), the point \(((1-\lambda)t_1+\lambda t_2,(1-\lambda)y_1+\lambda y_2)\) also belongs to \(D\).

Existence and Uniqueness

Theorem

Suppose \(f(t,y)\) is defined on a convex set \(D\subset \mathbb{R}^2\). If a constant \(L>0\) exists with \[ \begin{aligned} \left|\frac{\partial f}{\partial y}(t,y)\right|\le L,\quad \text{for all }(t,y)\in D, \end{aligned} \] then \(f\) satisfies a Lipschitz condition on \(D\) in the variable \(y\) with Lipschitz constant \(L\).

Theorem

Suppose that \(D=\{(t,y)\ |\ a\le t \le b, -\infty < y < \infty \}\) and that \(f(t,y)\) is continuous on \(D\). If \(f\) satisfies a Lipschitz condition on \(D\) in the variable \(y\), then the initial-value problem \[ \begin{aligned} y'(t)=f(t,y),\quad a\le t\le b,\quad y(a)=\alpha, \end{aligned} \] has a unique solution \(y(t)\) for \(a\le t\le b\).

Well-Posedness

Definition

The initial-value problem \[ \begin{aligned} \frac{dy}{dt} = f(t,y),\quad a\le t\le b,\quad y(a)=\alpha, \end{aligned} \] is said to be a well-posed problem if:

• A unique solution, \(y(t)\), to the problem exists, and
• There exist constants \(\varepsilon_0>0\) and \(k>0\) such that for any \(\varepsilon\), with \(\varepsilon_0>\varepsilon>0\), whenever \(\delta(t)\) is continuous with \(|\delta(t)|<\varepsilon\) for all \(t\) in \([a,b]\), and when \(|\delta_0|<\varepsilon\), the initial-value problem \[ \begin{aligned} \frac{dz}{dt}=f(t,z)+\delta(t),\quad a\le t \le b,\quad z(a)=\alpha+\delta_0, \end{aligned} \] has a unique solution \(z(t)\) that satisfies \[ \begin{aligned} |z(t)-y(t)|<k\varepsilon\quad\text{for all }t\text{ in }[a,b]. \end{aligned} \]

Well-Posedness

Theorem

Suppose \(D=\{(t,y)\ |\ a\le t \le b\text{ and }-\infty<y<\infty\}\). If \(f\) is continuous and satisfies a Lipschitz condition in the variable \(y\) on the set \(D\), then the initial-value problem \[ \begin{aligned} \frac{dy}{dt} = f(t,y),\quad a\le t\le b,\quad y(a)=\alpha \end{aligned} \] is well-posed.

Euler’s Method

Suppose a well-posed initial-value problem is given: \[ \begin{aligned} \frac{dy}{dt} = f(t,y),\quad a\le t \le b,\quad y(a)=\alpha \end{aligned} \] Distribute mesh points equally throughout \([a,b]\): \[ \begin{aligned} t_i=a+ih,\quad\text{for each }i=0,1,2,\ldots,N. \end{aligned} \] The step size \(h=(b-a)/N = t_{i+1}-t_i\).

Euler’s Method

Use Taylor’s Theorem for \(y(t)\): \[ \begin{aligned} y(t_{i+1}) = y(t_i) + (t_{i+1}-t_i)y'(t_i) +\frac{(t_{i+1}-t_i)^2}{2}y''(\xi_i) \end{aligned} \] for \(\xi_i\in(t_i,t_{i+1})\). Since \(h=t_{i+1}-t_i\) and \(y'(t_i)=f(t_i,y(t_i))\), \[ \begin{aligned} y(t_{i+1}) = y(t_i) + h f(t_i,y(t_i)) + \frac{h^2}{2}y''(\xi_i). \end{aligned} \] Neglecting the remainder term gives Euler’s method for \(w_i\approx y(t_i)\): \[ \begin{aligned} w_0&=\alpha \\ w_{i+1} &= w_i + h f(t_i,w_i),\qquad i=0,1,\ldots,N-1 \end{aligned} \] The well-posedness implies that \[ \begin{aligned} f(t_i,w_i) \approx y'(t_i) = f(t_i,y(t_i)) \end{aligned} \]

Error Bound

Theorem

Suppose \(f\) is continuous and satisfies a Lipschitz condition with constant \(L\) on \[ \begin{aligned} D = \{(t,y)\ |\ a\le t \le b, -\infty<y<\infty \} \end{aligned} \] and that a constant \(M\) exists with \[ \begin{aligned} |y''(t)|\le M,\quad \text{for all }t\in[a,b]. \end{aligned} \] Let \(y(t)\) denote the unique solution to the initial-value problem \[ \begin{aligned} y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha, \end{aligned} \] and \(w_0,w_1,\ldots,w_n\) as in Euler’s method. Then \[ \begin{aligned} |y(t_i)-w_i|\le \frac{hM}{2L}\left[ e^{L(t_i-a)}-1\right]. \end{aligned} \]

Local Truncation Error

Definition

The difference method \[ \begin{aligned} w_0 &= \alpha \\ w_{i+1} &= w_i + h\phi(t_i,w_i) \end{aligned} \] has local truncation error \[ \begin{aligned} \tau_{i+1}(h) = \frac{y_{i+1}-(y_i+h\phi(t_i,y_i))}{h} = \frac{y_{i+1}-y_i}{h} - \phi(t_i,y_i), \end{aligned} \] for each \(i=0,1,\ldots,N-1\).

Higher-Order Taylor Methods

Consider initial-value problem \[ \begin{aligned} y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha. \end{aligned} \] Expand \(y(t)\) in \(n\)th Taylor polynomial about \(t_i\), evaluated at \(t_{i+1}\): \[ \begin{aligned} y(t_{i+1}) &= y(t_i)+hy'(t_i)+\frac{h^2}{2}y''(t_i) + \cdots \\ &+\frac{h^n}{n!}y^{(n)}(t_i) + \frac{h^{n+1}}{(n+1)!} y^{(n+1)}(\xi_i) \\ &= y(t_i) + hf(t_i,y(t_i)) + \frac{h^2}{2}f'(t_i,y(t_i)) + \cdots \\ &+ \frac{h^n}{n!}f^{(n-1)}(t_i,y(t_i))+ \frac{h^{n+1}}{(n+1)!}f^{(n)}(\xi_i,y(\xi_i)) \end{aligned} \] for some \(\xi_i\in (t_i,t_{i+1})\). Delete remainder term to obtain the Taylor method of order \(n\).

Higher-Order Taylor Methods

Taylor Method of Order \(n\)

\[ \begin{aligned} w_0 &= \alpha \\ w_{i+1} &= w_i + h T^{(n)} (t_i,w_i),\qquad i=0,1,\ldots,N-1 \end{aligned} \] where \[ \begin{aligned} T^{(n)}(t_i,w_i) &=f(t_i,w_i) + \frac{h}{2} f'(t_i,w_i) + \cdots + \frac{h^{(n-1)}}{n!} f^{(n-1)}(t_i,w_i) \end{aligned} \]

Higher-Order Taylor Methods

Theorem

If Taylor’s method of order \(n\) is used to approximate the solution to \[ \begin{aligned} y'(t) = f(t,y(t)),\quad a\le t\le b,\quad y(a)=\alpha, \end{aligned} \] with step size \(h\) and if \(y\in C^{n+1}[a,b]\), then the local truncation error is \(O(h^n)\).

Taylor’s Theorem in Two Variables

Theorem

Suppose \(f(t,y)\) and partial derivatives up to order \(n+1\) continuous on \(D=\{(t,y)\ |\ a\le t \le b, c\le y \le d\}\), let \((t_0,y_0)\in D\). For \((t,y)\in D\), there is \(\xi\in[t,t_0]\) and \(\mu\in[y,y_0]\) with \[ \begin{aligned} f(t,y) &= P_n(t,y)+ R_n(t,y) \\ P_n(t,y) &=f(t_0,y_0) + \left[(t-t_0)\frac{\partial f}{\partial t} (t_0,y_0) + (y-y_0)\frac{\partial f}{\partial y}(t_0,y_0)\right] \\ &+\left[\frac{(t-t_0)^2}{2} \frac{\partial^2 f}{\partial t^2}(t_0,y_0) +(t-t_0)(y-y_0)\frac{\partial^2 f}{\partial t\partial y} (t_0,y_0)\right. \\ &+\left.\frac{(y-y_0)^2}{2} \frac{\partial^2 f}{\partial y^2} (t_0,y_0)\right] + \cdots \\ &+ \left[ \frac{1}{n!} \sum_{j=0}^n \binom{n}{j} (t-t_0)^{n-j} (y-y_0)^j \frac{\partial^n f}{\partial t^{n-j}\partial y^j}(t_0,y_0)\right] \end{aligned} \]

Taylor’s Theorem in Two Variables

Theorem

(cont’d) \[ \begin{aligned} R_n(t,y) &= \frac{1}{(n+1)!} \sum_{j=0}^{n+1} \binom{n+1}{j} (t-t_0)^{n+1-j} (y-y_0)^j \cdot \\ &\cdot\frac{\partial^{n+1} f}{\partial t^{n+1-j}\partial y^j}(\xi,\mu) \end{aligned} \] \(P_n(t,y)\) is the \(n\)th Taylor polynomial in two variables.

Runge-Kutta Methods

• Obtain high-order accuracy of Taylor methods without knowledges of derivatives of \(f\)
• Determine \(a_1,\alpha_1,\beta_1\) such that \[ \begin{aligned} a_1 f(t+\alpha_1,y+\beta_1) \approx f(t,y) + \frac{h}{2} f'(t,y) = T^{(2)}(t,y). \end{aligned} \] with \(O(h^2)\) error.
• Since \[ \begin{aligned} f'(t,y) = \frac{df}{dt}(t,y) = \frac{\partial f}{\partial t}(t,y) +\frac{\partial f}{\partial y}(t,y)\cdot y'(t) \end{aligned} \] and \(y'(t)=f(t,y)\), we have \[ \begin{aligned} T^{(2)}(t,y)=f(t,y)+\frac{h}{2}\frac{\partial f}{\partial t}(t,y) + \frac{h}{2}\frac{\partial f}{\partial y}(t,y)\cdot f(t,y) \end{aligned} \]

Runge-Kutta Methods

• Expand \(f(t+\alpha_1,y+\beta_1)\) in 1st degree Taylor polynomial: \[ \begin{aligned} a_1f(t+\alpha_1,y+\beta_1)&=a_1f(t,y)+a_1\alpha_1\frac{\partial f}{\partial t}(t,y) \\ &+ a_1\beta_1\frac{\partial f}{\partial y}(t,y)+ a_1\cdot R_1(t+\alpha_1,y+\beta_1) \end{aligned} \]
• Matching coefficients gives \[ \begin{aligned} a_1=1 \quad a_1\alpha_1 = \frac{h}{2},\quad a_1\beta_1=\frac{h}{2}f(t,y) \end{aligned} \] with unique solution \[ \begin{aligned} a_1=1,\quad \alpha_1=\frac{h}{2},\quad \beta_1=\frac{h}{2}f(t,y) \end{aligned} \]

Runge-Kutta Methods

• This gives \[ \begin{aligned} T^{(2)}(t,y) = f\left(t+\frac{h}{2},y+\frac{h}{2}f(t,y)\right) -R_1 \left(t+\frac{h}{2},y+\frac{h}{2}f(t,y)\right) \end{aligned} \] with \(R_1(\cdot,\cdot)=O(h^2)\)

Midpoint Method

\[ \begin{aligned} w_0 &= \alpha, \\ w_{i+1} &= w_i + hf\left(t+\frac{h}{2}, w_i+\frac{h}{2}f(t_i,w_i)\right),\qquad i=0,1,\ldots,N-1 \end{aligned} \] Local truncation error of order two.

Runge-Kutta Methods

Runge-Kutta Order Four

\[ \begin{aligned} w_0 &= \alpha \\ k_1 &= hf(t_i,w_i) \\ k_2 &= hf\left(t_i+\frac{h}{2},w_i+\frac{1}{2}k_1\right) \\ k_3 &=hf\left(t_i+\frac{h}{2},w_i+\frac{1}{2}k_2\right) \\ k_4 &= hf(t_{i+1},w_i+k_3) \\ w_{i+1} &= w_i +\frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} \] Local truncation error \(O(h^4)\)

Runge-Kutta Order Four

MATLAB Implementation

function [t, w] = rk4(f, a, b, alpha, N)
% Solve ODE y'(t) = f(t, y(t)) using Runge-Kutta 4.

h = (b-a) / N;
t = (a:h:b)';
w = zeros(N+1, length(alpha));
w(1,:) = alpha(:)';
for i = 1:N
    k1 = h*f(t(i),       w(i,:));
    k2 = h*f(t(i) + h/2, w(i,:) + k1/2);
    k3 = h*f(t(i) + h/2, w(i,:) + k2/2);
    k4 = h*f(t(i) + h,   w(i,:) + k3);
    w(i+1,:) = w(i,:) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

Multistep Methods

Definition

An \(m\)-step multistep method for solving the initial-value problem \[ \begin{aligned} y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha, \end{aligned} \] has a difference equation for approximate \(w_{i+1}\) at \(t_{i+1}\): \[ \begin{aligned} w_{i+1} &= a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0w_{i+1-m} \\ &+ h[b_m f(t_{i+1},w_{i+1}) + b_{m-1} f(t_i,w_i) + \cdots \\ &+ b_0 f(t_{i+1-m},w_{i+1-m})], \end{aligned} \] where \(h=(b-a)/N\), and starting values are specified: \[ \begin{aligned} w_0=\alpha,\quad w_1=\alpha_1,\quad \ldots,\quad w_{m-1}=\alpha_{m-1} \end{aligned} \] Explicit method if \(b_m=0\), implicit method if \(b_m\ne 0\).

Multistep Methods

Fourth-Order Adams-Bashforth Technique

\[ \begin{aligned} w_0 &= \alpha, \quad w_1 = \alpha_1, \quad w_2=\alpha_2, \quad w_3=\alpha_3, \\ w_{i+1} &= w_i + \frac{h}{24}[55 f(t_i,w_i) - 59 f(t_{i-1},w_{i-1}) \\ & \qquad\qquad\ \ \ +37 f(t_{i-2},w_{i-2}) - 9 f(t_{i-3},w_{i-3})] \end{aligned} \]

Fourth-Order Adams-Moulton Technique

\[ \begin{aligned} w_0 &= \alpha, \quad w_1 = \alpha_1, \quad w_2=\alpha_2, \\ w_{i+1} &= w_i + \frac{h}{24}[9 f(t_{i+1},w_{i+1}) + 19 f(t_i,w_i) \\ & \qquad\qquad\ \ \ -5 f(t_{i-1},w_{i-1}) + f(t_{i-2},w_{i-2})] \end{aligned} \]

Derivation of Multistep Methods

Integrate the initial-value problem \[ \begin{aligned} y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha \end{aligned} \] over \([t_i,t_{i+1}]\): \[ \begin{aligned} y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} f(t,y(t))\,dt \end{aligned} \] Replace \(f\) by polynomial \(P(t)\) interpolating \((t_0,w_0),\ldots,(t_i,w_i)\), and approximate \(y(t_i)\approx w_i\): \[ \begin{aligned} y(t_{i+1}) \approx w_i + \int_{t_i}^{t_{i+1}} P(t)\,dt \end{aligned} \]

Derivation of Multistep Methods

Adams-Bashforth explicit \(m\)-step: Newton backward-difference polynomial through \((t_i,f(t_i,y(t_i))),\ldots,(t_{i+1-m},f(t_{i+1-m},y(t_{i+1-m})))\). \[ \begin{aligned} \int_{t_i}^{t_{i+1}} f(t,y(t))\,dt &\approx \int_{t_i}^{t_{i+1}} \sum_{k=0}^{m-1} (-1)^k \binom{-s}{k} \nabla^k f(t_i,y(t_i))\,dt \\ &=\sum_{k=0}^{m-1} \nabla^k f(t_i,y(t_i)) h (-1)^k \int_0^1 \binom{-s}{k}\,ds \end{aligned} \]

\[ \begin{array}{c|cccccc} k & 0 & 1 & 2 & 3 & 4 & 5 \\ \hline (-1)^k\int_0^1 \binom{-s}{k}\,ds & 1 & \frac12 & \frac{5}{12} & \frac{3}{8} & \frac{251}{720} & \frac{95}{288} \end{array} \]

Derivation of Multistep Methods

The method is: \[ \begin{aligned} w_0 &= \alpha,\quad w_1=\alpha_1, \quad w_2=\alpha_2, \\ w_{i+1} &= w_i + \frac{h}{12} [ 23 f(t_i,w_i) -16 f(t_{i-1},w_{i-1}) +5 f(t_{i-2},w_{i-2})] \end{aligned} \]

Local Truncation Error

Definition

If \(y(t)\) solves \[ \begin{aligned} y'=f(t,y),\quad a\le t\le b,\quad y(a)=\alpha, \end{aligned} \] and \[ \begin{aligned} w_{i+1} =& a_{m-1}w_i + \cdots + a_0 w_{i+1-m} \\ & + h[b_m f(t_{i+1},w_{i+1}) + \cdots + b_0 f(t_{i+1-m},w_{i+1-m})], \end{aligned} \] the local truncation error is \[ \begin{aligned} \tau_{i+1}(h) =& \frac{y(t_{i+1})-a_{m-1}y(t_i) - \cdots - a_0 y(t_{i+1-m})} {h} \\ & - [b_m f(t_{i+1},y(t_{i+1})) + \cdots + b_0 f(t_{i+1-m},y(t_{i+1-m}))]. \end{aligned} \]

High-Order Systems of Initial-Value Problems

An \(m\)th-order system of first-order initial-value problems has the form \[ \begin{aligned} \frac{du_1}{dt}(t) &= f_1(t,u_1,u_2,\ldots,u_m), \\ \frac{du_2}{dt}(t) &= f_2(t,u_1,u_2,\ldots,u_m), \\ & \vdots \\ \frac{du_m}{dt}(t) &= f_m(t,u_1,u_2,\ldots,u_m), \end{aligned} \] for \(a\le t \le b\), with the initial conditions \[ \begin{aligned} u_1(a)=\alpha_1, u_2(a)=\alpha_2,\ldots,u_m(a)=\alpha_m. \end{aligned} \]

Existence and Uniqueness

Definition

The function \(f(t,y_1,\ldots,y_m)\), defined on the set \[ \begin{aligned} D=\{(t,u_1,\ldots,u_m)\ |\ a\le t \le b, -\infty<u_i<\infty, i=1,2,\ldots,m\} \end{aligned} \] is said to satisfy a Lipschitz condition on \(D\) in the variables \(u_1,u_2,\ldots,u_m\) if a constant \(L>0\) exists with \[ \begin{aligned} |f(t,u_1,\ldots,u_m)-f(t,z_1,\ldots,z_m)| \le L\sum_{j=1}^m |u_j-z_j|, \end{aligned} \] for all \((t,u_1,\ldots,u_m)\) and \((t,z_1,\ldots,z_m)\) in \(D\).

Existence and Uniqueness

Theorem

Suppose \[ \begin{aligned} D=\{(t,u_1,\ldots,u_m)\ |\ a\le t \le b, -\infty<u_i<\infty, i=1,2,\ldots,m\} \end{aligned} \] and let \(f_i(t,u_1,\ldots,u_m)\), for each \(i=1,2,\ldots,m\), be continuous on \(D\) and satisfy a Lipschitz condition there. The system of first-order differential equations \[ \begin{aligned} \frac{du_k}{dt}(t) = f_k(t,u_1,\ldots,u_m),\quad u_k(a)=\alpha_k,\quad k=1,\ldots,m \end{aligned} \] has a unique solution \(u_1(t),\ldots,u_m(t)\) for \(a\le t \le b\).

Numerical Methods

Numerical methods for systems of first-order differential equations are vector-valued generalizations of methods for single equations.

Fourth order Runge-Kutta for systems

\[ \begin{aligned} \mathbf{w}_0 &= \mathbf{\alpha} \\ \mathbf{k}_1 &= h \mathbf{f}(t_i,\mathbf{w}_i) \\ \mathbf{k}_2 &= h \mathbf{f}(t_i+\frac{h}{2},\mathbf{w}_i+\frac12 \mathbf{k}_1) \\ \mathbf{k}_3 &= h \mathbf{f}(t_i+\frac{h}{2},\mathbf{w}_i+\frac12 \mathbf{k}_2) \\ \mathbf{k}_4 &= h \mathbf{f}(t_{i+1},\mathbf{w}_i+\mathbf{k}_3) \\ \mathbf{w}_{i+1} &= \mathbf{w}_i + \frac16 (\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4) \end{aligned} \] where \(\mathbf{w}_i=(w_{i,1},\ldots,w_{i,m})\) is the vector of unknowns.

Consistency and Convergence

Definition

A one-step difference-equation with local truncation error \(\tau_i(h)\) is said to be consistent if \[ \begin{aligned} \lim_{h\rightarrow 0} \max_{1\le i \le N} |\tau_i(h)| = 0 \end{aligned} \]

Definition

A one-step difference equation is said to be convergent if \[ \begin{aligned} \lim_{h\rightarrow 0} \max_{1\le i \le N} | w_i - y(t_i) | = 0, \end{aligned} \] where \(y_i=y(t_i)\) is the exact solution and \(w_i\) the approximation.

Convergence of One-Step Methods

Theorem

Suppose the initial-value problem \(y'=f(t,y)\), \(a\le t \le b\), \(\quad y(a)=\alpha\) is approximated by a one-step difference method in the form \(w_0=\alpha\), \(w_{i+1}=w_i + h \phi(t_i,w_i,h)\). Suppose also that \(h_0>0\) exists and \(\phi(t,w,h)\) is continuous with a Lipschitz condition in \(w\) with constant \(L\) on \(D\), then: \[ \begin{aligned} D=\{(t,w,h)\ |\ a\le t \le b, -\infty<w<\infty, 0\le h\le h_0 \}. \end{aligned} \]

1. The method is stable;
2. The method is convergent if and only if it is consistent: \[ \begin{aligned} \phi(t,y,0) = f(t,y) \end{aligned} \]
3. If \(\tau\) exists s.t. \(|\tau_i(h)|\le \tau(h)\) when \(0\le h \le h_0\), then \[ \begin{aligned} |y(t_i) - w_i| \le \frac{\tau(h)}{L} e^{L(t_i-a)}. \end{aligned} \]

Root Condition

Definition

Let \(\lambda_1,\ldots,\lambda_m\) denote the roots of the characteristic equation \[ \begin{aligned} P(\lambda) = \lambda^m - a_{m-1}\lambda^{m-1} - \cdots - a_1\lambda - a_0 = 0 \end{aligned} \] associated with the multistep method \[ \begin{aligned} w_{i+1} =& a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0 w_{i+1-m} \\ &+ hF(t_i,h,w_{i+1},w_i,\ldots,w_{i+1-m}). \end{aligned} \] If \(|\lambda_i|\le 1\) and all roots with absolute value 1 are simple, the method is said to satisfy the root condition.

Stability

Definition

1. Methods that satisfy the root condition and have \(\lambda=1\) as the only root of magnitude one are called strongly stable.
2. Methods that satisfy the root condition and have more than one distinct root with magnitude one are called weakly stable.
3. Methods that do not satisfy the root condition are unstable.

Theorem

A multistep method \[ \begin{aligned} w_{i+1}=&a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0 w_{i+1-m} \\ &+h F(t_i,h,w_{i+1},w_i,\ldots,w_{i+1-m}) \end{aligned} \] is stable if and and only if it satisfies the root condition. If it is also consistent, then it is stable if and only if it is convergent.

Stiff Equations

• A stiff differential equation is numerically unstable unless the step size is extremely small
• Large derivatives give error terms that are dominating the solution
• Example: The initial-value problem \[ \begin{aligned} y'=-30y, \quad 0\le t \le 1.5, \quad y(0)=\frac13 \end{aligned} \] has exact solution \(y=\frac13 e^{-30t}\). But RK4 is unstable with step size \(h=0.1\).

Euler’s Method for Test Equation

• Consider the simple test equation \[ \begin{aligned} y'=\lambda y,\quad y(0)=\alpha, \quad \text{where } \lambda<0 \end{aligned} \] with solution \(y(t)=\alpha e^{\lambda t}\).
• Euler’s method gives \(w_0 = \alpha\) and \[ \begin{aligned} w_{j+1} &= w_j + h(\lambda w_j) = (1+h\lambda)w_j = (1+h\lambda)^{j+1}\alpha. \end{aligned} \]
• The absolute error is \[ \begin{aligned} |y(t_j)-w_j| &= \left| e^{jh\lambda} - (1+h\lambda)^j \right| |\alpha| \\ &= \left| (e^{h\lambda})^j - (1+h\lambda)^j \right| |\alpha| \end{aligned} \]
• Stability requires \(|1+h\lambda|<1\), or \(h<2/|\lambda|\).

Multistep Methods

Apply a multistep method to the test equation: \[ \begin{aligned} w_{j+1} =& a_{m-1}w_j + \cdots + a_0 w_{j+1-m} \\ &+h\lambda(b_mw_{j+1} + b_{m-1} w_j + \cdots + b_0 w_{j+1-m}) \end{aligned} \] or \[ \begin{aligned} (1-h\lambda b_m)w_{j+1} - (a_{m-1}+h\lambda b_{m-1}) w_j - \cdots - (a_0 + h\lambda b_0) w_{j+1-m} = 0 \end{aligned} \] Let \(\beta_1,\ldots,\beta_m\) be the zeros of the characteristic polynomial \[ \begin{aligned} Q(z,h\lambda) = (1-h\lambda b_m)z^m - (a_{m-1}+h\lambda b_{m-1}) z^{m-1} -\cdots - (a_0+h\lambda b_0) \end{aligned} \] Then \(c_1,\ldots,c_m\) exist with \[ \begin{aligned} w_j = \sum_{k=1}^m c_k (\beta_k)^j \end{aligned} \] and \(|\beta_k|<1\) is required for stability.

Region of Stability

Definition

The region \(R\) of absolute stability for a one-step method is \(R=\{ h\lambda\in\mathcal{C}\ | \ |Q(h\lambda)|<1\}\), and for a multistep method, it is \(R= \{h\lambda\in\mathcal{C}\ |\ |\beta_k|<1,\text{ for all zeros }\beta_k\text{ of }Q(z,h\lambda)\}\).

A numerical method is said to be A-stable if its region \(R\) of absolute stability contains the entire left half-plane.