计算物理（一）基础运算与微分方程求解

本篇文章是西安交大《计算物理与程序设计》课程的笔记，程序语言为 MATLAB。内容主要包括 MATLAB 基础语法，数值积分方法，方程求根，常微分方程的初值问题，常微分方程的边值问题和偏微分方程的边值问题。

MATLAB基础语法

标量向量矩阵

1. 复数 $\mathrm{i}$ 在MATLAB中应当写作 1i 。
2. 注意向量是行向量还是列向量。
3. 矩阵的非共轭转置为 transpose ，共轭转置为 ctranspose 。

函数

匿名函数：f = @(x) expr(x) ，这样创建的 f 是函数句柄 (function handle) 。

普通函数：function y = f(x) ，这样创建的 f 是函数本身，调用函数句柄需要使用 @f 。

数值积分

内置命令

MATLAB内置的一维数值积分命令为 integral ，语法为

I = integral(f,xmin,xmax)

integral 使用全局自适应积分和默认误差容限在 x_\min 至 x_\max 间以数值形式为函数 $f$ 求积分，积分限为 $+\infty$ 时可用 Inf 表示，积分限为 $-\infty$ 时可用 -Inf 表示。

类似的二维数值积分和三维数值积分命令为 integral2 和 integral3 ，语法为

integral2(fun,xmin,xmax,ymin,ymax)
integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)

梯形公式

以下讨论中，总是将积分区间 $(a,b)$ 切割为 $N$ 份，并设初始节点 $x_1=a$，末端节点 $x_{N+1}=b$，积分步长 $h$ 为

$$h = \frac{b-a}{N}$$

梯形公式来源于线性插值。在区间 $[x_i,x_{i+1}]$ 上梯形公式给出的数值积分为

$$I_{\text{trapz}} = \frac{h}{2} \left[ f(x_i)+f(x_{i+1})\right]$$

函数 $f(x)$ 在区间 $(a,b)$ 上的梯形数值积分公式为

$$I_{\text{trapz}} = \frac{h}{2} f(a)+ h \sum_{i=2}^N f(x_i) +\frac{h}{2} f(b)$$

MATLAB代码为

function numint = numint_trapz(f,x_min,x_max,n)
% numint_trapz.m
h = (x_max-x_min)/n;
x = x_min:h:x_max;
fx = f(x);
numint = h/2*(f(x_min)+f(x_max)+2*sum(fx(2:n)));
end

Simpson公式

Simpson公式来源于二次插值，涉及连续三个节点，三个节点夹着两个区间，因此 $N$ 必须是 $2$ 的倍数。

在区间 $[x_i,x_{i+2}]$ 上Simpson公式给出的数值积分为

$$I_{\text{simpon}} = \frac{h}{3} \left[ f(x_i)+ 4f(x_{i+1}) + f(x_{i+2} )\right]$$

函数 $f(x)$ 在区间 $(a,b)$ 上的Simpson数值积分公式为

$$I_{\text{simpson}} = \frac{h}{3} f(a)+ \frac{4h}{3} \sum_{\scriptstyle i=2 \atop \scriptstyle \text{sep}=2}^N f(x_i) + \frac{2h}{3} \sum_{\scriptstyle i=3 \atop \scriptstyle \text{sep}=2}^{N-1} f(x_i) + \frac{h}{2} f(b)$$

MATLAB代码为

function numint = simpson(f,x_min,x_max,n)
% simpson.m
% f is a known function, n must be odd
h = (x_max-x_min)/n;
x = x_min:h:x_max;
fx = f(x);
numint = h/3*(f(x_min)+f(x_max)+4*sum(fx(2:2:n))+2*sum(fx(3:2:n-1)));
end

方程的数值解

内置命令

x = fzero(f,x0) 从 $x_0$ 开始，尝试求出 $f(x) = 0$ 的解并赋给 $x$ 。

命令要求根必须是 $f(x)$ 变号的位置，因此 fzero 无法求解函数零点同时为函数极值点的问题，例如 $x^2=0$ 。

x = fsolve(F,x0) 从 $\bm{x}_0$ 开始尝试求解向量方程/ $n$ 元非线性方程组 $\boldsymbol{F}(\boldsymbol{x}) = \boldsymbol{0}$ 的解 $\bm{x} = [x_1, x_2, \cdots, x_n]^T$ 。

例如求解非线性方程组对应的命令为

$$\begin{cases} e^{-\left(x_1+x_2\right)}-x_2\left(1+x_1^2\right)=0 \\ x_1 \cos \left(x_2\right)+x_2 \sin \left(x_1\right)-\frac{1}{2}=0 \end{cases}$$

x = fsolve(F,[0,0]);
F(1) = exp(-(x(1)+x(2))) - x(2)*(1+x(1)^2);
F(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;

二分法

二分法利用二分查找逐步缩小根 $x_0$ 所在的区间 $[a,b]$ ，直到 $| a-b | <\varepsilon$ 。

function [root,deviation] = bisection(f,a,b,epsilon)
if f(a)*f(b)>0
    warning('No root in the interval.');
else
    while abs(b-a) > epsilon
        mid = (a+b)/2;
        if f(mid)*f(a)<0
            b = mid;
        else
            a = mid;
        end
    end
end
root = (a+b)/2;
deviation = abs(f(a)-f(b));
end

牛顿法

牛顿法需要知晓 $f(x)$ 的导数 $f^\prime(x)$ 的解析表达式，并给定一个启动值 $x_0$ 。Newton - Raphson 迭代公式为

$$x_{i+1} = x_i - \frac{f(x_i)}{f^\prime (x_{i+1})}$$

当 $| x_{i+1} - x_i | <\varepsilon$ 且 $f(x_i)<\delta$ 时视为找到了根。

MATLAB代码为

function root = newton(f,f_derivative,x0,delta,epsilon)
% newton.m
x1 = x0 - f(x0)/f_derivative(x0);
while abs(x1-x0)>epsilon || abs(f(x1))>delta
    x2 = x1 - f(x1)/f_derivative(x1);
    x0 = x1;
    x1 = x2;
end
root = x1;
end

弦割法

弦割法迭代公式需要两个相近的启动值 $x_0,x_1$，并用割线的斜率近似代替导数的斜率。迭代公式为

$$x_{i+1} = x_i - f(x_i) \frac{x_i - x_{i-1}}{f(x_i) - f(x_{i-1})}$$

function root = secant(f,x0,x1,epsilon,delta)
% secant.m
while abs(x1-x0)>epsilon || abs(f(x1))>delta
    x2 = x1 - (x1-x0)/(f(x1)-f(x0))*f(x1);
    x0 = x1;
    x1 = x2;
end
root = x1;
end

预搜索弦割

当恰当的启动值未知时，可从某点开始进行预搜索，找到有根区间后再应用上述三种方法求数值解，程序以弦割法为例。

简单搜索结合弦割法求根：$x_0$ 为启动值，step 为搜索步长，搜索到异号区间 $[x_i,x_i + \text{step}]$ 时应用弦割法求解区间内的根。

function root = presecant(f,x0,step,epsilon,delta)
% presecant.m
while f(x0+step)*f(x0)>0
    x0 = x0 + step;
end
root = secant(f,x0,x0+step,epsilon,delta);
end

常微分方程的初值问题

内置命令

一阶非线性常微分方程的普遍形式为

$$\left\{ \begin{matrix} \dfrac{\mathrm{d}y}{\mathrm{d}x} = f(x,y) \\ y(x_0) = y_0 \end{matrix} \right. \implies y=y(x)$$

找出一阶非线性常微分方程的解法后，可以很容易将其推广至常微分方程组，只需将待求解函数改为向量函数即可

$$\left\{ \begin{matrix} \dfrac{\mathrm{d}\boldsymbol{y}}{\mathrm{d}x} = f(x,\boldsymbol{y}) \\ \boldsymbol{y}(x_0) = \boldsymbol{y}_0 \end{matrix} \right. \implies \boldsymbol{y}=\boldsymbol{y}(x)$$

对于高阶常微分方程，总可以通过降阶化为一阶常微分方程组，例如牛顿第三定律

$$m \frac{\mathrm{d}^2 x}{\mathrm{d} t^2}=F(x, t) \iff \begin{cases} \frac{\mathrm{d} x}{\mathrm{d} t} = \frac{p}{m} \\ \frac{\mathrm{d}p}{\mathrm{d} t} = F(x,t) \end{cases} \iff \dfrac{\mathrm{d}\boldsymbol{y}}{\mathrm{d}t} = f(t,\boldsymbol{y})$$

其中 $\boldsymbol{y}(t)=[x(t),p(t)], \ f(t,\boldsymbol{y}) = [p(t)/m, F(x,t)]$ 。因此，关键在于求解一阶非线性常微分方程。

Taylor级数方法

Taylor级数法适用于 $f(x,y)$ 的解析表达式已知且简单，只是方程非线性导致需要数值计算。方法就是将区间离散化，对 $y(x)$ 求Taylor级数然后迭代。设区间步长为 $h=x_{i+1}-x_i$ ，则有

$$y(x_{i+1}) = y(x_i) + h \frac{\mathrm{d}y}{\mathrm{d}x}\bigg|_{x_i} + \frac{h^2}{2} \frac{\mathrm{d}^2 y}{\mathrm{d}x^2}\bigg|_{x_i} + \cdots$$

其中

$$\frac{\mathrm{d}y}{\mathrm{d}x}\bigg|_{x_i} = f(x_i,y_i)$$

$$\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}\bigg|_{x_i} = \frac{\mathrm{d}}{\mathrm{d}x}f(x,y(x))\bigg|_{x_i} = \left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} \frac{\mathrm{d}y}{\mathrm{d}x} \right)\bigg|_{x_i} = \left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} f\right) \bigg|_{x_i}$$

数值方法综述

其最基本的数值方法是，将第一个方程的 $\mathrm{d}x$ 移到右边去，在 $[x_i,x_{i+1}]$ 上积分得到

$$y_{i+1} - y_i = \int_{x_i}^{x_{i+1}} f(x,y(x)) \,\mathrm{d}x$$

这样就能得到从已知的 $x_i$ 处的函数值 $y_i$ ，推得下一点 $x_{i+1}$ 的函数值 $y_{i+1}$ 。

但问题在于 $y(x)$ 在区间 $[x_i,x_{i+1}]$ 上的表达式未知，需要用合适的方法进行估计，从而给出被积函数 $f(x,y(x))$ 的表达式。我们可以利用前几点的函数值 $(x_{i-1},y_{i-1}),(x_{i-2},y_{i-2}),\cdots$ 作线性外插拟合 $[x_i,x_{i+1}]$ 上 $y(x)$ 的表达式，这就是所谓的「多步法」，在此不介绍。

或者，也可以利用中值定理近似求积分。根据中值定理，必然存在某点 $(x_i+\theta h,y_i + \theta h)$ ，使得被积函数可用该点的函数值代替

$$\int_{x_i}^{x_{i+1}} f(x,y(x)) \,\mathrm{d}x = h f(x_i+\theta h,y_i + \theta h)$$

如何找出最优的 $f(x_i+\theta h,y_i + \theta h)$ 呢？这就是所谓的Runge-Kutta法。

Euler法

最简单的方法是直接用起点 $x_i$ 处函数斜率 $K_1=f(x_i,y_i)$ 去当中值，也就是 $\theta=0$ ，这就是Euler法：

$$y_{i+1} - y_i = h f(x_i,y_i)$$

其几何意义是用斜率为 $K_1=f(x_i,y_i)$ 的斜线拟合了 $y(x)$ ，即设 $y(x)$ 在 $[x_i,x_{i+1}]$ 上的表达式是 $y(x) = K_1 (x - x_i)+y_i$ 。Euler法的精度是 $\omicron(h^2)$ ，因为 $y(x)$ 在 $x_i$ 处的Taylor级数为 $y(x_{i+1}) = y(x_i) + h f(x_i,y_i) + \omicron(h^2)$ 。

Euler法的MATLAB代码为

function y = odeEuler(f,tspan,y0,n)
% odeEuler.m
dim = length(y0);  % dimension
h = (tspan(2) - tspan(1)) / (n-1); % step length
y = zeros(dim, n);
y(:, 1) = y0(:); % initial value
t = linspace(tspan(1), tspan(2), n);

for i = 1:n-1
    K = f(t(i), y(:,i));
    y(:,i+1) = y(:,i) + h*K;
end
end

Euler法可以改进为梯形公式

梯形公式的MATLAB代码为

function y = odeTrapz(f,tspan,y0,n)
% odeEuler.m
dim = length(y0);  % dimension
h = (tspan(2) - tspan(1)) / (n-1); % step length
y = zeros(dim, n);
y(:, 1) = y0(:); % initial value
t = linspace(tspan(1), tspan(2), n);

for i = 1:n-1
    K1 = f(t(i), y(:,i));
    K2 = f(t(i)+h, y(:,i)+h*K1);
    y(:,i+1) = y(:,i) + h/2*(K1+K2);
end
end

Runge-Kutta方法

Runge-Kutta法通过求不同节点上的斜率 $f(x_i+\theta_1 h,y_i + \theta_1 h),f(x_i+\theta_2 h,y_i + \theta_2 h),\cdots$ 作线性组合，构造中值 $K$ 的近似公式，把近似公式和解的Taylor展开式相比较，使其前面若干项相吻合，从而得到相应的高阶方法。其中，节点处的函数值 $y_i + \theta_n h$ 以及相应的 $f(x_i+\theta_n h,y_i + \theta_n h)$ 由前面 $n-1$ 个节点处的斜率 $K_1,K_2,\cdots$ 线性组合后作为区间平均斜率得到的一次函数近似而得。

例如：已有 $K_1$，欲求区间中点处的函数值，可得 $y_{1/2} = y_i + \frac{h}{2} K_1$；已有 $K_1,K_2$，欲求 $2/3$ 节点处的函数值，可设为 $y_{2/3} = y_i + \frac{2}{3}h (\lambda_1 K_1 +\lambda_2 K_2)$ ，其中 $\lambda_1$ 和 $\lambda_2$ 是待定系数。

二阶Runge-Kutta法使用区间中点处的斜率 $K_2$ 作中值，表达式为

$$\begin{cases} y_{i+1} = y_i + h K_2 + \omicron(h^3) \\ K_1 = f(x_i,y_i) \\ K_2 = f(x_i+\frac{1}{2}h,y_i + \frac{1}{2}hK_1) \end{cases}$$

四阶Runge-Kutta法由区间起点、中点、终点处的斜率组合而成，其中中点斜率经过两次修正分别得到 $K_2,K_3$ ，均要使用

$$\begin{cases} y_{i+1} = y_i + h \cdot\frac{1}{6} (K_1 + 2 K_2 + 2K_3 +K_4) + \omicron(h^5) \\ K_1 = f(x_i,y_i) \\ K_2 = f(x_i+\frac{1}{2}h,y_i + \frac{1}{2}hK_1) \\ K_3 = f(x_i+\frac{1}{2}h,y_i + \frac{1}{2}hK_2) \\ K_4 = f(x_i+h,y_i + h K_3) \end{cases}$$

经典四阶Runge-Kutta法的MATLAB代码为

function y = odeRK4(f,tspan,y0,n)
% odeRK4.m
dim = length(y0);  % dimension
h = (tspan(2) - tspan(1)) / (n-1); % step length
y = zeros(dim, n);
y(:, 1) = y0(:); % initial value
t = linspace(tspan(1), tspan(2), n);

for i = 1:n-1
    K1 = f(t(i), y(:,i));
    K2 = f(t(i)+h/2, y(:,i)+K1*h/2);
    K3 = f(t(i)+h/2, y(:,i)+K2*h/2);
    K4 = f(t(i)+h, y(:,i)+K2*h/2);
    y(:,i+1) = y(:,i) + h/6 * (K1+2*K2+2*K3+K4);
end
end

Numerov算法

Numerov算法仅适用于求解不含一阶导数项的二阶常微分方程初值问题

$$\begin{cases} \displaystyle \frac{\mathrm{d}^2 y}{\mathrm{~d} x^2}+k(x) y=S(x) \\ \,y(x_0)=y_0 , \, y^\prime (x_0)= K \end{cases}$$

其迭代公式为

$$y_{n+1}\left(1+\frac{h^2}{12} k_{n+1}\right)=2 y_n\left(1-\frac{5 h^2}{12} k_n\right)-y_{n-1}\left(1+\frac{h^2}{12} k_{n-1}\right)+\frac{h^2}{12}\left(S_{n+1}+10 S_n+S_{n-1}\right)+\mathcal{O}\left(h^6\right)$$

Numerov算法是二步法，需要2个启动值。我们已经有一个 $y_0$ 作为启动值，第二个启动值可利用Taylor展开获得

$$y_1 = y_0 +h y^\prime (x_0)=y_0+Kh$$

MATLAB代码为

function y = numerov(kfun,Sfun,xspan,y0,K,n)
% numerov.m
% y satisfies y'' + k(x)*y = S(x)

dim = length(y0);  % dimension
h = (xspan(2) - xspan(1)) / (n-1); % step length
y = zeros(dim, n);
x = linspace(xspan(1), xspan(2), n);
y(:, 1) = y0(:); % first initial value
y(:, 2) = y0(:)+K*h; % second initial value
k = arrayfun(kfun,x);
S = arrayfun(Sfun,x);

for i = 2:n-1
    y(:,i+1) = (2*y(:,i)*(1-5*h*h/12*k(i))-y(:,i-1)*(1+h*h/12*k(i-1))+h*h/12*(S(:,i+1)+10*S(:,i)+S(:,i-1)))/(1+h*h/12*k(i+1));
end

end

常微分方程的边值问题

二阶线性常微分方程

二阶线性常微分方程是一般二阶ODE的一个特例，其普遍形式为

$$\begin{cases} y^{\prime \prime}+p(x) y^{\prime}+q(x) y=f(x) \\ y(a)=\alpha, \,y(b)=\beta \end{cases}$$

可转化为两个初值问题的叠加

$$\begin{cases} y^{\prime \prime}_1 +p(x) y^{\prime}_1 +q(x) y_1 =f(x) \\ y_1 (a)=\alpha, \,y^\prime_1 (a)=0 \end{cases} \quad + \quad \begin{cases} y^{\prime \prime}_2 +p(x) y^{\prime}_2 +q(x) y_2 =0\\ y_2 (a)=0, \, y^\prime_2 (a)=1 \end{cases}$$

方程的解为

$$y(x)=y_1(x)+\frac{\beta-y_1(b)}{y_2(b)} y_2(x)$$

其中 $y_2(b)\neq 0$ 。

打靶法

第一类边值条件

第一类边值条件(Dirichlet condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime \prime} = f(x,y,y^\prime) \\ y(a)=\alpha, \,y(b)=\beta \end{cases}$$

打靶法的基本思想是：遍历可能的导数初值 $y^\prime(a) = s$ ，求解对应的初值问题

$$\begin{cases} y^{\prime \prime} = f(x,y,y^\prime) \\ y(a)=\alpha, \,y^\prime(a) = s \end{cases}$$

直到寻找到一个 $s=s_0$ ，使得初值问题的解满足 $y(b)=\beta$ 为止。即寻找方程 $y(b,s)-\beta = 0$ 关于变量 $s$ 的根，其中 $y$ 是满足初值 $y(a)=\alpha, \,y^\prime(a)=s$ 的函数。

在具体实现方法上，可应用弦割法求解。因此打靶法求解第一类边值条件的步骤为

1. 给定初始点的斜率猜测值 $s_1$；
2. 用常微分方程的初值问题解法求解 $y(b,s_1)$ ；
3. 给出另一个斜率猜测值 $s_2$ 作为弦割法的第二个启动点;
4. 用常微分方程的初值问题解法求解 $y(b,s_2)$ ；
5. 将 $s_1,s_2$ 作为启动点，弦割法求解方程 @(s) y(b,s) - beta = 0

第二类边值条件

第二类边值条件(Neumann condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime \prime} = f(x,y,y^\prime) \\ y^\prime(a)=\alpha, \,y^\prime(b)=\beta \end{cases}$$

此时靶心是 $y^\prime(b)=\beta$ ，子弹是 $y(a)$ ，待求解的初值问题是

$$\begin{cases} y^{\prime \prime} = f(x,y,y^\prime) \\ y(a)=s, \,y^\prime(a) = \alpha \end{cases}$$

遍历所有可能的 $s$ 值，直到找到 $s=s_0$ 使得 $y^\prime(b,s)-\beta = 0$ 。但要注意，此处的 $y^\prime(b,s)$ 需通过差分来计算

$$y^\prime(b,s) = \frac{y^\prime(b,s)-y^\prime(b-h,s)}{h}$$

第三类边值条件

第三类边值条件(Robin condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime \prime} = f(x,y,y^\prime) \\ \alpha_1 y(a) + \beta_1 y^\prime(a) = r_1 \\ \alpha_2 y(b) + \beta_2 y^\prime(b) = r_2 \end{cases}$$

注意到求解由 $y^{\prime \prime} = f(x,y,y^\prime)$ 描述的初值问题时，常常利用降阶法令 $z_1 = y, z_2 = y^\prime$，这时二阶ODE将会化成关于向量 $\boldsymbol{z}$ 的一阶ODE，进而可应用Runge-Kutta方法求解，同时初值条件也将被转写为

$$\begin{cases} \alpha_1 z_1(a) +\beta_1 z_2(a) = r_1 \\ \alpha_2 z_1(b) +\beta_2 z_2(b) = r_2 \end{cases}$$

此时子弹为 $\boldsymbol{z}(a) = \left[s,\dfrac{r_1-\alpha_1 s}{\beta_1}\right]^T$，$s$ 为可调参数，靶心为 $\alpha_2 z_1(b) +\beta_2 z_2(b) = r_2$ 。

差分法

差分法适用于线性边值问题，即要求微分方程是线性的，形如 $y^{\prime\prime} + p(x) y^{\prime} + q(x) y = f(x)$ 。

第一类边值条件

第一类边值条件(Dirichlet condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime\prime} + p(x) y^{\prime} + q(x) y = f(x) \\ y(a)=\alpha, \,y(b)=\beta \end{cases}$$

将待求解区间 $[a,b]$ 离散化为 $n$ 个节点，$y_1 = y(a), \, y_n = y(b)$ ，即向量的第 $1$ 个和第 $n$ 个分量为边界。

将二阶导数和一阶导数分别差分为

$$\frac{y_{k+1}-2 y_k+y_{k-1}}{h^2}+p_k \frac{y_{k+1}-y_{k-1}}{2 h}+q_k y_k=f_k \, , \, k=2,\cdots,n-1$$

$$\left(2-h p_k\right) y_{k-1}+\left(2 h^2 q_k-4\right) y_k+\left(2+h p_k\right) y_{k+1}=2 h^2 f_k \, , \, k=2,\cdots,n-1$$

结合边界条件，可得矩阵方程组

$$\begin{bmatrix} 2 h^2 q_2-4 & 2+h p_2 & & & \\ 2-h p_3 & 2 h^2 q_3-4 & 2+h p_3 & & \\ & \ddots & \ddots & \ddots & \\ & & 2-h p_{n-2} & 2 h^2 q_{n-2}-4 & 2+h p_{n-2} \\ & & & 2-h p_{n-1} & 2 h^2 q_{n-1}-4 \end{bmatrix}\begin{bmatrix} y_2 \\ y_3 \\ \vdots \\ y_{n-2} \\ y_{n-1} \end{bmatrix}=\begin{bmatrix} 2 h^2 f_1-\left(2-h p_1\right) \alpha \\ 2 h^2 f_2 \\ \vdots \\ 2 h^2 f_{n-2} \\ 2 h^2 f_{n-1}-\left(2+h p_{n-1}\right) \beta \end{bmatrix}$$

第二类边值条件

第二类边值条件(Neumann condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime\prime} + p(x) y^{\prime} + q(x) y = f(x) \\ y^\prime(a)=\alpha, \,y^\prime(b)=\beta \end{cases}$$

离散化后差分方程仍然为

$$\left(2-h p_k\right) y_{k-1}+\left(2 h^2 q_k-4\right) y_k+\left(2+h p_k\right) y_{k+1}=2 h^2 f_k \, , \, k=2,\cdots,n-1$$

离散化后边界条件为

$$y_2 - y_1 = h \alpha \, , \, y_n - y_{n-1} = h \beta$$

可得矩阵方程组

$$\begin{bmatrix} -1 & 1 & & & & \\ 2-h p_2 & 2 h^2 q_2-4 & 2+h p_2 & & & \\ & 2-h p_3 & 2 h^2 q_3-4 & 2+h p_3 & & \\ & & \ddots & \ddots & \ddots & \\ & & 2-h p_{n-2} & 2 h^2 q_{n-2}-4 & 2+h p_{n-2} & \\ & & & 2-h p_{n-1} & 2 h^2 q_{n-1}-4 & 2+h p_{n-1} \\ & & & & -1 & 1 \end{bmatrix}\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-2} \\ y_{n-1} \\ y_n \end{bmatrix}=\begin{bmatrix} h \alpha \\ 2 h^2 f_1 \\ 2 h^2 f_2 \\ \vdots \\ 2 h^2 f_{n-2} \\ 2 h^2 f_{n-1} \\ h \beta \end{bmatrix}$$

第三类边值条件

第三类边值条件(Robin condition)下的二阶非线性常微分方程可以写作

$$\begin{cases} y^{\prime\prime} + p(x) y^{\prime} + q(x) y = f(x) \\ \alpha_1 y(a) + \beta_1 y^\prime(a) = r_1 \\ \alpha_2 y(b) + \beta_2 y^\prime(b) = r_2 \end{cases}$$

边界条件可差分为

$$\alpha_1 y_1 + \beta_1 \frac{y_2 - y_1}{h} = r_1 \, , \,\alpha_2 y_n + \beta_2 \frac{y_n - y_{n-1}}{h} = r_2$$

整理得到

$$\begin{cases} (\alpha_1 h - \beta_1) y_1 + \beta_1 y_2 = r_1 h \\ -\beta_2 y_{n-1} + (\alpha_2 h + \beta_2 ) y_n = r_2 h \end{cases}$$

可得矩阵方程组

$$\begin{bmatrix} \alpha_1 h - \beta_1 & \beta_1 & & & & \\ 2-h p_2 & 2 h^2 q_2-4 & 2+h p_2 & & & \\ & 2-h p_3 & 2 h^2 q_3-4 & 2+h p_3 & & \\ & & \ddots & \ddots & \ddots & \\ & & 2-h p_{n-2} & 2 h^2 q_{n-2}-4 & 2+h p_{n-2} & \\ & & & 2-h p_{n-1} & 2 h^2 q_{n-1}-4 & 2+h p_{n-1} \\ & & & & -\beta_2 & \alpha_2 h + \beta_2 \end{bmatrix}\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-2} \\ y_{n-1} \\ y_n \end{bmatrix}=\begin{bmatrix} r_1 h \\ 2 h^2 f_1 \\ 2 h^2 f_2 \\ \vdots \\ 2 h^2 f_{n-2} \\ 2 h^2 f_{n-1} \\ r_2 h \end{bmatrix}$$

偏微分方程的边值问题

二阶椭圆型偏微分方程

二阶椭圆型PDE形如

$$-\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right) \phi=S(x, y)$$

求解区域为长方形 $0 \le x \le a, 0\le y \le b$ ，辅以边界条件

$$\begin{gathered} \phi(x,0) = f_1(x) \quad \phi(x,b) = f_2(x) \\ \phi(0,y) = g_1(y) \quad \phi(a,y) = g_2(y) \end{gathered}$$

为便于计算 $x$ 和 $y$ 方向应当选取相同差分步长 $h$ ，这样二阶椭圆型PDE的差分形式可以写作

$$-\left(\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{h^2} + \frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{h^2}\right) = S_{ij}$$

$$\phi_{ij} = \frac{1}{4} \left(\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}+h^2 S_{ij}\right)$$

Jacobi迭代法

Jacobi迭代法的基本思路是：猜测一个初始解 $\phi^0$，第 $n$ 次迭代解 $\phi_{ij}^n$ 由前一次 $\phi_{ij}^{n-1}$ 迭代解代入PDE的差分方程计算得到，直到收敛为止，即前后两次迭代解的差别不大：$|\phi_{ij}^n - \phi_{ij}^{n-1}| < \varepsilon$ 。

$$\phi_{ij}^n = \frac{1}{4} \left(\phi_{i+1,j}^{n-1}+\phi_{i-1,j}^{n-1}+\phi_{i,j+1}^{n-1}+\phi_{i,j-1}^{n-1}+h^2 S_{ij}\right)$$

Gauss-Seidel迭代法

Gauss-Seidel迭代法在Jacobi迭代法的基础上，在计算第 $n$ 次迭代解的 $(i,j)$ 元时，将周围已经算出的 $i-1$ 和 $j-1$ 元，从使用第 $n-1$ 次的结果替换为使用第 $n$ 次的结果。

$$\phi_{ij}^n = \frac{1}{4} \left(\phi_{i+1,j}^{n-1}+\phi_{i-1,j}^{n}+\phi_{i,j+1}^{n-1}+\phi_{i,j-1}^{n}+h^2 S_{ij}\right)$$

松弛法

松弛法将第 $n$ 次迭代解的值取为新值和旧值的线性组合，即

$$\phi_{ij}^n = \omega \tilde{\phi}_{ij}^n + (1-\omega) \phi_{ij}^n$$

其中 $\tilde{\phi}_{ij}^n = \frac{1}{4} \left(\phi_{i+1,j}^{n-1}+\phi_{i-1,j}^{n}+\phi_{i,j+1}^{n-1}+\phi_{i,j-1}^{n}+h^2 S_{ij}\right)$ 是本应当取的新值，$\omega$ 是松弛因子，取值范围 $0<\omega<2$ 。

• 当 $0<\omega<1$ 时称为低松弛，新的迭代解中旧值占主导，收敛速度变慢，计算稳定；
• 当 $1<\omega<2$ 时称为超松弛，新的迭代解中新值占主导，收敛速度变快，计算加速。

因此Gauss-Seidel迭代法结合松弛法的迭代公式为

$$\phi_{ij}^{n+1} = \frac{\omega}{4} \left(\phi_{i+1,j}^{n}+\phi_{i-1,j}^{n+1}+\phi_{i,j+1}^{n}+\phi_{i,j-1}^{n+1}+h^2 S_{ij}\right) +(1-\omega)\phi_{ij}^{n}$$

二阶抛物型偏微分方程

抛物型PDE形如

$$\frac{\partial \phi}{\partial t} = \nabla \cdot (D \nabla \phi) + S$$

对二阶抛物型PDE则有

$$\frac{\partial \phi}{\partial t} = D \frac{\partial^2 \phi}{\partial x^2} + S(x,t)$$

将空间 $x$ 间隔 $h$ 差分为 $\left\{x_i\right\} , i=1,2\cdots,N$，$x_1$ 为左端点，$x_N$ 为右端点。

时间 $t$ 间隔 $\Delta t$ 差分为 $\left\{t_j\right\} , j=1,2\cdots,M$，$t_1$ 为起始时刻，$t_M$ 为终止时刻。

显式差分法

显式差分法中，$\frac{\partial^2 \phi}{\partial x^2}$ 使用当前时刻 $j$ 的解 $\phi^j$ 进行差分，即

$$\frac{\phi_i^{j+1}-\phi_i^j}{\Delta t} = D \frac{\phi_{i+1}^j-2\phi_i^j+\phi_{i-1}^j}{h^2} + S_i^j$$

$$\phi_i^{j+1} = \phi_i^j + \frac{D \Delta t}{h^2} \left(\phi_{i+1}^j-2\phi_i^j+\phi_{i-1}^j\right) + S_i^j \Delta t$$

定义 $r = \frac{D \Delta t}{h^2}$，由数值分析可得 $r<1/2$ 时算法稳定。当我们选用更小的空间步长 $h$ 时，迫使我们不得不取一个更小的时间步长 $\Delta t$，这将付出很大的代价，特别是这个时间步长比恰当地描述系统演化所需要的自然时间尺度要小的多时。因此显式差分法需要改进。

隐式差分法

隐式差分法中，$\frac{\partial^2 \phi}{\partial x^2}$ 使用下一时刻 $j+1$ 的解 $\phi^{j+1}$ 进行差分，即

$$\frac{\phi_i^{j+1}-\phi_i^j}{\Delta t} = D \frac{\phi_{i+1}^{j+1}-2\phi_i^{j+1}+\phi_{i-1}^{j+1}}{h^2} + S_i^j$$

$$\phi_i^{j+1} = \phi_i^j + \frac{D \Delta t}{h^2} \left(\phi_{i+1}^{j+1}-2\phi_i^{j+1}+\phi_{i-1}^{j+1}\right) + S_i^j \Delta t$$

这是一个关于 $\phi^{j+1}_i$ 的三对角方程。若定义 $r = \frac{D \Delta t}{h^2}$，由数值分析可得隐式算法无条件稳定。方程可写为

$$-r\phi_{i-1}^{j+1} + (1+2r)\phi_i^{j+1} -r\phi_{i+1}^{j+1}= \phi_i^j + S_i^j \Delta t \ (i=2,3,\cdots,N-1)$$

对于左边界 $\phi_1^j$ 和右边界 $\phi_N^j$ ，它们由给定的边值条件决定，因此待求解的是一个 $N-2$ 维的三对角方程

$$\begin{bmatrix} 1+2r & -r & & & \\ -r & 1+2r & -r & & \\ & \ddots & \ddots & \ddots & \\ & & -r & 1+2r & -r \\ & & & -r & 1+2r \end{bmatrix}\begin{bmatrix} \phi^{j+1}_2 \\ \phi^{j+1}_3 \\ \vdots \\ \phi^{j+1}_{N-2} \\ \phi^{j+1}_{N-1} \end{bmatrix}=\begin{bmatrix} \phi^j_2 + S_2^j\Delta t + r \phi^{j+1}_1\\ \phi^j_3 + S_3^j \Delta t\\ \vdots \\ \phi^j_{N-2} + S_{N-2}^j \Delta t \\ \phi^j_{N-1} + S_{N-1}^j \Delta t + r\phi^{j+1}_N \end{bmatrix}$$

Crank-Nicolson算法

Crank-Nicolson算法又称为平均隐式差分法，$\frac{\partial^2 \phi}{\partial x^2}$ 使用当前时刻 $j$ 的解 $\phi^j$ 和下一时刻 $j+1$ 的解 $\phi^{j+1}$ 之平均进行差分

$$\frac{\phi_i^{j+1}-\phi_i^j}{\Delta t} = \frac{D}{2} \left(\frac{\phi_{i+1}^j-2\phi_i^j+\phi_{i-1}^j}{h^2} + \frac{\phi_{i+1}^{j+1}-2\phi_i^{j+1}+\phi_{i-1}^{j+1}}{h^2} \right)+ S_i^j$$

$$-\frac{D \Delta t}{2h^2} \phi_{i-1}^{j+1}+\left( 1+\frac{D \Delta t}{h^2} \right) \phi_{i}^{j+1} - \frac{D \Delta t}{2h^2}\phi_{i+1}^{j+1} = \frac{D \Delta t}{2h^2} \phi_{i+1}^{j}+\left( 1-\frac{D \Delta t}{h^2} \right) \phi_{i}^{j} + \frac{D \Delta t}{2h^2}\phi_{i-1}^{j} + S_i^j \Delta t$$

若仍定义 $r = \frac{D \Delta t}{h^2}$，由数值分析可得隐式算法无条件稳定。方程可写为

$$-\frac{r}{2} \phi_{i-1}^{j+1}+(1+r) \phi_{i}^{j+1} - \frac{r}{2}\phi_{i+1}^{j+1} = \frac{r}{2} \phi_{i+1}^{j}+(1-r)\phi_{i}^{j} + \frac{r}{2}\phi_{i-1}^{j} + S_i^j \Delta t$$

对于左边界 $\phi_1^j$ 和右边界 $\phi_N^j$ ，它们由给定的边值条件决定，因此待求解的是一个 $N-2$ 维的三对角方程

$$\begin{bmatrix} 1+r & -r/2 & & & \\ -r/2 & 1+r & -r/2 & & \\ & \ddots & \ddots & \ddots & \\ & & -r/2 & 1+r & -r/2 \\ & & & -r/2 & 1+r \end{bmatrix}\begin{bmatrix} \phi^{j+1}_2 \\ \phi^{j+1}_3 \\ \vdots \\ \phi^{j+1}_{N-2} \\ \phi^{j+1}_{N-1} \end{bmatrix}=\begin{bmatrix} \frac{r}{2}\phi^j_3 + (1-r)\phi^j_2+ \frac{r}{2}\phi^j_1+ S_2^j\Delta t + \frac{r}{2} \phi^{j+1}_1\\ \frac{r}{2}\phi^j_4 + (1-r)\phi^j_3+ \frac{r}{2}\phi^j_2 + S_3^j \Delta t\\ \vdots \\ \frac{r}{2}\phi^j_{N-1} + (1-r)\phi^j_{N-2}+ \frac{r}{2}\phi^j_{N-3} + S_{N-2}^j \Delta t \\ \frac{r}{2}\phi^j_{N} + (1-r)\phi^j_{N-1}+ \frac{r}{2}\phi^j_{N-2} + S_{N-1}^j \Delta t + \frac{r}{2} \phi^{j+1}_{N} \end{bmatrix}$$

二阶双曲型偏微分方程

二阶双曲型PDE形如

$$\frac{\partial u^2}{\partial t^2} - a^2 \frac{\partial^2 u}{\partial x^2} = S(x,t)$$

求解区域为 $a \le x \le b$ ，辅以边界条件

$$\begin{gathered} u|_{t=0} = \varphi(x) \quad \frac{\partial u}{\partial t}\bigg|_{t=0} = \psi(x) \\ u(a,t) = f(t) \quad u(b,t) = g(t) \end{gathered}$$

隐式差分法

隐式差分法中，$\frac{\partial^2 u}{\partial x^2}$ 使用下一时刻 $j+1$ 的解 $u^{j+1}$ 进行差分，即

$$\frac{u_i^{j+1}-2u_i^j+u_i^{j-1}}{(\Delta t)^2} - a^2 \frac{u_{i+1}^{j+1}-2u_i^{j+1}+u_{i-1}^{j+1}}{h^2} = S_i^j$$

若仍定义 $r = \frac{a^2 (\Delta t)^2}{h^2}$，由数值分析可得隐式算法无条件稳定。方程可写为

$$-ru_{i-1}^{j+1} + (1+2r)u_i^{j+1} -ru_{i+1}^{j+1}= 2u_i^j - u_i^{j-1}+ S_i^j (\Delta t)^2$$

对于左边界 $u_1^j$ 和右边界 $u_N^j$ ，它们由给定的边值条件决定，因此待求解的是一个 $N-2$ 维的三对角方程

$$\begin{bmatrix} 1+2r & -r & & & \\ -r & 1+2r & -r & & \\ & \ddots & \ddots & \ddots & \\ & & -r & 1+2r & -r \\ & & & -r & 1+2r \end{bmatrix}\begin{bmatrix} u^{j+1}_2 \\ u^{j+1}_3 \\ \vdots \\ u^{j+1}_{N-2} \\ u^{j+1}_{N-1} \end{bmatrix}=\begin{bmatrix} 2u^j_2 -u_2^{j-1}+ S_2^j(\Delta t)^2 + r u^{j+1}_1\\ 2u^j_3 -u_3^{j-1}+ S_3^j(\Delta t)^2 \\ \vdots \\ 2u^j_{N-2} -u_{N-2}^{j-1}+ S_{N-2}^j(\Delta t)^2 \\ 2u^j_{N-1} -u_{N-1}^{j-1}+ S_{N-1}^j(\Delta t)^2 + ru^{j+1}_N \end{bmatrix}$$