數(shù)值ODE（二）—— 歐拉差分

$\textbf{初值問題：}\frac{\mathrmdph1pf7y}{\mathrm9pbx3blx}=f(x, y), \quad y(x_0)=y_0$

離散形式

初值問題的數(shù)值解法兑牡，其實就是就是尋求 $y(x)$ 在一系列離散點 ${a}={x}_{0}<x_{1}<\cdots<{x}_{N}=$ 上的近似值税灌，相鄰兩點的間距 $h_n=x_{n+1}-x_n$ 均函，一般步長都取相等的值亿虽，也就是有 $h_{n} \equiv h \quad(n=0,1,2, \cdots)$ ，于是 $x_n=x_0+nh$ 苞也。我們用數(shù)值求解的函數(shù) $y$ 的近似洛勉，就是若干段分段折線。

所以接下來要面臨的問題就是如迟，如何離散化原微分方程坯认？

一種最顯然的方式，把常數(shù) $f(x_0, y_0)$ 當(dāng)作 $y$ 在 $(x_0, x_1)$ 上的導(dǎo)函數(shù)近似值氓涣，這樣就有 $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)劳吠。是數(shù)值ODE里可以說是最簡單的方法吧。
直覺告訴我們巩趁，當(dāng)離散間隔 $h$ 取的足夠小的時候痒玩，數(shù)值解是能對原函數(shù)有好的近似的。

執(zhí)行如下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$ 用不同的步長做數(shù)值近似议慰，其結(jié)果如下圖：

圖1：顯式歐拉

顯式歐拉公式是用向前差商來代替導(dǎo)數(shù)的： $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)$

如果使用向后差商代替導(dǎo)數(shù)： $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.$

顯式歐拉是可以一步一步迭代進(jìn)行計算的蠢古。隱式歐拉，需要注意的是别凹，已知 $(x_n, y_n)$ 時求 $y_{n+1}$ 時是通過方程 $y_{n+1}=y_{n}+h f\left(x_{n+1}, y_{n+1}\right)$ 來完成的草讶！所以隱式歐拉要比顯式歐拉求解起來要困難一些。

如果 $f$ 是比較好的函數(shù)炉菲，那么可以直接寫出 $y_{n+1}$ 與 $y_n, x_{n+1}$ 之間的關(guān)系堕战。否則可能就要解一個非線性方程，幸運的是拍霜，在 $hL<1$ 的時候嘱丢，（ $L$ 是 $f(x,y)$ 的Lipschitz常數(shù)），可以使用迭代法來求解這個方程：

$\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.$

隱式歐拉和顯式歐拉是最基礎(chǔ)的兩種離散形式祠饺。在求解 $y_{n+1}$ 的時候只用到了 $y_n$ 的信息越驻，這樣的方法叫做單步法。

如果使用中心差商代替導(dǎo)數(shù)： $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īng)用兩點歐拉公式的時候勺鸦，除了給出的初值條件里的 $y_0$ 诵棵，還需要用顯式歐拉或者隱式歐拉方法求出 $y_1$ ，然后再進(jìn)行迭代計算祝旷。因為兩點歐拉公式應(yīng)用到了兩步的信息履澳，所以通常比單步的方法要精確嘶窄。

以下作出類似上面顯式歐拉公式的圖：

圖2：隱式歐拉

圖3：兩點法

$\frac{dy}{dx}=ay, \; (a\neq 0)$ 常用來作為測試函數(shù)，用來測試各種數(shù)值格式的性能距贷。

誤差分析

給定一個具體的方法柄冲，我們還要去研究這個數(shù)值方法產(chǎn)生的誤差有多大。關(guān)于要研究什么樣的誤差忠蝗，看這里现横。

兩大類誤差，一種是離散形式來近似的截斷誤差阁最，另一種是迭代格式中不斷積累的計算誤差戒祠。

局部截斷誤差是在假設(shè) $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}$ 配阵，可以證明馏颂，在有限的閉區(qū)間內(nèi)，顯式歐拉公式和隱式歐拉公式的整體截斷誤差與步長 $h$ 呈線性關(guān)系： $|e_{n+1}|\leq Ch$ 就是說棋傍，在 $[0,T]$ 這個區(qū)間上救拉，如果我們把步長縮小為原來的一半，那么得到的數(shù)值解在 $x=T$ 這個位置的截斷誤差也縮小為原來的一半瘫拣，這條性質(zhì)可以在圖1和圖2中得到驗證亿絮。

so，那么問題來了麸拄，是不是把步長取的越小越好呢壹无？
現(xiàn)在我們把之前數(shù)值模擬的那個ODE做一個小小的改動，考慮 $\frac{dy}{dx}=2y, \;\;y|_{x=0}=1$ 這個方程感帅。