计算物理复习[ep.2|数值微分与积分、常&偏微分方程数值解]

本文含有大量数学公式，请使用电脑或平板电脑浏览以获得最佳阅读体验

第三章 数值微分与积分

数值微分

算微分时，我们一般用两点差商法。这种算法无论是理解还是实现起来都非常简单，只需要记住这几个公式即可：

• 一阶向前差商公式

  $$f'(x_i)\approx\frac{f(x_i+h)-f(x_i)}{h}\to\Delta_f f_i=f_{i+1}-f_i \tag{1}$$

• 一阶向后差商公式

  $$f'(x_i)\approx\frac{f(x_i)-f(x_i-h)}{h}\to\Delta_b f_i=f_i-f_{i-1} \tag{2}$$

• 一阶中心差商公式

  $$f'(x_i)\approx\frac{f(x_i+h/2)-f(x_i-h/2)}{h}\to\Delta_c f_i=f_{i+1/2}-f_{i-1/2} \tag{3}$$

• 二阶中心差商公式

  $$f''_i(x_i)\approx\frac{f(x_i+h/2)-2f(x_i)+f(x_i-h/2)}{h^2} \tag{4}$$

可以看到，以上几种差商公式都是以微分的定义为基础构建起来的，它们之间唯一的区别就是差商的方向，即对微分结果产生影响的方式，它们的计算方式也会因此有些不同，我们将会在后文提到。

数值积分

程序求积分原理简介

我们平时所用到的最简单的积分算法是把函数 $f(x)$ 对应的积分区域 $[a,b]$ 等分成 $n$ 份，再以步长乘每一段左端或者右端的函数值来构造矩形，最后把这些矩形的面积加在一起就能得到 $f(x)$ 在 $[a,b]$ 上积分的近似值，如下图：

从矩形的形状来看，这种积分是非常粗糙的，为了提高近似的程度，我们会用梯形来代替矩形进行求和：

显然，在对同一个积分区间分相同段数时，使用梯形的积分结果将更加接近真实值。从这一步优化可以看出，我们对积分精度的追求实质上是对更好的积分模型的追求，下面将依次为基础展开讨论。

牛顿-科茨公式（插值型求积公式）

牛顿-科茨公式通过构造插值函数的方式来使得积分模型更加贴近真实值，具体来说是根据一系列求积节点 $a\leq x_0<x_1<\cdots<x_{n-1}<x_n\leq b$ 和其对应的被积函数值利用拉格朗日插值法来构造插值函数完成积分 $\int_a^b f(x)$。

$$\int_{a}^{b}f(x)d x\approx\int_{a}^{b}L_{n}(x)d x=\int_{a}^{b}\sum_{i}^{n}l_{i}(x)f(x_{i})d x=\sum_{a}^{n}A_{k}f(x_{k}),A_{k}=\int_{a}^{b}l_{k}(x)d x$$

上式表明，随着 $n$ 的改变，插值函数的形式也会改变，其中 $n=1$ 和 $n=2$ 的情况分别对应我们下面要提到的梯形公式和辛普森公式

梯形公式（n = 1 时的牛顿-科茨公式）

先看看 $n=1$ 时牛顿-科茨公式的形式：

$$\int_{a}^{b}f(x)d x\approx\sum_{k=0}^{1}A_{k}f(x_{k})=A_{0}f(x_{0})+A_{1}f(x_{1}) \tag{5}$$

其中:

$$\begin{cases} A_{0}=\int_{a}^{b}l_{0}(x)d x=\int_{a}^{b}{\frac{x-x_{1}}{x_{0}-x_{1}}}d x={\frac{1}{2}}(b-a) \\ A_{1}=\int_{a}^{b}l_{1}(x)d x=\int_{a}^{b}{\frac{x-x_{0}}{x_{1}-x_{0}}}d x={\frac{1}{2}}(b-a) \end{cases} \tag{6}$$

把（6）代入（5）得：

$$\int_{a}^{b}f(x)\mathrm{d}x\approx{\frac{b-a}{2}}[f(a)+f(b)] \tag{7}$$

你会发现，等号的右侧是一个梯形的面积（这也是该公式名称的由来），它直接把整个积分区域当成一整个梯形来算了，说实话这算的也够粗糙的。

辛普森公式（n = 2 时的牛顿-科茨公式）

比起梯形公式，辛普森公式稍显复杂：

$$\int_a^b f(x)dx\approx\sum_{k=0}^2A_k f(x_k)=A_0f(x_0)+A_1f(x_1)+A_2f(x_2) \tag{8}$$

其中：

$$\begin{cases} A_0 =\int_a^b l_0(x)dx=\int_a^b\frac{x-x_1}{x_0-x_1}\cdot\frac{x-x_2}{x_0-x_2}dx=\frac16(b-a) \\ A_1 =\int_a^b l_1(x)dx=\int_{a}^{b}\frac{x-x_{0}}{x_{1}-x_{0}}\cdot\frac{x-x_{2}}{x_{1}-x_{2}}d x=\frac{4}{6}(b-a) \\ A_2=\int_{a}^{b}l_{2}(x)d x=\int_{a}^{b}{\frac{x-x_{0}}{x_{2}-x_{0}}}\cdot{\frac{x-x_{1}}{x_{2}-x_{1}}}d x={\frac{1}{6}}(b-a) \end{cases} \tag{9}$$

把（9）代入（8）得：

$$\int_{a}^{b}f(x)\mathrm{d}x\approx{\frac{b-a}{6}}\bigg[f(a)+f(b)+4f\bigg({\frac{a+b}{2}}\bigg)\bigg] \tag{10}$$

虽然麻烦一点，但是由于 $n = 2$ 时的拉格朗日基函数都是二次函数，所以结果也会更接近真实值。

复化求积公式

无论是梯形公式还是辛普森公式都是把积分区间当成一整块函数进行计算，为了得到更精确的解，我们可以把积分区间等分为若干个子区间，然后在子区间中使用梯形公式或辛普森公式进行求积，就像本节的第一小节中提到的那样。

还记得 1.2.1 里第二个用梯形求积分的图吗？那个实际上就是梯形公式用于复化求积的示意图，使用复化求积往往会获得比单独用牛顿-科茨公式更准确的结果，下面是一段示例代码：

% 复化梯形公式法求积分
f = @(x)(1./(1 + x.^2));
a = 0;
b = 10;
n = 100;
h = (b - a)/n;
xx = linspace(a, b, n + 1);
T = h*sum(f(xx)) - h/2*(f(b) + f(a));
fprintf('%.4f\n',T)

第四章 常微分方程数值解

显式欧拉与隐式欧拉法

对于一阶常微分方程来说，显式欧拉与隐式欧拉的区别就在于使用向前差商或向后差商来表示一阶导数项：

以求解初值问题为例：

$$\begin{cases}y'(x)=f(x,y),\quad a\le x\le b\\ y(0)=y_0\end{cases}$$

• 向前差商 – 显式欧拉法

  $$\begin{array}{c} y'(x_k) = \frac{y(x_{k+1})-y(x_k)}{h} \\ \Downarrow \\ y'(x_k)=f(x_k,y_k)\Rightarrow\frac{y(x_{k+1})-y(x_k)}{h}=f(x_k,y_k) \\ \Downarrow \\ y(x_{k+1})=y(x_k)+hf(x_k,y(x_k)) \end{array} \tag{11}$$

• 向后差商 – 隐式欧拉法

  $$\begin{array}{c} y'(x_k) = \frac{y(x_k)-y(x_{k-1})}{h} \\ \Downarrow \\ y'(x_k)=f(x_k,y_k)\Rightarrow\frac{y(x_k)-y(x_{k-1})}{h}=f(x_k,y_k) \\ \Downarrow \\ y(x_{k})=y(x_{k-1})+hf(x_k,y(x_k)) \end{array} \tag{12}$$

显式欧拉法的递推公式 $y(x_{k+1})=y(x_k)+hf(x_k,y(x_k))$ 依靠前一项的结果，利用初值即可完成推导；隐式欧拉法需要“后一项”的结果，要解方程才能得到，麻烦了点但是好处是每次迭代结果受到前一步的迭代值的影响较小。

龙格-库塔（Runge-Kutta）方法

龙格-库塔方法（以下简称 RK 法）的核心思想是：用 $(x_k,y_k)$ 附近点的函数值的线性组合对 $f(x_k,y_k)$ 修正，近似各级导数项的加和。其一般形式为：

$$y_{k+1}=y_k+h\sum_{i=1}^N c_ik_i \tag{13}$$

四阶 RK 法公式：

$$\begin{aligned} &k_1 =f(x_k,y_k) \\ &k_2 =f(x_k+h/2,y_k+hk_1/2) \\ &k_3 =f(x_k+h/2,y_k+hk_2/2) \\ &k_4 =f(x_k+h,y_k+hk_3) \\ &y_k =y_k+h(k_1+2k_2+2k_3+k_4)/6 \end{aligned} \tag{14}$$

其中 $h$ 是步长。

以求解常微分初值问题

$$\left\{\begin{array}{l} \frac{\mathrm{d}y}{\mathrm{d}x} = xy^2 + 2y \space(0<x<5)\\ y(0) = -5 \end{array}\right.$$

为例，有以下代码：

h = 0.43; % 步长
x = 0:h:5; % 定义域
N = size(x,2);
y = zeros(1,N);
y(1) = -5; % 设置初值
y3 = Runge_Kutta(x,y,N,h);

%Runge-Kutta法
function z = Runge_Kutta(x,y,N,h)
    f = @(x,y)(x*y^2 + 2*y);
    for i = 1:N - 1     %Runge-Kutta
        k1 = f(x(i), y(i));
        k2 = f(x(i) + h/2, y(i) + h*k1/2);
        k3 = f(x(i) + h/2, y(i) + h*k2/2);
        k4 = f(x(i) + h, y(i) + h*k3);
        y(i + 1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4);
    end
    z = y;
end

常微分方程组初值问题

对于常微分方程组，无非是将一系列常微分问题用向量的方式表示，求解过程相似。以下面这个方程组为例：

$$\begin{cases}y_1'=f_1(x,y_1,y_2,\cdots,y_n)&y_1(x_0)=y_{10}\\ y_2'=f_2(x,y_1,y_2,\cdots,y_n)&y_2(x_0)=y_{20}\\.......\\ y_n'=f_n(x,y_1,y_2,\cdots,y_n)&y_n(x_0)=y_{n0}\end{cases} \tag{15}$$

我们引入：

$$\vec{y}=\begin{bmatrix}y_1\\ y_2\\ \vdots\\ y_n\end{bmatrix}\quad \vec{f}(x,\vec{y})= \begin{bmatrix}f_1(x,\vec{y})\\ f_2(x,\vec{y})\\ \vdots\\ f_n(x,\vec{y})\end{bmatrix}\quad \vec{y}_0= \begin{bmatrix}y_1(x_n)\\ y_2(x_n)\\ \vdots\\ y_n(x_n)\end{bmatrix}=\begin{bmatrix}y_{1a}\\ y_{2a}\\ \vdots\\ y_{n0} \end{bmatrix} \tag{16}$$

则（15）可以表示为：

$$\begin{cases}\vec{y}'=\vec{f}(x,\vec{y})\\ \vec{y}(x_0)=\vec{y}_0\end{cases} \tag{17}$$

下面是一个关于野兔和狐狸数量的例题：

在一个生物圈中，只有野兔和狐狸两种动物,记 $y_1$ 为野兔数量, $y_2$ 为狐狸数量。设有足够的食物供野兔生存,而狐狸只以野兔为食物。模型如下：

$$\begin{aligned} \\ \begin{cases} \frac{dy_1}{dx}=y_1-0.01y_1y_2\\ \frac{dy_2}{dx}=-y_2+0.02y_1y_2 \end{cases}& 0 < x < 15\\ \\ y_1(0)=20,y_2(0)=20 \end{aligned}$$

代码部分：

%Runge-Kutta法解微分方程组（野兔和狐狸数量问题）
clear, clc
h = 0.1;
x = 0:h:15;
N = size(x,2);
y = zeros(N,2);
y(1,:) = [20 20];

for i = 1:N - 1     %Runge-Kutta
    k1 = f(x(i), y(i,:));
    k2 = f(x(i) + h/2, y(i,:) + h*k1/2);
    k3 = f(x(i) + h/2, y(i,:) + h*k2/2);
    k4 = f(x(i) + h, y(i,:) + h*k3);
    y(i + 1,:) =y(i,:) + h/6*(k1 + 2*k2 + 2*k3 + k4);
end

hold on
plot(x,y(:,1))
plot(x,y(:,2))
grid on

function z = f(x,y)
    z = zeros(1,2);
    z(1) = y(1) - 0.01*y(1)*y(2);
    z(2) = -y(2) + 0.02*y(1)*y(2);
end

显然，这段代码中的 $y$ 和2.2节里代码中的 $y$ 略显不同，因为2.2节的 $y_i$ 每迭代一次都会产生一个标量，而对于方程组，因为是二元的，所以这里的 $\mathbf{y}_i$ 是一个 $2\times1$ 的向量，其他计算部分都是相同的。

高阶常微分方程初值问题

既然我们已经攻克了一阶常微分方程组的初值问题，那么高阶常微分方程初值问题也将不在话下，因为核心思路就是将求解高阶常微分方程转换为解一阶常微分方程组。

对于一个 n 阶常微分方程：

$$\begin{cases} y^{(n)}=f(x,y,y',\cdots,y^{(n-1)}) \\ y(x_0)=y_0,y'(x_0)=y_0',\cdots,y^{(n-1)}(x_0)=y_0^{(n-1)} \end{cases} \tag{18}$$

引入新向量来表示各阶导数：

$$\boldsymbol{y}=\begin{bmatrix} y_1=y \\ y_2=y' \\ \vdots \\ y_n=y^{(n-1)} \end{bmatrix}\Rightarrow \boldsymbol{y}'=\begin{bmatrix} y_1'=y_2 \\ y_2'=y_3 \\ \vdots \\ y_n'=f(x,y,y',\cdots,y^{(n-1)}) \end{bmatrix} \tag{19}$$

和各阶初值：

$$\boldsymbol{y}(x_0)=\begin{bmatrix} y_{1}{\Big(}x_{0}{\Big)}=y_{0} \\ {y_{2}\left(x_{0}\right)=y_{0}'}\\ {\vdots}\\ {y_{n-1}\left(x_{0}\right)=y_{0}^{(n-2)}}\\ {y_{n}\left(x_{0}\right)=y_{0}^{(n-1)}} \end{bmatrix} \tag{20}$$

有：

$$\begin{cases}\boldsymbol{y}^{'}=\boldsymbol{f}(x,\boldsymbol{y})\\ \boldsymbol{y}(x_0)=\boldsymbol{y}_0\end{cases} \tag{21}$$

接下来的步骤就和解常微分方程组一样了

边值问题打靶法

打靶法的基本思路：将边值问题转化为初值问题考虑。或者说适当选择初始值使初值问题的解满足边值条件。然后用求解初值问题的任一种有效的数值方法求解。

从思路来看，求解边值问题实际上就是一个不断“猜”的过程，我们可以大概估算初值的范围，然后以边界条件为判定依据不断在这个范围内寻找，最后找到一个满足特定精度误差的初值，具体参考下面的例题：

为了使烟花从地面发射 5s 后在距离地面 40m 的空中爆炸，初始的发射速度应该是多大？考虑空气阻力与速度成正比，阻力系数 $\alpha= 0.01$。

列出微分方程：

$$\begin{gathered} \frac{\mathrm{d}^{2}y}{\mathrm{d}t^{2}}=-9.8-0.01\frac{\mathrm{d}y}{\mathrm{d}t} \\ y(0)=0,\quad y(5)=40 \end{gathered}$$

确定初值 $y(0)=0$ 另一个初值 $y'(0)$ 待定，猜测 $y'(0)\in[10,1000]$ ，以 $y(5)=40$ 为目标值，误差要求小于 $\mathrm{eps} = 10^{-5}$，有以下代码：

% 打靶法
clear,clc;

h = 0.05; % 步长
x = 0:h:5;
N = size(x,2);
eps = 1e-5; % 精度
err = 10; % 误差初值
a = 10; % 初速度范围下限
b = 1000; % 初速度范围上限

while (err > eps)
    y = zeros(N,2);
    y(1,:) = [0,(a+b)/2]; % 二分法查找初速度

    for i = 1:N-1
        % R-K
        k1 = f(x(i), y(i,:));
        k2 = f(x(i)+h/2,y(i,:)+h*k1/2);
        k3 = f(x(i)+h/2,y(i,:)+h*k2/2);
        k4 = f(x(i)+h,y(i,:)+h*k3);
        y(i+1,:) = y(i,:)+h/6*(k1+2*k2+2*k3+k4);
    end
    
    err = y(N,1) - 40; % 计算RK法算出来的y(5)和目标值的差距
    if err > 0
        b = (a+b)/2; % 高了，降低初速度
    else
        a = (a+b)/2; % 低了，增加初速度
    end
    err = abs(err);
end

plot(x,y(:,1));
hold on
hold off

function z = f(x,y)
        z = zeros(1,2);
        z(1) = y(2);
        z(2) = -9.8-0.01*y(2);
end

第五章 偏微分方程数值解

常见数值求解PDE的方法

常见数值求解PDE的方法有这~~~~么多种 我都不会！（骄傲）

点击打开列表

• 有限差分法（Finite Difference Method, FDM）：这是一种基于离散化空间和时间的方法，将偏导数用有限差分近似表示。FDM 是求解椭圆、抛物线和双曲线型偏微分方程的常用方法。
• 有限元法（Finite Element Method, FEM）：FEM 是一种将复杂的几何形状划分成简单的子域（通常为三角形或四边形）的方法。在每个子域上，用多项式近似方程的解。FEM 可以求解各种类型的偏微分方程，包括结构力学、流体力学和传热问题。
• 有限体积法（Finite Volume Method, FVM）：FVM 是一种基于守恒定律的方法，将计算域划分为有限数量的控制体积。通过求解在每个控制体积内的流量平衡，可以近似求解偏微分方程。FVM 常用于流体力学和传热问题。
• 谱方法（Spectral Method）：这种方法使用正交多项式（如切比雪夫多项式、Legendre 多项式等）作为基函数来表示方程的解。谱方法在解空间和频率域具有优秀的性能，特别适用于周期性和光滑问题。
• 边界元法（Boundary Element Method, BEM）：BEM 是一种基于边界积分方程的方法，只需要在边界上离散点而无需在整个计算域内进行离散。BEM 特别适用于无穷大区域、多层介质和多物体问题。
• 网格自适应方法（Meshless Method）：这些方法不依赖于传统的网格划分，而是使用无网格的点集来表示计算域。常见的无网格方法有光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）、移动最小二乘法（Moving Least Squares, MLS）等。这些方法适用于具有复杂边界和变形的问题。

本章主要讲解其中的有限差分法，解题的三个步骤为：

1. 选取差分格式，将微分方程离散化为差分方程
2. 确定网格步长，确保差分方程的解收敛于微分方程的解
3. 求解相应的迭代方程或代数方程组

有限差分法蛙跳格式

蛙跳格式的时间、空间均取中心差商，以对流方程为例：

$$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0 \tag{22}$$

将时间、空间一阶中心差商代入，整理得到迭代公式：

$$\begin{array}{c} \frac{u_i^{k+1}-u_i^{k-1}}{2\Delta t}+\frac{a}{2\Delta x}\left(u_{i+1}^{k}-u_{i-1}^{k}\right)=0\\ \Downarrow\\ u_i^{k+1}=u_i^{k-1}+\frac{a\Delta t}{\Delta x}\left(u_{i+1}^{k}-u_{i-1}^{k}\right) \end{array} \tag{23}$$

其中 $k$ 是有关时间的下标， $i$ 是和空间 $x$ 有关的下标，递推过程如图所示：

FTCS、BTCS 格式

FTCS 是 Forward Time Central Space（时间前向差商，空间中心差商）的简写，同理可以衍生出 BTCS, FTFS, FTBC, CTCS etc. 好家伙，搁这排列组合呢

话不多说，直接用例题解释：

设杆一端 $x=0$ 的温度始终为0度，另一端 $x=1m$ 的温度始终为100度，初始时杆整体为0度，$\kappa = 0.835$。采用 FTCS 格式数值计算 $0.1s$ 内的一维热传导方程：

$$\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial x^2} \tag{24}$$

FTCS 求解过程：

首先，用差商替换等号两端的导数项，写出 $(i,k)$ 处的差分方程：

$$\frac{T_i^{k+1}-T_i^k}{\Delta t}=\kappa\frac{T_{i+1}^k-2T_i^k+T_{i-1}^k}{(\Delta x)^2} \tag{25}$$

得到递推公式：

$$T_i^{k+1}=T_i^k+\frac{\kappa\Delta t}{(\Delta x)^2}(T_{i+1}^k-2T_i^k+T_{i-1}^k) \tag{26}$$

然后设置初值和精度 $\Delta t=0.0001,\Delta x=0.02,\mathrm{eps} = 10^{-5}$ 在程序里用一个二重循环进行迭代，每次迭后进行判敛，以下是实现代码：

clear, clc
%有限差分法FTCS格式解偏微分方程u(t,x)

dt = 0.0001;
dx = 0.02;

t = 0:dt:0.1;
x = 0:dx:1;
K = size(t,2);
I = size(x,2);

u = zeros(I,K);
u(I,:) = 100;

eps = 1e-5;
err = 10;
r = 0.835*dt/dx^2;
while (err > eps)
    u1 = u;
    for i = 2:I - 1
        for k = 1:K - 1
            u(i,k + 1) = u(i,k) + r*(u(i + 1,k) - 2*u(i,k) + u(i - 1,k));
        end
    end
    err = norm(u - u1,1);
end

heatmap(u,'GridVisible','off');

FTCS 格式虽然容易理解和实现，但是该方法在步长过大的情况下很容易发散，所以我们再考虑 BTCS 格式。

BTCS 求解过程：

同样首先写出 $(i,k)$ 处的差分方程，但是时间改为向后差商：

$$\frac{T_i^k-T_i^{k-1}}{\Delta t}=\kappa\frac{T_{i+1}^k-2T_i^k+T_{i-1}^k}{(\Delta x)^2} \tag{27}$$

令 $r=\frac{\kappa\Delta t}{(\Delta x)^2}$ 得到递推公式：

$$T_i^{k-1}=(1+2r)T_i^k-rT_{i+1}^k-rT_{i-1}^k \tag{28}$$

由于向后差商是隐式的，所以我们需要解方程来得到下一步迭代的值，由（28）得到的方程组我们可以用下面这个矩阵表示，其中 $T_1^k=0$ 和 $T_N^k=100$ 是边界条件：

$$\begin{bmatrix} 1 + 2r&-r\\ -r&1 + 2r&-r\\ &\ddots&\ddots&\ddots\\ &&-r&1 + 2r&-r \\&&& -r & 1 + 2r \end{bmatrix} \begin{bmatrix} T_2^k\\ T_3^k\\ \vdots\\ T_{N-2}^k \\T_{N-1}^k \end{bmatrix}= \begin{bmatrix} T_2^{k-1}+ rT_1^k\\ T_3^{k-1}\\ \vdots\\ T_{N-2}^{k-1}\\ T_{N-1}^{k-1}+ rT_N^k \end{bmatrix} \tag{29}$$

显然，系数矩阵是一个三对角矩阵，我们可以用追赶法快速求解，实现代码如下：

%有限差分法BTCS格式解偏微分方程u(t,x)
clear, clc
dt = 0.001;
dx = 0.02;
x = 0:dx:1;
t = 0:dt:0.1;
I = size(x, 2);
K = size(t, 2);

u = zeros(I, K);
u(I, :) = 100;
kappa = 0.835;
r = kappa*dt/dx^2;

a = -r*ones(1, I - 2);
a(1) = 0; % 初始化a向量
b = (1 + 2*r)*ones(1, I - 2); % 初始化b向量
c = -r*ones(1, I - 2);
c(I - 2) = 0; % 初始化c向量
eps = 1e-5;
err = 10;

while (err > eps)
    u1 = u;
    for k = 2:1:K
        d = zeros(1, I - 2);
        d(1) = u(2, k - 1) + r*u(1, k);
        d(I - 2) = u(I - 1, k - 1) + r*u(I, k);
        d(2:I - 3) = u(3:I - 2, k - 1)'; % 计算d向量
        T = chase(a, b, c, d, I - 2);
        u(2:I - 1, k) = T';
    end
    err = norm(u - u1, 1);
end
heatmap(u,'GridVisible','off');

function x = chase(a,b,c,d,n) % 追赶法函数
    w = zeros(1, n);
    g = zeros(1, n);
    
    w(1) = c(1)/b(1);
    g(1) = d(1)/b(1);
    
    for i=2:n
        w(i) = c(i)/(b(i)-a(i)*w(i-1));
        g(i) = (d(i)-a(i)*g(i-1))/(b(i)-a(i)*w(i-1));
    end
    
    x = zeros(1,n);
    x(n) = g(n);
    for i = n-1:-1:1
        x(i) = g(i)-w(i)*x(i+1);
    end
end

Crank-Nicolson 格式

Crank-Nicolson 格式（以下简称 CN 法）结合了 FTCS 和 BTCS 的计算方式，仍然用上一个小结中的例子来解释：

我们结合 FTCS 和 BTCS 的差分方程：

$$\begin{cases} \frac{T_i^{k+1}-T_i^k}{\Delta t}=\kappa\frac{T_{i+1}^k-2T_i^k+T_{i-1}^k}{(\Delta x)^2} \\ \frac{T_i^k-T_i^{k-1}}{\Delta t}=\kappa\frac{T_{i+1}^k-2T_i^k+T_{i-1}^k}{(\Delta x)^2} \end{cases}$$

得到 CN 法的差分方程：

$$\frac{T_i^{k+1}-T_i^k}{\Delta t}=\frac{\kappa}{2(\Delta x)^2}\left(T_{i+1}^{k+1}-2T_i^{k+1}+T_{i-1}^{k+1}+T_{i+1}^k-2T_i^k+T_{i-1}^k\right) \tag{30}$$

其中上标为 $k+1$ 的都是未知信息，上标为 $k$ 的是已知信息。由于有未知信息（隐式），显然我们又需要通过解方程组来获得这些未知信息了。

整理（30）得：

$$-BT^{k+1}_{i+1}+CT^{k+1}_i-BT^{k+1}_{i-1}=BT^{k}_{i+1}+AT^{k}_i+BT^{k}_{i-1} \tag{31}$$

其中：

$$\begin{cases} A=\frac{1}{\Delta t}+\frac{\kappa}{(\Delta x)^2}\\ B = \frac{\kappa}{(2\Delta x)^2} \\ C=\frac{1}{\Delta t}-\frac{\kappa}{(\Delta x)^2} \end{cases} \tag{32}$$

得到矩阵：

$$\begin{bmatrix} A&-B\\ -B&A&-B\\ &\ddots&\ddots&\ddots\\ &&-B&A&-B\\ &&& -B&A \end{bmatrix} \begin{bmatrix} T_2^{k+1}\\ T_3^{k+1}\\ \vdots\\ T_{N-2}^{k+1} \\ T_{N-1}^{k+1} \end{bmatrix}= \begin{bmatrix} BT_1^{k+1}\\ BT_3^{k}+CT_2^{k}+BT_1^{k}\\ \vdots\\ BT_N^{k}+CT_{N-1}^{k}+BT_{N-2}^{k}\\ BT_N^{k+1} \end{bmatrix} \tag{33}$$

最后追赶法求解即可。

終了

写完辣！这篇博客我已经尽我所能精简内容了，比起上一篇两个章节八千来字，这一篇三个章节才五千七百多字，嗯嗯。。。很不错了。。。

我发现用例题和代码来讲解是一个不错的偷懒手段，对于理解一个问题或者算法，让我看代码比单纯看数学公式要清楚许多，所以这篇博客有很多长长的代码块

感谢你能看完这篇文章，有误之处欢迎反馈！