12-13. Introduction to Numerical Methods for Initial-Value Problems

First-Order ODE

Let $I=[x_0,X]\subset \mathbb{R}$ be a time interval
Let $D\subset {\mathbb{R}}^d$ be an open subset where $d\in {\mathbb{N}}^{+}$ denote space dimension

First Order Differential Equation (ODE) is an equation in form

$${\mathbf{y}}^{\prime }(x)=\mathbf{f}(x,\mathbf{y})$$

where $\mathbf{f}:I\times D\to {\mathbb{R}}^d$

Link to original

Initial Value Problem

Let $I=[x_0,X]\subset \mathbb{R}$ be a time interval
Let $D\subset {\mathbb{R}}^d$ be an open subset where $d\in {\mathbb{N}}^{+}$ denote space dimension

Initial Value Problem (IVP) is an ODE with initial condition with

$${\mathbf{y}}^{\prime }(x)=\mathbf{f}(x,\mathbf{y}),\mathbf{y}(x_0)={\mathbf{y}}_0$$

Link to original

Picard's Theorem

Let $\mathbf{f}$ be continuous in neighbourhood $U\subset {\mathbb{R}}^{1+d}$ of $(x_0,{\mathbf{y}}_0)$ containing closed cylinder

$$R=\{(x,\mathbf{y}):x_0\le x\le X_m,||\mathbf{y}-{\mathbf{y}}_0||\le Y_M\}$$

where $X_M>x_0$ and $Y_M>0$ are constants

Suppose exists positive constant $L$ such that

$$||\mathbf{f}(x,\mathbf{y})-\mathbf{f}(x,\mathbf{z})||\le L||\mathbf{y}-\mathbf{z}||$$

holds whenever $(x,\mathbf{y})$ and $(x,\mathbf{z})$ lie in $R$

Let

$$M:=\max \{||\mathbf{f}(x,\mathbf{y})||:(x,\mathbf{y})\in R\}$$

Suppose

$$M(X_m-x_0)\le Y_M$$

Then there exists unique continuously differentiable function

$$[x_0,X_M]\ni nix\mapsto \mathbf{y}(x)\in {\mathbb{R}}^d$$

which is the solution to the IVP

Integral Form for Solution

If solution exists then, solution can be expressed by integral equation in form

$$\mathbf{y}(x)=\mathbf{y}(x_0)+{\int }_{x_0}^x\mathbf{f}(t,\mathbf{y}(t))dx$$

Existence of Solution for Finite Time

Picard’s Theorem only guarantees existence of solution only up to finite time $X_M$

Stability of Solutions with Respect to Initial Conditions

If IVP satisfies assumptions of Picard’s Theorem then

Solution is stable on bounded interval $[x_0,X]$
So if $\mathbf{y}:[x_0,X]\to D$ solves IVP

$${\mathbf{y}}^{\prime }(x)=\mathbf{f}(x,\mathbf{y}),\mathbf{y}(x_0)={\mathbf{y}}_0$$

and $\tilde{\mathbf{y}}:[x_0,X]\to D$ solves same ODE with perturbed initial condition ${\tilde{\mathbf{y}}}_0$ so

$${\tilde{\mathbf{y}}}^{\prime }(x)=\mathbf{f}(x,\tilde{\mathbf{y}}),\tilde{\mathbf{y}}(x_0)=\tilde{\mathbf{y}}0$$

Then

$$|||\mathbf{y}|(x)-\tilde{\mathbf{y}}(x)||\le e^{L(X-x_0)}||{\mathbf{y}}_0-{\tilde{\mathbf{y}}}_0||\forall x\in [x_0,X]$$

Link to original

Discretisation of an IVP

Dividing $[x_0,X]$ into $N\in {\mathbb{N}}^{+}$ subintervals defined by equidistance points

$$x_n=x_0+nh,n=0,\cdots ,N$$

where step size $h$ is defined by

$$h=\frac{X-x_0}{N}$$

with for each step $x_n$ then associate approximation

$${\mathbf{y}}_n\approx \mathbf{y}(x_n)$$

Link to original

Explicit Euler Method

$${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\mathbf{f}(x_n,{\mathbf{y}}_n)$$

Link to original

Implicit Euler Method

$${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\mathbf{f}(x_{n+1},{\mathbf{y}}_{n+1})$$

Link to original

Implicit Midpoint Method

$${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\mathbf{f}\left(x_n+\frac{h}{2},\frac{{\mathbf{y}}_n+{\mathbf{y}}_{n+1}}{2}\right)$$

Link to original

One-Step Method

One-step method is function $\Psi$ which takes
Triplet $(s,\mathbf{y},h)\subset \mathbb{R}\times {\mathbb{R}}^d\times \mathbb{R}$ and function $\mathbf{f}$ to compute approximation

$$\mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})\text{ of }\mathbf{y}(s+h)$$

which is solution at $s+h$ of

$${\mathbf{y}}^{\prime }(x)=\mathbf{f}(x,\mathbf{y}),\mathbf{y}(s)=\mathbf{y}$$

Link to original

Consistent One-Step Method

One-Step Method $\Psi$ is consistent if

$$\mathbf{\Psi }(s,\mathbf{y},0,\mathbf{f})=\mathbf{y}$$

and

$$\frac{d}{dh}\mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f}){\mid }_{h=0}=\mathbf{f}(s,\mathbf{y})$$

Link to original

Consistency Error

Consistency error (aka local truncation error) $\tau$ is defined as

$$\begin{array}{rl}\tau (s,\mathbf{y},h,\mathbf{f})& :=\frac{\mathbf{y}(s+h)-\mathbf{y}}{h}-\frac{\mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})-\mathbf{y}}{h}\\ & =\frac{\mathbf{y}(s+h)-\mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})}{h}\end{array}$$

where $\mathbf{y}(s+h)$ is the solution at $s+h$ of IVP

Link to original

Characterisation of Consistency lemma

Assume $h\mapsto \mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})$ is continuously differentiable in neighbourhood of $0$ then

$\mathbf{\Psi }$ is consistent if and only if for any fixed $\mathbf{f}$

$$||\tau (\tilde{s},\tilde{\mathbf{y}},h,\mathbf{f})||\to 0\text{ as }h\to 0$$

locally uniformly in $(\tilde{s},\tilde{\mathbf{y}})\in R$ where $R$ is the cylinder from Picard’s Theorem

Link to original

Increment Function Form of Consistency lemma

Assume $h\mapsto \mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})$ is continuously differentiable in neighbourhood of $0$ then

$\mathbf{\Psi }$ is consistent if and only if there is continuous incremenent function

$$h\mapsto \psi (s,\mathbf{y},h,\mathbf{f})$$

such that

$$\mathbf{\Psi }(s,\mathbf{y},h,\mathbf{f})=\mathbf{y}+h\psi (s,\mathbf{y},h,\mathbf{f}),\psi (s,\mathbf{y},0,\mathbf{f})=\mathbf{f}(s,\mathbf{y})$$

Link to original

Global Error

Global Error is defined by

$$e_n:=\mathbf{y}(x_n)-{\mathbf{y}}_n$$

Link to original

Global Error Bound for One-Step Methods

Let $\mathbf{\Psi }$ be a consistent one-step method

Assume increment function $\psi$ is Lipschitz continuous with respect to $\mathbf{y}$
So exists positive constant $L_{\psi }$ such that

$$||\psi (x,\mathbf{y},h,\mathbf{f})-\psi (x,\mathbf{z},h,\mathbf{f})||\le L_{\psi }||\mathbf{y}-\mathbf{z}||\text{ for }(x,\mathbf{y}),(x,\mathbf{z})\text{ in }R$$

for all $0\le h\le h_0$ and for the same $R$ from Picard’s Theorem

Then assuming $(x_n,{\mathbf{y}}_n)$ remains $R$ then

$$e\le \left(\frac{\exp (L_{\psi }(x_N-x_0))-1}{L_{\psi }}\right)\underset{n=0,\cdots ,N-1}{\max }||\tau (x_n,{\mathbf{y}}_n,h,\mathbf{f})||$$

Proof

For generic $n\in \{1,\cdots ,N-1\}$

$$\begin{array}{rl}{e}_{n+1}& =||\mathbf{y}(x_{n+1})-{\mathbf{y}}_{n+1}||\\ & =||\mathbf{y}(x_{n+1})-\mathbf{\Psi }(x_n,{\mathbf{y}}_n,h,\mathbf{f})||\\ & =||\mathbf{y}(x_{n+1})-\mathbf{\Psi }(x_n,\mathbf{y}(x_n),h,\mathbf{f})+\mathbf{\Psi }(x_n,\mathbf{y}(x_n),h,\mathbf{f})-\mathbf{\Psi }(x_n,{\mathbf{y}}_n,h,\mathbf{f})||\\ & \le ||\mathbf{y}(x_{n+1})-\mathbf{\Psi }(x_n,\mathbf{y}(x_n),h,\mathbf{f})||+||\mathbf{\Psi }(x_n,\mathbf{y}(x_n),h,\mathbf{f})-\mathbf{\Psi }(x_n,{\mathbf{y}}_n,h,\mathbf{f})||\\ & =h||\tau (x_n,\mathbf{y}(x_n),h,\mathbf{f})||+||(\mathbf{y}(x_n)+h\psi (x,\mathbf{y}(x_n),h,\mathbf{f}))-({\mathbf{y}}_n+h\psi (x,{\mathbf{y}}_n,h,\mathbf{f}))||\\ & \le h||\tau (x_n,\mathbf{y}(x_n),h,\mathbf{f})||+||\mathbf{y}(x_n)-{\mathbf{y}}_n||+h||\psi (x,\mathbf{y}(x_n),h,\mathbf{f})-\psi (x,{\mathbf{y}}_n,h,\mathbf{f})||\\ & =h||\tau (x_n,\mathbf{y}(x_n),h,\mathbf{f})||+e_n+h||\psi (x,\mathbf{y}(x_n),h,\mathbf{f})-\psi (x,{\mathbf{y}}_n,h,\mathbf{f})||\\ & \le h||\tau (x_n,\mathbf{y}(x_n),h,\mathbf{f})||+e_n+h{L}_{\psi }||\mathbf{y}(x_n)-{\mathbf{y}}_n||,\\ & =h||\tau (x_n,\mathbf{y}(x_n),h,\mathbf{f})||+(1+h{L}_{\psi })e_n\end{array}$$

Iterating recursively then implying (with $e_0=0$)

$$\begin{array}{rl}{e}_{n+1}& \le (1+h{L}_{\psi }{)}^{n+1}{e}_0+h\sum _{k=0}^n(1+h{L}_{\psi }{)}^k\underset{m=0,\cdots ,n}{\max }||\tau (x_m,\mathbf{y}(x_m),h,\mathbf{f})||\\ & =\frac{(1+h{L}_{\psi }{)}^{n+1}-1}{L_{\psi }}\underset{m=0,\cdot ,sn}{\max }||\tau (x_m,\mathbf{y}(x_m),h,\mathbf{f})||\end{array}$$

And as $1+h{L}_{\psi }\le \exp h{L}_{\psi }$ then result follows

Link to original