对一些有限差分方法的数学推导

简介

本文默认读者对有限差分方法有一定了解，因此不做过多的前置知识介绍。考虑如下初值问题:

{\begin{matrix} 𝑢^{'} = 𝑓 (𝑡, 𝑢), 𝑡 \in [0, 𝑇] \\ 𝑢 (0) = 𝑢_{0} \end{matrix}

将 $[0, 𝑇]$ 划分为 $𝑀$ 个等距区间， $ℎ = 𝑇 / 𝑀$ ，考虑如下几种差分方法：

向前 Euler 格式 (FE): $𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .$
向后 Euler 格式 (BE): $𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ 𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)}) .$
跃点格式 (LF): $𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛 - 1)} + 2 ℎ 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .$
梯形格式 (Tz): $𝑢^{(𝑛 + 1)} = 𝑢^{(𝑛)} + ℎ \frac{𝑓 (𝑓^{(𝑛)}, 𝑢^{(𝑛)}) + 𝑓 (𝑡^{(𝑛 + 1)}, 𝑢^{(𝑛 + 1)})}{2} .$

接下来，对以上四种有限差分格式，我们将进行其准确性 (Accuracy) 与稳定性 (Stability) 的数学推导。

准确性推导

首先定义局部截断误差 (Local Truncation Error) $𝑅^{(𝑛)}$ ，其定义¹¹ 此处的局部截断误差与通常定义的局部截断误差相差一个 $ℎ$ ，参考时请注意。为：其他点处的精确解在 $𝑡^{(𝑛)}$ 处按照差分格式计算得到的结果与 $𝑡^{(𝑛)}$ 处的精确解的差值。

向前 Euler 格式

向前 Euler 格式的截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛 + 1)}) - 𝑢 (𝑡^{(𝑛)})}{ℎ} - 𝑢^{'} (𝑡^{(𝑛)}) .

对 $𝑢 (𝑡^{(𝑛 + 1)})$ 进行 Taylor 展开，我们有

𝑢 (𝑡^{(𝑛 + 1)}) = 𝑢 (𝑡^{(𝑛)}) + ℎ 𝑢^{'} (𝑡^{(𝑛)}) + \frac{ℎ^{2}}{2} 𝑢^{″} (𝜉^{(𝑛)}), 𝜉^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

于是有

𝑅^{(𝑛)} = \frac{ℎ}{2} 𝑢^{″} (𝜉^{(𝑛)}) = 𝑂 (ℎ), 𝜉^{(𝑛)} \in [𝑡^{(𝑛)}, 𝑡^{(𝑛 + 1)}] .

误差 $𝐸^{(𝑛)}$ 可表示为如下形式：

\frac{𝐸^{(𝑛 + 1)} - 𝐸^{(𝑛)}}{ℎ} + 𝑅^{(𝑛)} = 𝑓 (𝑡^{(𝑛)}, 𝑢 (𝑡^{(𝑛)})) - 𝑓 (𝑡^{(𝑛)}, 𝑢^{(𝑛)}) .

令 $𝑅 = \max_{𝑛} | 𝑅^{(𝑛)} |$ ，对误差 $𝑒^{(𝑛)}$ 有如下估计

\begin{matrix} | 𝐸^{(𝑛 + 1)} | & \leq | 𝐸^{(𝑛)} | + ℎ 𝑅 + ℎ 𝐿 | 𝑢 (𝑡^{(𝑛)}) - 𝑢^{(𝑛)} | = (1 + ℎ 𝐿) | 𝐸^{(𝑛)} | + ℎ 𝑅 \\ \leq {(1 + ℎ 𝐿)}^{𝑛 + 1} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} ({(1 + ℎ 𝐿)}^{𝑛 + 1} - 1) \\ \leq 𝑒^{𝐿 𝑇} | 𝐸^{(0)} | + \frac{𝑅}{𝐿} (𝑒^{𝐿 𝑇} - 1) \leq 𝑂 (ℎ), 𝑛 = 0, 1, \dots, 𝑀 - 1 . \end{matrix}

向后 Euler 格式

向后 Euler 格式的截断误差 $𝑅^{(𝑛)}$ 可表示为如下等式

𝑅^{(𝑛)} = \frac{𝑢 (𝑡^{(𝑛)}) - 𝑢 (𝑡^{(𝑛 - 1)})}{ℎ} - 𝑢^{'} (𝑡^{(𝑛)}) .

对 $𝑢 (𝑡^{(𝑛 - 1)})$ 进行 Taylor 展开，我们有

TODO: 完成全部迁移工作

数值试验

考虑如下初值问题

{\begin{matrix} 𝑢^{'} (𝑡) = 𝑓 (𝑡, 𝑢 (𝑡)) = - \cos (5 𝜋 𝑡) \tan (5 𝜋 𝑢), \\ 𝑢 (0) = 1 / 60 . \end{matrix}

其中 $𝑡 \in [0, 1]$ 。上述初值问题的解析解为

𝑢 (𝑡) = \frac{1}{5 𝜋} \arcsin (\frac{\sqrt{6} - \sqrt{2}}{4 \exp (\sin (5 𝜋 𝑡))}) .

代码实现

向前欧拉格式：

1
function [x, result] = forward_euler(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        result(k) = result(k - 1) + h * f(x(k - 1), result(k - 1));
8
    end
9
end

向后欧拉格式：

1
function [x, result] = backward_euler(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        temp_func = @(uk) result(k - 1) + h * f(x(k), uk) - uk;
8
        result(k) = fzero(temp_func, result(k - 1));
9
    end
10
end

跃点格式：

1
function [x, result] = leap_frog(T, step, u0, u1, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    result(2) = u1;
7
    for k = 3:step + 1
8
        result(k) = result(k - 2) + 2 * h * f(x(k - 1), result(k - 1));
9
    end
10
end

梯形格式：

1
function [x, result] = trapezoidal(T, step, u0, f)
2
    x = linspace(0, T, step + 1);
3
    h = T / step;
4
    result = zeros(1, step + 1);
5
    result(1) = u0;
6
    for k = 2:step + 1
7
        temp_func = @(uk) result(k - 1) - uk + ...
8
            h * (f(x(k), uk) + f(x(k - 1), result(k - 1))) / 2;
9
        result(k) = fzero(temp_func, result(k - 1));
10
    end
11
end