Numerical solution of Differential Equation

Reference

微分方程数值解：有限差分理论方法与数值计算, 张文生

微分方程数值解, 陈文斌

In ordinary differential equation

\[ \frac{du}{dt}f(t,u),\quad u(0)=u_0 \]

in region \([a,b]\). Note that it is hard to know the function, but we know its derivatives at all points on the region. That is, for any dicrete point \(t_n\in [a,b]\), we have

\[ \frac{du}{dt}\Bigg|_{t=t_n}=f(t_n,u(t_n)) \]

The following part is focused on using there infomation to get a solution of ODE.

数值微分 | Numerical Differentiation

Denote \(h\) as the time step forward, then according to Taylor expansion evaluated at \(t_n\) and take another forward point \(t_{n+1}\)

\[ \begin{equation} u(t_{n+1})=u(t_n)+h u'(t_n)+\frac{h^2}{2!}u''(t_n)+\frac{h^3}{3!}u'''(t_n)+o(h^4) \label{forward} \end{equation} \]

backward at point \(t_{n-1}\), we substitute \(h\) with \(-h\) and get

\[ \begin{equation} u(t_{n-1})=u(t_n)-h u'(t_n)+\frac{h^2}{2!}u''(t_n)-\frac{h^3}{3!}u'''(t_n)+o(h^4)\label{backward} \end{equation} \]

So we have an approximation of the first order derivative from equation \(\ref{forward}\), that is, the forward form of divided difference

\[ \begin{align} u'(t_n)=\frac{u(t_{n+1})-u(t_n)}{h}-\frac{h}{2!}u''(t_n)-\frac{h^2}{3!}u'''(t_n)+o(h^3) \label{Forward} \end{align} \]

or from equation \(\ref{backward}\), that is, the backward form of divided difference

\[ \begin{equation} u'(t_n)=\frac{u(t_{n})-u(t_{n-1})}{h}+\frac{h}{2!}u''(t_n)-\frac{h^2}{3!}u'''(t_n)+o(h^3) \label{Backward} \end{equation} \]

Combine equation \(\ref{Forward}\) and equation \(\ref{Backward}\) we have the center divided difference form

\[ \begin{equation} u'(t_n)=\frac{u(t_{n+1})-u(t_{n-1})}{2h}-\frac{h^2}{3!}u'''(t_n)+o(h^3) \label{center} \end{equation} \]

which is more accurate. Now we can use the above equation to get Euler' method.

Three Point Midpoint Formula

\[ f'(x_0) = \frac{1}{2h}[f(x_0+h)-f(x_0-h)] -\frac{h^2}{6}f'''(\xi) \]

Use Lagrange Polynomial to approximate a function. Then the \(n\)th derivative of Lagrange Poly approximate the \(n\)th derivative of the original function.

Euler's Method

Here we introduce Euler's one step method

\[ \begin{cases} \omega_0=\alpha,\\ \omega_{i+1}=\omega_i+hf(t_i,\omega_i),\quad i=0,1,\cdots,n-1 \end{cases} \]

From forward form of divided difference \(\ref{Forward}\) we have its rounding error

\[ \tau_n=\frac{w(t_{n+1})-w(t_n)}{h}=\frac{h}{2}f'(\eta),\quad t_n\leq \eta\leq t_{n+1} \]

First analyze its convergence.

Convergence of Euler's Method

Assume solution \(\phi(t)\in C^2[t_0,b]\), then the error of the solution \(w_n\) acquired from Euler's method satisfies

\[ \max_{a\leq t_n \leq b}|w(t_n)-w_n|\leq e^{(b-t_0)L}|e_0|+\frac{e^{(b-t_0)L}-1}{L}\tau(h) \]

where

\[ \tau(h)=\frac{h}{2}\|w''\|_{\infty},\quad e_0=w(t_0)-w_0 \]

If \(h\rightarrow 0\), s.t.

\[ |w(t_0)-w_0|\leq c_1h,\quad c_1>0 \]

then \(\exists c>0\), such that

\[ \max_{t_0\leq t_n\leq b}|w(t_n)-w_n|\leq ch \]

Then we analyze its asymptotic stability.

Definition of asymptotic stability

We call Euler's one step method is asymptotic stable, if \(\exists h_0>0\), \(\exists C>0\), such that \(\forall h\in (0,h_0]\),

\[ |z_n-w_n|\leq C\varepsilon,\quad \forall 0\leq n\leq N(h) \]

where \(w_n\) and \(z_n\) denotes the solution from before disturbance(initial value \(w_0\)) and after disturbance(initial value \(z_0\)) using Euler's one step method.

The above shows that Euler's method is asymptotic stable.

Implicite Euler's method

Implicite means we have to guess the number before iterating. Usually it is used after Explicite Euler's method to improve stability.

\[ \begin{cases} \omega_0=\alpha, \\ \omega_{i+1}=\omega_i+hf(t_i,\omega_{i+1}),\quad i=0,1,\cdots,n-1 \end{cases} \]

From Backward form of Numerical Derivative(Divided Difference) \(\ref{Backward}\), we can express its rounding error

\[ \tau_n=\frac{w(t_{n+1})-w(t_{n})}{h}=-\frac{h}{2}f'(\eta),\quad t_n\leq \eta\leq t_{n+1} \]

Modified Euler's Method

So we can choose higher order items to improve accuracy.

This method is also called Tranpezoid Method.

\[ w(t_{n+1})=w(t_n)+\frac{h}{2}\left[f(t_n,w(t_n))+f(t_{n+1},w(t_{n+1}))\right] \]

We have this inspiration from tranpezoidal integration.

From Tranpezoidal integration

From integration the ODE on \([t_n, t_{n+1}]\) we have

\[ w(t_{n+1})=w(t_n)+\int_{t_n}^{t_{n+1}}f(\tau, w(\tau))d\tau \]

Using Tranpezoidal Integration Formula with its Error, we have

\[ w(t_{n+1})=w(t_n)+\frac{h}{2}\left[f(t_n,w(t_n))+f(t_{n+1},w(t_{n+1}))\right]-\frac{h^3}{12}f'''(\eta_n),\quad t_n\leq \eta_n\leq t_{n+1} \]

which means a higher accuracy for method

\[ w(t_{n+1})=w(t_n)+\frac{h}{2}\left[f(t_n,w(t_n))+f(t_{n+1},w(t_{n+1}))\right] \]

Double-step Method

We have this method inspired by the center form of divided difference.

According to center form of divided difference \(\ref{center}\) we have

\[ \begin{cases} \omega_0=\alpha,\\ \omega_1=\omega_0+hf(t_0,\omega_0), \\ \omega_{i+1}=\omega_{i-1}+2hf(t_i,\omega_i), \quad i=1,2,\cdots,N-1 \end{cases} \]

with its rounding error expressed by

\[ \tau_n=\frac{w_{t_{n+1}}-w(t_n)}{h}=-\frac{h^2}{3!}f''(\eta), \quad t_n\leq \eta\leq t_{n+1} \]

Taylor method of order n

To improve precision, we can choose more items of Taylor's expansion for the one order derivative.

\[ \begin{cases} \omega_0=\alpha, \\ \omega_{i+1}=\omega_i+hT^{(n)}(t_i,\omega_i),\quad i=0,1,\cdots,N-1 \end{cases} \]

where

\[ T^{(n)}(t_i,\omega_i)=f(t_i,\omega_i)+\frac{h}{2}f'(t_i,\omega_i)+\cdots+\frac{h^{n-1}}{n!}f^{(n-1)}(t_i,\omega_i) \]

Euler's method is Tarlor's method of order one.

Runge-Kutta Method

For Taylor's method of higher order, it is trivial to calculate high order derivatives, especially when \(p>3\) and \(f\) has a complex expression. So the following Runge-Kutta Method aims to calculate function \(f\)'s values at different points to avoid high order derivatives, reducing calculation time.

• Runge-Kutta method with order 4.

\[ \begin{cases} w_{i+1}=w_i+\frac{1}{6}(k_1+2K_2+2k_3+k_4) w_0=\alpha \end{cases} \]

Multistep Methods

Derive from integration

\[ y(t_{i+1})-y(t_i)=\int_{t_i}^{t_{i+1}}f(t,y)dt \]

use interpolation polynomial to replace f(t,y).

• Adams-Bashforth Four-step Explicit Method

\[ w_{i+1}=w_i + \frac{h}{24}(55f_i-59f_{i-1}+37f_{i-2}-9f_{i -3}) \]

• Adams-Bashforth Three-step Implicit Method

\[ w_{i+1}=w_i +\frac{h}{24}(9f_{i+1}+19f_{i}-5f_{i-1}+f_{i-2}) \]

Or we derive from Taylar expansion.

ODEs