常微分方程初值问题: Euler 法

常微分方程初值问题又称为柯西问题, 在物理、化学、金融、社会科学等领域具有重要应用. Euler 方法是求解一阶常微分方程最基础的经典数值方法, 本文给出了 Euler 对 Euler 方法的描述和误差估计, 最后给出了一些数值结果.

前言

我们考虑如下的一阶常微分方程 $$ \begin{equation} \tag{Eq. 1} \begin{split} y'(t) &= f(t,y(t)),\ t\in (t_0,T]\\ y(t_0) &= y_0. \end{split} \end{equation} $$ 1768 年出版的 Euler 三卷本著作《Institutiones Calculi Integralis》中收录了著名的求解常微分方程的 Euler 方法 — 一种显式的单步方法. 它基于一个非常简单的数学思想: 当时间间隔非常短时, 短到运动物体的速度可以看成是恒定的, 那么这段时间的位移就等于当前的速度乘以时间间隔: $$ \Delta y = v \cdot \Delta t, \quad \text{i.e.}, \quad y(t_0+\Delta t) – y(t_0) = y'(t_0) \cdot \Delta t. $$ 如果我们将上面的公式进行简单的变形, \(y'(t_0)=\frac{y(t_0+\Delta t) – y(t_0)}{\Delta t} \), 可以看出这就是一阶导数的近似. 其实, 数值分析中的很多近似方法都植根于数学对象的数学定义, “真理极其复杂以至于除了近似之外什么都不允许. — John von Neumann”.

Euler 对 Euler 法的描述

了解 Euler 的描述可以加强我们对 Euler 方法的理解. 最初的 Euler 方法求解的问题实际上就是形如 (Eq. 1) 的带有初值的微分方程, 不过在当时, 数学家更喜欢使用另外一种等价形式, $$ y = \int f(t,y(t)) dt . $$ Euler 的目标是找到当变量 \(t\) 从 \(t_0\) 微小变化到 \(t_0+\Delta t\) 时对应的 \(y(t_0+\Delta t) \) 的近似值. Euler 假设 \(f\) 在这段微小的时间间隔内是常量, 即 \(y'(t) \equiv F := f(t_0, y_0) \). 从而, 存在一个常数 \(C\) 使得 $$ y(t)=F\cdot t+ C,\quad t\in (t_0,t_0+\Delta t). $$ 因为初值已知 \(y(t_0)=y_0\), 代入上面的等式可得 \(C = y_0-F \cdot t_0\). 因此, 可以得到近似值 $$ y(t_0+\Delta t) = y_0 + F\cdot \Delta t,\quad \forall t\in (t_0,t_0+\Delta t). \tag{Eq. 2} $$ 上面这个迭代过程其实就是现代数值分析教材中给出的 Euler 方法 \(y_{k+1}=y_k + \Delta t f(x_k, y_k) \). 如果继续上面的过程, 将 \(t_0+\Delta t\) 看作新的 \(t_0\), \(y(t_0+\Delta t) \) 看作新的 \(y_0\), 就可以迭代得到任意点的近似值.

导出 (Eq. 2) 的另一种方式. 通过对等式 \(y'(t) \equiv F\) 在微小区间 \((t_0, t_0+\Delta t) \) 上进行积分可以得到
\(y(t_0+\Delta t)-y_0 = \int_{t_0}^{t_0+\Delta t}F dt=F\cdot \Delta t\), 即 (Eq. 2).

Euler 法的定义

考虑时间节点 \(t_0, t_1, \ldots, t_N=T\) 对应的数值解分别为 \(y_0, y_1, \ldots, y_N\), 那么 Euler 方法可以表述如下 $$ y_{k} = y_{k-1}+ h_k f(x_{k-1}, y_{k-1}), \quad k=1,2, \ldots, N, \tag{Eq. 3} $$ 其中, \(h_k=t_k-t_{k-1}.\)

误差分析

数值方法的提出一般会利用误差估计作为可靠性的保证. 本节我们考虑 Euler 法的误差估计. 对于一个二阶连续可微的光滑函数 \(y(t)\), 泰勒定理表明如下等式 (\(t_k=t_{k-1}+h_{k}\)) $$ y(t_k) = y(t_{k-1})+y'(t_{k-1})h_k+\frac{1}{2}y^{\prime\prime}(s)h_k^2, \quad s\in (t_{k-1}, t_k). \tag{Eq. 4} $$ 如果我们定义误差 \(e_k=y(t_{k})-y_k\), 那么通过比较 (Eq.3)(Eq. 4), 并利用 \(y'(t_{k-1})=f(t_{k-1},y(t_{k-1}))\) 可以得到第 \(k\) 步的局部误差方程为 $$ e_k = e_{k-1}+h_k \left[ f(t_{k-1},y(t_{k-1})) – f(t_{k-1}, y_{k-1}) \right] + \frac{1}{2}y^{\prime\prime}(s)h_k^2 . $$

假设 1.

1. \(f\) 关于 \(y\) 是 Lipschitz 连续的, Lipschitz 常数记为 \(L\).
2. \(y^{\prime\prime}\) 一致有界, 上界记为 \(M\).

为了得到一个简洁的误差估计, 我们令 \(h = \max_{k}h_k\). 则误差有如下估计 $$ |e_k| \leq |e_{k-1}| + h_kL e_{k-1} + \frac{M}{2}h_k^2 = (1+h_kL)|e_{k-1}|+\frac{Mh_k}{2}h. $$ 不等式两边同时加上 \(\frac{Mh}{2L}\), 化简可得 $$ \begin{align*} |e_k| + \frac{Mh}{2L} &\leq (1+h_kL) \left( |e_{k-1}|+\frac{Mh}{2L} \right) && \\ &\leq e^{h_kL} \left( |e_{k-1}|+\frac{Mh}{2L} \right) && (\text{by } 1+x< e^x, x>0)\\ & \leq e^{(h_k+h_{k-1}+\ldots+h_1)L} \left( |e_{0}|+\frac{Mh}{2L} \right) && (\text{by } h_k=t_k-t_{k-1})\\ &= e^{(t_k-t_0)L}\frac{Mh}{2L} .&& (\text{by } e_0=y(t_0)-y_0=0) \end{align*} $$

定理 1. 如果假设 1 成立, 那么当令 \(h = \max_{k}h_k\) 时, Euler 方法的误差满足 $$ |y(t_{k})-y_k| \leq e^{(t_k-t_0)L}\frac{Mh}{2L}. \tag{Eq. 5}$$

舍入误差的影响

由于计算机精度所限, 舍入误差对数值结果的影响不可避免, 直接的结果就是时间步长不能无下限的减小, 因为步长减小到一定程度后, 误差就会从逐渐减小转变到逐渐增大. 具体地, 如果我们在误差方程 (Eq. 3) 中添加舍入误差项 \(\delta_k\) (并且令 \(\delta=\max_{k}\delta_i\)) $$ y_{k} = y_{k-1}+ h_k f(x_{k-1}, y_{k-1})+\delta_i, $$ 利用与前面相同的误差估计方法, 可以得到 $$ |e_k| \leq (1+h_kL)|e_{k-1}|+\frac{Mh_k}{2}h+\delta. $$ 不等式两边同时加上 \(\frac{Mh}{2L}+\frac{\delta}{h_kL}\), 化简可得 $$ |e_k| + \frac{Mh}{2L} +\frac{\delta}{h_kL} \leq e^{(t_k-t_0)L}\left( \frac{Mh}{2L}+ \frac{\delta}{h_kL}\right) . $$ 从而, 我们可以得到与定理 1 类似的一个结果 (一致步长时) $$ |e_k| \leq e^{(t_k-t_0)L} \left( \frac{Mh}{2L}+ \frac{\delta}{hL}\right) \tag{Eq. 6}. $$ 这个结果表明了, 在考虑舍入误差的影响时, Euler 方法的误差估计与时间步长 \(h\) 不再是线性的关系了; 容易验证, 当 \(h^2>\frac{2 \delta}{M}\) 时误差逐渐减小; 否则误差逐渐增大.

数值算例

本节利用数值算例验证 Euler 方法的收敛性和收敛速度, 实现过程见附录中的 Python 代码. 下面, 我们考虑区间 \([0,1]\) 上的单变量常微分方程 (Eq. 1).

例 1. 令 \(f(t,y)=y-t^2+1\), 解析解 \(y=(t+1)^2-0.5e^r.\)

这个数值例子中的右端项满足 Lipschitz 连续并且解析解的二阶导数在区间 [0,1] 上一致有界, 所以, 数值结果中 Euler 方法的收敛速度应该是 \(O(h)\). 我们在下表中列出了在 \(t=1\) 上的误差随步长变化的数值结果, 可以看到, 收敛速度正是如定理 1 给出的一阶收敛. 值得注意的是, 在这个数值例子中, 舍入误差的影响几乎可以忽略, 即使当步长等于 \(1/671088640\) 时, 收敛阶依然符合理论估计.

例 2. 令 \(f(t,y)=-\frac{ty}{1-t^2}\), 解析解 \(y=\sqrt{1-t^2}.\)

例 2 中给出的解析解并不满足误差分析中的条件, 首先 \(f\) 关于 \(y\) 并不是 Lipschitz 连续的, 其次解析解的二阶导数也并非一致有界. 对于这样的一个非光滑解, Euler 方法也是收敛的, 但是并不能达到最优的收敛速度, 见下表. 理论证明笔者并不了解, 希望感兴趣的读者可以查找或研究一下.

附录

# Euler method
import numpy as np
import matplotlib.pyplot as plt

class ode_sol:
    def __init__(self):
        self.y = None
        self.err = None

class my_ode_ivp:
    def __init__(self, fun, t_span, y0, t_eval=None, method='Euler'):
        self.N = 100    # if t_eval is none
        self.fun = fun
        self.t_span = t_span
        self.y0 = y0
        self.sol = ode_sol()

        if t_eval is None:
            self.t_eval = np.linspace(t_span[0], t_span[1], self.N)
        else:
            self.t_eval = t_eval

        self.sol.y = np.zeros(self.t_eval.shape[0])
        self.sol.y[0] = y0

        if method == 'Euler':
            self.euler()
        else:
            pass

    def euler(self):
        for i in range(1, self.t_eval.shape[0]):
            self.sol.y[i] = self.sol.y[i-1] + (self.t_eval[i]-self.t_eval[i-1])*self.fun(self.t_eval[i-1], self.sol.y[i-1])


# setting for ode: f(t,y), y0, t_span, t_eval
# two test examples: ex1, ex2
def f(t, y):
    # return y-t**2+1 #ex1
    return -t*y/(1-t*t) #ex2

# analytic solution
def yfun(t):
    # return (t+1)**2 - .5*np.exp(t) #ex1
    return np.sqrt(1-t*t) #ex2

# y0 = .5 #ex1
y0 = 1. #ex2
t_span = (0, 1.)

N = 60
t_eval = np.linspace(t_span[0], t_span[1], N)

# solve, compute error and plot
ode = my_ode_ivp(fun=f, t_span=t_span, y0=y0, t_eval=t_eval, method='Euler')

y_exact = yfun(ode.t_eval)
ode.sol.err = abs(ode.sol.y[-1]-y_exact[-1])
print(ode.sol.err)

plt.plot(ode.t_eval, y_exact,'.-', label='Exact sol')
plt.plot(ode.t_eval, ode.sol.y, 'o-', label='Euler sol')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()

参考

1. 常微分方程初值问题: 若干数学模型.
2. Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations. In Numerical Methods for Ordinary Differential Equations.
3. Jardine, D. (2011). Euler’s Method in Euler’s Words. In D. Jardine & A. Shell-Gellasch (Eds.), Mathematical Time Capsules: Historical Modules for the Mathematics Classroom (pp. 215-222). Mathematical Association of America.
4. Python 科学计算: NumPy.