Differential Equations

Say we have $f:[a,b]\to \mathbb{R}$, how do we approximate $f^{\prime }(x)$?

Of course we need Taylor’s theorem! Suppose $f\in C^{n_1}([a,b])$, for any $x,y\in [a,b]$ such that $x+y\in [a,b)$ we have $f(x+y)=f(x)+f^{\prime }(x)y+\frac{f^{\prime \prime }(x)}{2}{y}^2+\cdots +\frac{f^{(n)}(x)}{n!}{y}^n+\frac{f^{(n+1)}(\xi (x,y))}{(n+1)!}{y}^{n+1}$ for $x<\xi (x,y)<x+y$

Now assume $f\in C^2([a,b])$, $x\in [a,b]$, $x+h\in [a,b]$, $h>0$
$f(x+h)=f(x)+f^{\prime }(x)h+\frac{f^{\prime \prime }(\xi )}{2}{h}^2$, $x<\xi <x+h$

We can rewrite this to get a familiar formula,
$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}+\frac{f^{\prime \prime }(\xi )}{2}h$

In our methodology, we view this as an approximation that comes with an error term, and we arrive at the forward finite difference equation
$|f^{\prime }(x)-\frac{f(x+h)-f(x)}{h}|=\frac{|f^{\prime \prime }(\xi )|}{2}h\le \frac{h}{2}{M}_2$ where $M_2={\max }_{z\in [x,x+h]}|f^{\prime \prime }(z)|$

We can also write backward finite difference
$|f^{\prime }(x)-\frac{f(x)-f(x-h)}{h}|=\frac{f^{\prime \prime }(\xi )}{2}h\le \frac{h}{2}{M}_2$ where $M_2={\min }_{x\in [x,x-h]}|f^{\prime \prime }(z)|$

Can we do better than an error of $O(h)$? We can if we know more about the derivatives

If we assume $f\in C^3([a,b])$, we can combine higher order forward and backward equations to yield,
$|f^{\prime }(x)-\frac{f(x+h)-f(x-h)}{2h}|=\frac{h^2}{6}|f^{\prime \prime \prime }({\xi }_{+})-f^{\prime \prime \prime }({\xi }_{-})|\le \frac{h^2}{3}{M}_2$ where $M_2={\max }_{z\in [x-h,x+h]}|f^{\prime \prime \prime }(z)|$

The key point is this procedure has essentially equivalent complexity but gives us $O(h^2)$, much faster convergence, at the cost of requiring more regularity

Now what if we’d like to calculate the integral $I={\int }_a^{b}f(a)dx$? We use trapezoids

$T_1=\frac{f(a)+f(b)}{2}(b-a)$
The line for this is $p_1(x)=\frac{b-x}{b-a}f(a)+\frac{x-a}{b-a}f(b)$
$E_1\le {\int }_a^b|f(x)-p_1(x)|dx$
$\le \frac{M_2}{2}{\int }_a^b|(x-a)(x-b)|dx$ where $M_2={\max }_{z\in [a,b]}|f^{\prime \prime }(z)|$ (see Lagrange interpolation)
$=\frac{(b-a{)}^3}{12}{M}_2$

We can divide the function into strips to improve this approximation,
$I={\int }_a^{b}f(x)dx={\sum }_{j=1}^m{\int }_{x_{j-1}}^{x_j}f(x)dx$
$\approx {\sum }_{j=1}^m\frac{f(x_{j-1})+f(x_j)}{2}h$,
$=\frac{h}{2}[f(x_0)+f(x_m)]+h{\sum }_{j=1}^{m-1}f(x_j)$
$=T_m$
where $x_j=a+\frac{b-a}{m}j$ and $h=\frac{b-a}{m}$

We can view $I$ as an inner product between values of $f$ and weights, $w^{T}f$

Observe that this is pretty similar to the equations we derived for differentiation; this is ultimately because they are both linear operators

ODEs

Overview

Now say we have a first-order ODE with an initial condition,
$y^{\prime }(t)=f(t,y(t))$
$y(t_0)=y_0$

To guarantee a well-behaved solution, we assume Lipschitz continuity $|f(t,y_1)-f(t,y_2)|\le L|y_1-y_2|$

Say we are given a physical system, such as the mass-spring oscillator,
$m{x}^{\prime \prime }=-kx+g(t)$

We can convert this into a system of first-order equations,
$x^{\prime }(t)=v(t)$
$v^{\prime }(t)=-(k/m)x(t)+g(t)/m$

Expressed in vector form,
$y^{\prime }(t)=\left[\begin{array}{c}v(t)\\ -(k/m)x+g(t)/m\end{array}\right]$ where $y(t)=\left[\begin{array}{c}x(t)\\ v(t)\end{array}\right]$

Most ODEs don’t have closed-form solutions, so numerical methods are very important here

Our general idea is to create a series of approximations at evenly spaced time steps, but there are multiple techniques to do this

Forward Euler uses the forward finite difference formula for the derivative,
$y^{\prime }(t_n)\approx \frac{y(t_{n+1})-y(t_n)}{h}+O(h)$
$y(t_{n+1})\approx y(t_n)+h\cdot f(t_n,y(t_n))$

This is as simple as it gets, but requires a “small enough” step size

Backward Euler uses the backward finite difference formula,
$y^{\prime }(t_n)\approx \frac{y(t_n)-y(t_{n-1})}{h}+O(h)$
$y(t_n)\approx y(t_{n-1})+h\cdot f(t_n,y(t_n))$

This cannot be solved explicitly because we have a potentially complex equation in terms of $y(t_n)$

Instead, we must solve for $y(t_n)$ using a root-finding method,
$y(t_n)-y(t_{n-1})-h\cdot f(t_n,y(t_n))=0$

This method is more complex but stable for any step size

Crank-Nicolson uses the trapezoidal rule for integration,
${\int }_{t_{n-1}}^{t_n}{y}^{\prime }(t)dt={\int }_{t_{n-1}}^{t_n}f(t,y(t)dt$
$y(t_n)-y(t_{n-1})\approx \frac{h}{2}[f(t_n,y(t_n))-f(t_{n-1},y(t_{n-1}))]+O(h^3)$
$y(t_n)\approx y(t_{n-1})+\frac{h}{2}[f(t_n,y(t_n))-f(t_{n-1},y(t_{n-1}))]$

Again, this is an implicit equation that must be solved with a root-finding method

This method is unconditionally stable and ends up with an error of $O(h^2)$!

Theory

How do we actually analyze the errors? Let’s be more precise

Definition: A scheme is convergent if $|e_n|\le \psi (h)$, $\psi (h)\to 0$ as $h\to 0$
A scheme converges with order q if $\psi (h)=O(h^q)$ as $h\to 0$

$e_n=y_n-U_n$ is the error between the exact value $y_n$ and the approximated value $U_n$

In practice, we split our error

$e_n=y_n-U_n=(y_n-U_n^{\ast })+(U_n^{\ast }-U_n)$ where $U_n^{\ast }$ is the approximation made when $U_{n-1}$ is exact (the ideal estimate)

$|e_n|\le |E_{1,n}|+|E_2,n|$
$E_{1,n}=y_n-U_n^{\ast }$ measures local consistency
$E_{2,n}=U_n^{\ast }-U_n$ measures global consistency

Definition: A scheme is consistent if $\tau (h)\to 0$ as $h\to 0$ where $\tau (h)={\max }_{n=1}^N\left|\frac{E_{1,n}}{h}\right|$ is called the global truncation error (the maximum of the local truncations)

Let’s put this into practice

Theorem: If $y(t)$ of the Cauchy problem is in $C^2([t_0,T])$ then Forward-Euler is convergent with order 1

$y(t_n)=y(t_{n-1})+h{y}^{\prime }(t_{n-1})+\frac{h^2}{2}{y}^{\prime \prime }({\xi }_{n-1})$ for some ${\xi }_{n-1}\in (t_{n-1},t_n)$
$|E_{1,n}|=\frac{h^2}{2}|y^{\prime \prime }({\xi }_{n-1})|\le M_2\frac{h^2}{2}$, where $M_2={\max }_t|y^{\prime \prime }(t)|$

So $\tau (h)\le M_2\frac{h}{2}$

$E_{2,n}=y_{n-1}-U_{n-1}+h[f(t_{n-1},y_{n-1})-f(t_{n-1},U_{n-1})]$
$|E_{2,n}|\le |e_{n-1}|+h|f(t_{n-1},y_{n-1})-f(t_{n-1},U_{n-1})|$
$\le (1+hL)|e_{n-1}|$, using the Lipschitz continuity condition

We combine $E_{1,n}$ and $E_{2,n}$

$|e_n|\le O(h^2)+(1+hL)|e_{n-1}|$
$\le O(h^2)\cdot [1+(1+hL)+(1+hL{)}^2+\cdots +(1+hL{)}^{n-1}]$
$=O(h^2)\cdot \frac{(1+hL{)}^n-1}{hL}=O(h^2)\cdot \frac{(1+hL{)}^{T/h}-1}{hL}\le O(h^2)\cdot \frac{e^{LT}-1}{hL}$
$=O(h)$

If we run through similar steps with the other schemes, we find Backward Euler is also order 1 and Crank-Nicolson is order 2

We can go deeper with our analysis by analyzing stability, which quantifies how a scheme responds to perturbations

Definition: A scheme is zero-stable if there exists $C,{ϵ}_0,h_0>0$ such that if $|p_n|<ϵ$ for all $n$ then $|U_n-z_n|<C_{ϵ}$, where $ϵ<{ϵ}_0$, $h<h_0$, and $z_n$ is the estimate given a perturbation $p_n$ to $U_{n-1}$

Theorem: A consistent scheme is convergent if and only if it is zero-stable

Definition: A scheme is absolutely stable if ${\lim }_{n\to \infty }|U_n|=0$ for the simple system $y^{\prime }(t)=\lambda y(t)$, $y(0)=1$

For Forward Euler, $U_n=U_{n-1}+hf(t_{n-1},U_{n-1})=(1+h\lambda {)}^n{U}_0$
So it is stable if $|1+h\lambda |<1$, or $h<\frac{2}{|\lambda |}$, which we call conditional absolute stability

Backwards Euler leads to the equation $U_n={\left(\frac{1}{1-h\lambda }\right)}^n{U}_0$, which is unconditional

Crank-Nicolson leads to the equation $U_n={\left(\frac{1+h\lambda /2}{1-h\lambda /2}\right)}^n{U}_0$, which is also unconditionally true

This kind of analysis can be used for more complex equations too