MATH661

Lecture 21: Ordinary Differential Equations - Single Step Methods

1.Ordinary differential equations

An $n^{th}$ -order ordinary differential equation given in explicit form

y^{(n)} = f (t, y, y^{'}, \dots, y^{(n - 1)})

(1)

is a statement of equality between the action of two operators. On the left hand side the linear differential operator

ℒ = \frac{d}{d t^{n}}

acts upon a sufficiently smooth function, $y \in C^{(n)} (ℝ)$ , $ℒ : C^{(n)} (ℝ) \to C (ℝ)$ . On the right hand side, a nonlinear operator $ℱ$ acts upon the independent variable $t$ and the first $n - 1$ derivatives

ℱ : ℝ \times C (ℝ) \times \dots \times C^{(n - 1)} (ℝ) .

An associated function $f : ℝ^{n + 1} \to ℝ$ has values given by

f (t) = f (t, y (t), y^{'} (t), \dots, y^{(n - 1)} (t)) .

The numerical solution of (1) seeks to find an approximant of $y$ through:

Approximation of the differentiation operator $ℒ$ ;
Approximation of the nonlinear operator $ℱ;$
Approximation of the equality between the effect of the two operators
$ℒ (y) = ℱ (t, y, \dots, y^{(n - 1)}) .$

These approximation problems shall be considered one-by-one, starting with approximation of $ℒ$ assuming that the action of $ℱ$ is exactly represented through knowledge of $f$ .

Note that an $n^{th}$ -order differential equation can be restated as a system of $n$ first-order equations

𝒛^{'} = 𝑭 (t, 𝒛)

(2)

by introducing

𝒛 = {[\begin{array}{lllll} z_{1} & z_{2} & \dots & z_{n - 1} & z_{n} \end{array}]}^{T} = {[\begin{array}{lllll} y & y^{'} & \dots & y^{(n - 2)} & y^{(n - 1)} \end{array}]}^{T},

𝑭 (t, 𝒛) = {[\begin{array}{lllll} z_{2} (t) & z_{3} (t) & \dots & z_{n} (t) & f (t, z_{1} (t), \dots, z_{n} (t)) \end{array}]}^{T} .

Approximation of the differentiation operator for the problem

y^{'} = f (t, y)

(3)

can readily be extended to the individual equations of system (2).

Construction of approximants to (3) is first considered for the initial value problem (IVP)

y^{'} = f (t, y), y (0) = y_{0} .

(4)

The two procedures are:

Approximation of the differentiation operator;
Differentiation of an approximation of $y$ .

Often the two approaches leads to the same algorithm. The problem (4) has a unique solution over some rectangle $R = [0, T] \times [y_{1}, y_{2}]$ in the $t y$ -plane if $f$ is Lipschitz-continuous, stated as the existence of $K \in ℝ_{+}$ such that

| f (t, y_{2}) - f (t, y_{1}) | ⩽ K | y_{2} - y_{1} | .

Note that Lipschitz continuity is a stronger condition than standard continuity in that it states $| f (t, y_{2}) - f (t, y_{1}) | = 𝒪 (| y_{2} - y_{1} |)$ . Differentiability implies Lipschitz continuity.

Consider approximation of $d / d t$ through forward finite differences

\frac{d}{d t} = \frac{1}{h} (Δ - \frac{1}{2} Δ^{2} + \frac{1}{3} Δ^{3} - \dots -),

(5)

and denote by $y_{i}$ the approximation of $y (t)$ , $y_{i} ≅ y (t_{i})$ at the equidistant sample points $t_{i} = i h$ . Evaluation of (2) with a $k^{th}$ order truncation of (5) then gives

f (t_{i}, y (t_{i})) = (\frac{d y}{d t}) (t_{i}) ≅ \frac{1}{h} (Δ - \frac{1}{2} Δ^{2} + \frac{1}{3} Δ^{3} - \dots - {(- 1)}^{k} \frac{1}{k} Δ^{k}) .

2.Differentiation operator approximation from polynomial interpolants

2.1.Euler forward scheme

For $k = 1$ , the resulting scheme is

\frac{1}{h} Δ y_{i} = \frac{y_{i + 1} - y_{i}}{h} = f (t_{i}, y_{i}) = f_{i} \Rightarrow y_{i + 1} = y_{i} + h f_{i},

where $f_{i} ≅ f (t_{i}, y (t_{i}))$ , and is known as the Euler forward scheme. New values are obtained from previous values. Such methods are said to be explicit schemes. As to be expected from the truncation of (5) to the first term in the series, the scheme is first-order accurate. This can be formally established by evaluation of the error at step $i$

e_{i} = y (t_{i}) - y_{i} .

At the next step, $e_{i + 1} = y (t_{i + 1}) - y_{i}$ , and subtraction of the two errors gives upon Taylor-series expansion

e_{i + 1} - e_{i} = y (t_{i + 1}) - y (t_{i}) - (y_{i + 1} - y_{i}) = y (t_{i}) + h y^{'} (t_{i}) + \frac{h^{2}}{2} y^{''} (ξ_{i}) - y (t_{i}) - h f_{i} .

Since $f_{i} = f (t_{i}, y (t_{i}))$ , the one-step error is given by

τ_{i} = e_{i + 1} - e_{i} = \frac{h^{2}}{2} y^{''} (ξ_{i}) .

After $N$ steps,

e_{N} - e_{0} = \frac{h^{2}}{2} \sum_{i = 1}^{N} y^{''} (ξ_{i}) .

Assuming $e_{0} = 0$ (exact representation of the initial condition),

e_{N} ⩽ \frac{N h^{2}}{2} {|| y^{''} ||}_{\infty} .

Numerical solution of the initial value problem is carried out over some finite interval $[0, T]$ , with $T = N h$ , hence

e_{N} ⩽ h \frac{T}{2} {|| y^{''} ||}_{\infty} = 𝒪 (h),

(6)

indeed with first-order convergence.

Alternatively, one could use the backward or centered finite difference approximations of the derivative

\frac{d}{d t} = \frac{1}{h} (\nabla + \frac{1}{2} \nabla^{2} + \frac{1}{3} \nabla^{3} + \dots) = \frac{1}{h} (δ - \frac{1}{24} δ^{3} + \frac{3}{640} δ^{5} - \dots -) .

(7)

2.2.Backward Euler scheme

Truncation of the backward operator at first order gives

f (t_{i}, y (t_{i})) = {(\frac{d y}{d t})}_{i} ≅ \frac{1}{h} {(\nabla y)}_{i} = \frac{y_{i} - y_{i - 1}}{h} \Rightarrow y_{i} = y_{i - 1} + h f_{i} = y_{i - 1} + h f (t_{i}, y_{i}) .

Note now that the unknown value $y_{i}$ appears as an argument to $f$ , with $f_{i} = f (t_{i}, y_{i})$ , the approximation of the exact slope $f (t_{i}, y (t_{i}))$ . Some procedure to solve the equation

y_{i} - y_{i - 1} - h f (t_{i}, y_{i}) = 0,

must be introduced in order to advance the solution from $t_{i - 1}$ to $t_{i}$ . Such methods are said to be implicit schemes. The same type of error analysis as in the forward Euler case again leads to the conclusion that the one-step error is $𝒪 (h^{2})$ , while the overall error over a finite interval $[0, T]$ satisfies (6), and is first-order.

2.3.Leapfrog scheme

Truncation of the centered operator at first order gives

f (t_{i}, y (t_{i})) = {(\frac{d y}{d t})}_{i} ≅ \frac{1}{h} {(δ y)}_{i} = \frac{y_{i + 1 / 2} - y_{i - 1 / 2}}{h} \Rightarrow y_{i + 1 / 2} = y_{i - 1 / 2} + h f_{i} = y_{i - 1 / 2} + h f (t_{i}, y_{i}) .

The higher-order accuracy of the centered finite differences leads to a more accurate numerical solution of the problem (4). The one-step error is third-order accurate,

e_{i + 1 / 2} - e_{i - 1 / 2} = y (t_{i + 1 / 2}) - y (t_{i - 1 / 2}) + h f (t_{i}, y_{i}) = \frac{h^{3}}{3} y^{'''} (ξ_{i}),

and the overall error over interval $[0, T = N h]$ is second-order accurate

e_{N} ⩽ \frac{h^{2}}{3} T {|| y^{'''} ||}_{\infty} .