数值ODE（二）—— 欧拉差分

$\textbf{初值问题：}\frac{\mathrm{d}y}{\mathrm{d}x}=f(x, y), \quad y(x_0)=y_0$

离散形式

初值问题的数值解法，其实就是就是寻求 $y(x)$ 在一系列离散点 ${a}={x}_{0}<x_{1}<\cdots<{x}_{N}={b}$ 上的近似值，相邻两点的间距 $h_n=x_{n+1}-x_n$ ，一般步长都取相等的值，也就是有 $h_{n} \equiv h \quad(n=0,1,2, \cdots)$ ，于是 $x_n=x_0+nh$ 。我们用数值求解的函数 $y$ 的近似，就是若干段分段折线。

所以接下来要面临的问题就是，如何离散化原微分方程？

一种最显然的方式，把常数 $f(x_0, y_0)$ 当作 $y$ 在 $(x_0, x_1)$ 上的导函数近似值，这样就有 $y_1=y_0+hf(x_0, y_0)$ ，就这样依次类推下去，计算公式就是： ${y}_{n+1}={y}_{n}+{h} f\left({x}_{n}, {y}_{n}\right), \quad {n}={0}, {1}, {2}, \cdots$

这种方法叫做显式欧拉(Explicit Euler)。是数值ODE里可以说是最简单的方法吧。
直觉告诉我们，当离散间隔 $h$ 取的足够小的时候，数值解是能对原函数有好的近似的。

执行如下Julia代码：

# using Plots; using LaTeXStrings
p = []
for h in [0.25, 0.2, 0.1, 0.05]
    x = 0 : h : 5
    y = [1.0]
    for i = 1 : length(x) - 1
        push!(y, y[i] + h * y[i])
    end
    plot(x, y, label="h=$h")
    plot!(exp, x, label=L"y=e^x")
    t = plot!(x, exp.(x) - y, label="error")
    push!(p, t)
end
plot(p[1],p[2],p[3],p[4], layout=(2, 2), legend=:topleft)

使用显式欧拉方法对 $\frac{dy}{dx}=y, \;\;y|_{x=0}=1$ 用不同的步长做数值近似，其结果如下图：

图1：显式欧拉

显式欧拉公式是用向前差商来代替导数的： $y^{\prime}\left(x_{n}\right) \approx \frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}\approx f(x_n, y_n)$

如果使用向后差商代替导数： $y^{\prime}\left(x_{n+1}\right) \approx \frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$ 就得到了隐式欧拉公式(Implicit Euler)：
$\left\{\begin{array}{l} y_{n+1}=y_{n}+h f\left(x_{n+1}, y_{n+1}\right) \\ y_{0}=y\left(x_{0}\right) \end{array} \quad n=0,1,2, \cdots\right.$

显式欧拉是可以一步一步迭代进行计算的。隐式欧拉，需要注意的是，已知 $(x_n, y_n)$ 时求 $y_{n+1}$ 时是通过方程 $y_{n+1}=y_{n}+h f\left(x_{n+1}, y_{n+1}\right)$ 来完成的！所以隐式欧拉要比显式欧拉求解起来要困难一些。

如果 $f$ 是比较好的函数，那么可以直接写出 $y_{n+1}$ 与 $y_n, x_{n+1}$ 之间的关系。否则可能就要解一个非线性方程，幸运的是，在 $hL<1$ 的时候，（ $L$ 是 $f(x,y)$ 的Lipschitz常数），可以使用迭代法来求解这个方程：

$\left\{\begin{array}{l} {y}_{n+1}^{(0)}={y}_{n}+{h} {f}\left({x}_{n}, {y}_{n}\right) \\ {y}_{n+1}^{(k+1)}={y}_{n}+{h} f\left({x}_{n+1}, {y}_{n+1}^{(k)}\right), \quad {k}={0}, {1}, {2}, \cdots \end{array}\right.$

隐式欧拉和显式欧拉是最基础的两种离散形式。在求解 $y_{n+1}$ 的时候只用到了 $y_n$ 的信息，这样的方法叫做单步法。

如果使用中心差商代替导数： $y^{\prime}\left(x_{n}\right) \approx \frac{y\left(x_{n+1}\right)-y\left(x_{n-1}\right)}{2 h}$ 就得到了两点欧拉公式 $y_{n+1}=y_{n-1}+2 h f\left(x_{n}, y_{n}\right), \quad n=0,1,2, \cdots$ 这个公式在计算 $y_{n+1}$ 的时候要用到 $y_n$ 和 $y_{n-1}$ ，也就是前两步的信息，这样的方法叫做多步法。在应用两点欧拉公式的时候，除了给出的初值条件里的 $y_0$ ，还需要用显式欧拉或者隐式欧拉方法求出 $y_1$ ，然后再进行迭代计算。因为两点欧拉公式应用到了两步的信息，所以通常比单步的方法要精确。

以下作出类似上面显式欧拉公式的图：

图2：隐式欧拉

图3：两点法

$\frac{dy}{dx}=ay, \; (a\neq 0)$ 常用来作为测试函数，用来测试各种数值格式的性能。

误差分析

给定一个具体的方法，我们还要去研究这个数值方法产生的误差有多大。关于要研究什么样的误差，看这里。

两大类误差，一种是离散形式来近似的截断误差，另一种是迭代格式中不断积累的计算误差。

局部截断误差是在假设 $y_n=y(x_n)$ 被精确计算的时候， $y(x_{n+1})$ 和 $y_{n+1}$ 的差。
对于显式欧拉方法： $\begin{aligned} y\left(x_{n+1}\right)=y\left(x_{n}+h\right) &=y\left(x_{n}\right)+h y^{\prime}\left(x_{n}\right)+\frac{h^{2}}{2 !} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right) \\ &=y\left(x_{n}\right)+h f\left(x_{n}, y_{n}\right)+\frac{h^{2}}{2 !} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right) \end{aligned}$ 于是： $y\left(x_{n+1}\right)-y_{n+1}=\frac{h^{2}}{2 !} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right)=O\left(h^{2}\right)$ 同理，对于隐式欧拉方法： $y\left(x_{n+1}\right)-y_{n+1}=-\frac{h^{2}}{2 !} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right)=O\left(h^{2}\right)$

整体截断误差是在不考虑舍入误差的情况下, $y(x_{n+1})$ 与 $y_{n+1}$ 之间的误差。记 $e_{n+1}=y\left(x_{n+1}\right)-y_{n+1}$ ，可以证明，在有限的闭区间内，显式欧拉公式和隐式欧拉公式的整体截断误差与步长 $h$ 呈线性关系： $|e_{n+1}|\leq Ch$ 就是说，在 $[0,T]$ 这个区间上，如果我们把步长缩小为原来的一半，那么得到的数值解在 $x=T$ 这个位置的截断误差也缩小为原来的一半，这条性质可以在图1和图2中得到验证。

so，那么问题来了，是不是把步长取的越小越好呢？
现在我们把之前数值模拟的那个ODE做一个小小的改动，考虑 $\frac{dy}{dx}=2y, \;\;y|_{x=0}=1$ 这个方程。