计算物理(一)基础运算与微分方程求解
fengxiaot Lv4

本篇文章是西安交大《计算物理与程序设计》课程的笔记,内容比较散乱,主打作为一个考前速查手册。

MATLAB基础语法

标量向量矩阵

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

函数

匿名函数:f = @(x) expr(x) ,这样创建的 f 是函数句柄 (function bundle) 。

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


数值积分

内置命令

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

1
I = integral(f,xmin,xmax)

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

类似的二维数值积分和三维数值积分命令为 integral2integral3 ,语法为

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

梯形公式

以下讨论中,总是将积分区间 (a,b)(a,b) 切割为 NN 份,并设初始节点 x1=ax_1=a,末端节点 xN+1=bx_{N+1}=b,积分步长 hh

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

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

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

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

Itrapz=h2f(a)+hi=2Nf(xi)+h2f(b)I_{\text{trapz}} = \frac{h}{2} f(a)+ h \sum_{i=2}^N f(x_i) +\frac{h}{2} f(b)

MATLAB代码为

1
2
3
4
5
6
7
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公式来源于二次插值,涉及连续三个节点,三个节点夹着两个区间,因此 NN 必须是 22 的倍数。

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

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

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

Isimpson=h3f(a)+4h3i=2sep=2Nf(xi)+2h3i=3sep=2N1f(xi)+h2f(b)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代码为

1
2
3
4
5
6
7
8
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)x0x_0 开始,尝试求出 f(x)=0f(x) = 0 的解并赋给 xx

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

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

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

{e(x1+x2)x2(1+x12)=0x1cos(x2)+x2sin(x1)12=0\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}

1
2
3
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;

二分法

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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(x) 的导数 f(x)f^\prime(x) 的解析表达式,并给定一个启动值 x0x_0 。Newton - Raphson 迭代公式为

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

xi+1xi<ε| x_{i+1} - x_i | <\varepsilonf(xi)<δf(x_i)<\delta 时视为找到了根。

MATLAB代码为

1
2
3
4
5
6
7
8
9
10
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

弦割法

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

xi+1=xif(xi)xixi1f(xi)f(xi1)x_{i+1} = x_i - f(x_i) \frac{x_i - x_{i-1}}{f(x_i) - f(x_{i-1})}

1
2
3
4
5
6
7
8
9
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

预搜索弦割

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

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

1
2
3
4
5
6
7
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

常微分方程的初值问题

内置命令

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

{dydx=f(x,y)y(x0)=y0    y=y(x)\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)

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

{dydx=f(x,y)y(x0)=y0    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)

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

md2xdt2=F(x,t)    {dxdt=pmdpdt=F(x,t)    dydt=f(t,y)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})

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

Taylor级数方法

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

y(xi+1)=y(xi)+hdydxxi+h22d2ydx2xi+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

其中

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

d2ydx2xi=ddxf(x,y(x))xi=(fx+fydydx)xi=(fx+fyf)xi\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}

数值方法综述

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

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

这样就能得到从已知的 xix_i 处的函数值 yiy_i ,推得下一点 xi+1x_{i+1} 的函数值 yi+1y_{i+1}

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

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

xixi+1f(x,y(x))dx=hf(xi+θh,yi+θ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(xi+θh,yi+θh)f(x_i+\theta h,y_i + \theta h) 呢?这就是所谓的Runge-Kutta法。

Euler法

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

yi+1yi=hf(xi,yi)y_{i+1} - y_i = h f(x_i,y_i)

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

Euler法的MATLAB代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
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代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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(xi+θ1h,yi+θ1h),f(xi+θ2h,yi+θ2h),f(x_i+\theta_1 h,y_i + \theta_1 h),f(x_i+\theta_2 h,y_i + \theta_2 h),\cdots 作线性组合,构造中值 KK 的近似公式,把近似公式和解的Taylor展开式相比较,使其前面若干项相吻合,从而得到相应的高阶方法。其中,节点处的函数值 yi+θnhy_i + \theta_n h 以及相应的 f(xi+θnh,yi+θnh)f(x_i+\theta_n h,y_i + \theta_n h) 由前面 n1n-1 个节点处的斜率 K1,K2,K_1,K_2,\cdots 线性组合后作为区间平均斜率得到的一次函数近似而得。

例如:已有 K1K_1,欲求区间中点处的函数值,可得 y1/2=yi+h2K1y_{1/2} = y_i + \frac{h}{2} K_1;已有 K1,K2K_1,K_2,欲求 2/32/3 节点处的函数值,可设为 y2/3=yi+23h(λ1K1+λ2K2)y_{2/3} = y_i + \frac{2}{3}h (\lambda_1 K_1 +\lambda_2 K_2) ,其中 λ1\lambda_1λ2\lambda_2 是待定系数。

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

{yi+1=yi+hK2+ο(h3)K1=f(xi,yi)K2=f(xi+12h,yi+12hK1)\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法由区间起点、中点、终点处的斜率组合而成,其中中点斜率经过两次修正分别得到 K2,K3K_2,K_3 ,均要使用

{yi+1=yi+h16(K1+2K2+2K3+K4)+ο(h5)K1=f(xi,yi)K2=f(xi+12h,yi+12hK1)K3=f(xi+12h,yi+12hK2)K4=f(xi+h,yi+hK3)\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代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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算法仅适用于求解不含一阶导数项的二阶常微分方程初值问题

{d2y dx2+k(x)y=S(x)y(x0)=y0,y(x0)=K\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}

其迭代公式为

yn+1(1+h212kn+1)=2yn(15h212kn)yn1(1+h212kn1)+h212(Sn+1+10Sn+Sn1)+O(h6)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个启动值。我们已经有一个 y0y_0 作为启动值,第二个启动值可利用Taylor展开获得

y1=y0+hy(x0)=y0+Khy_1 = y_0 +h y^\prime (x_0)=y_0+Kh

MATLAB代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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的一个特例,其普遍形式为

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

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

{y1+p(x)y1+q(x)y1=f(x)y1(a)=α,y1(a)=0+{y2+p(x)y2+q(x)y2=0y2(a)=0,y2(a)=1\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)=y1(x)+βy1(b)y2(b)y2(x)y(x)=y_1(x)+\frac{\beta-y_1(b)}{y_2(b)} y_2(x)

其中 y2(b)0y_2(b)\neq 0

打靶法

第一类边值条件

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

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

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

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

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

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

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

第二类边值条件

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

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

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

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

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

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

第三类边值条件

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

{y=f(x,y,y)α1y(a)+β1y(a)=r1α2y(b)+β2y(b)=r2\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=f(x,y,y)y^{\prime \prime} = f(x,y,y^\prime) 描述的初值问题时,常常利用降阶法令 z1=y,z2=yz_1 = y, z_2 = y^\prime,这时二阶ODE将会化成关于向量 z\boldsymbol{z} 的一阶ODE,进而可应用Runge-Kutta方法求解,同时初值条件也将被转写为

{α1z1(a)+β1z2(a)=r1α2z1(b)+β2z2(b)=r2\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}

此时子弹为 z(a)=[s,r1α1sβ1]T\boldsymbol{z}(a) = \left[s,\dfrac{r_1-\alpha_1 s}{\beta_1}\right]^Tss 为可调参数,靶心为 α2z1(b)+β2z2(b)=r2\alpha_2 z_1(b) +\beta_2 z_2(b) = r_2

差分法

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

第一类边值条件

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

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

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

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

yk+12yk+yk1h2+pkyk+1yk12h+qkyk=fk,k=2,,n1\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

(2hpk)yk1+(2h2qk4)yk+(2+hpk)yk+1=2h2fk,k=2,,n1\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

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

[2h2q242+hp22hp32h2q342+hp32hpn22h2qn242+hpn22hpn12h2qn14][y2y3yn2yn1]=[2h2f1(2hp1)α2h2f22h2fn22h2fn1(2+hpn1)β]\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)下的二阶非线性常微分方程可以写作

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

离散化后差分方程仍然为

(2hpk)yk1+(2h2qk4)yk+(2+hpk)yk+1=2h2fk,k=2,,n1\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

离散化后边界条件为

y2y1=hα,ynyn1=hβy_2 - y_1 = h \alpha \, , \, y_n - y_{n-1} = h \beta

可得矩阵方程组

[112hp22h2q242+hp22hp32h2q342+hp32hpn22h2qn242+hpn22hpn12h2qn142+hpn111][y1y2y3yn2yn1yn]=[hα2h2f12h2f22h2fn22h2fn1hβ]\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)下的二阶非线性常微分方程可以写作

{y+p(x)y+q(x)y=f(x)α1y(a)+β1y(a)=r1α2y(b)+β2y(b)=r2\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}

边界条件可差分为

α1y1+β1y2y1h=r1,α2yn+β2ynyn1h=r2\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

整理得到

{(α1hβ1)y1+β1y2=r1hβ2yn1+(α2h+β2)yn=r2h\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}

可得矩阵方程组

[α1hβ1β12hp22h2q242+hp22hp32h2q342+hp32hpn22h2qn242+hpn22hpn12h2qn142+hpn1β2α2h+β2][y1y2y3yn2yn1yn]=[r1h2h2f12h2f22h2fn22h2fn1r2h]\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形如

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

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

ϕ(x,0)=f1(x)ϕ(x,b)=f2(x)ϕ(0,y)=g1(y)ϕ(a,y)=g2(y)\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}

为便于计算 xxyy 方向应当选取相同差分步长 hh ,这样二阶椭圆型PDE的差分形式可以写作

(ϕi+1,j2ϕi,j+ϕi1,jh2+ϕi,j+12ϕi,j+ϕi,j1h2)=Sij-\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}

ϕij=14(ϕi+1,j+ϕi1,j+ϕi,j+1+ϕi,j1+h2Sij)\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迭代法的基本思路是:猜测一个初始解 ϕ0\phi^0,第 nn 次迭代解 ϕijn\phi_{ij}^n 由前一次 ϕijn1\phi_{ij}^{n-1} 迭代解代入PDE的差分方程计算得到,直到收敛为止,即前后两次迭代解的差别不大:ϕijnϕijn1<ε|\phi_{ij}^n - \phi_{ij}^{n-1}| < \varepsilon

ϕijn=14(ϕi+1,jn1+ϕi1,jn1+ϕi,j+1n1+ϕi,j1n1+h2Sij)\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迭代法的基础上,在计算第 nn 次迭代解的 (i,j)(i,j) 元时,将周围已经算出的 i1i-1j1j-1 元,从使用第 n1n-1 次的结果替换为使用第 nn 次的结果。

ϕijn=14(ϕi+1,jn1+ϕi1,jn+ϕi,j+1n1+ϕi,j1n+h2Sij)\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)

松弛法

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

ϕijn=ωϕ~ijn+(1ω)ϕijn\phi_{ij}^n = \omega \tilde{\phi}_{ij}^n + (1-\omega) \phi_{ij}^n

其中 ϕ~ijn=14(ϕi+1,jn1+ϕi1,jn+ϕi,j+1n1+ϕi,j1n+h2Sij)\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<ω<20<\omega<2

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

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

ϕijn+1=ω4(ϕi+1,jn+ϕi1,jn+1+ϕi,j+1n+ϕi,j1n+1+h2Sij)+(1ω)ϕijn\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形如

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

对二阶抛物型PDE则有

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

将空间 xx 间隔 hh 差分为 {xi},i=1,2,N\left\{x_i\right\} , i=1,2\cdots,Nx1x_1 为左端点,xNx_N 为右端点。

时间 tt 间隔 Δt\Delta t 差分为 {tj},j=1,2,M\left\{t_j\right\} , j=1,2\cdots,Mt1t_1 为起始时刻,tMt_M 为终止时刻。

显式差分法

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

ϕij+1ϕijΔt=Dϕi+1j2ϕij+ϕi1jh2+Sij\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

ϕij+1=ϕij+DΔth2(ϕi+1j2ϕij+ϕi1j)+SijΔt\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=DΔth2r = \frac{D \Delta t}{h^2},由数值分析可得 r<1/2r<1/2 时算法稳定。当我们选用更小的空间步长 hh 时,迫使我们不得不取一个更小的时间步长 Δt\Delta t,这将付出很大的代价,特别是这个时间步长比恰当地描述系统演化所需要的自然时间尺度要小的多时。因此显式差分法需要改进。

隐式差分法

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

ϕij+1ϕijΔt=Dϕi+1j+12ϕij+1+ϕi1j+1h2+Sij\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

ϕij+1=ϕij+DΔth2(ϕi+1j+12ϕij+1+ϕi1j+1)+SijΔt\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

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

rϕi1j+1+(1+2r)ϕij+1rϕi+1j+1=ϕij+SijΔt (i=2,3,,N1)-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)

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

[1+2rrr1+2rrr1+2rrr1+2r][ϕ2j+1ϕ3j+1ϕN2j+1ϕN1j+1]=[ϕ2j+S2jΔt+rϕ1j+1ϕ3j+S3jΔtϕN2j+SN2jΔtϕN1j+SN1jΔt+rϕNj+1]\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算法又称为平均隐式差分法,2ϕx2\frac{\partial^2 \phi}{\partial x^2} 使用当前时刻 jj 的解 ϕj\phi^j 和下一时刻 j+1j+1 的解 ϕj+1\phi^{j+1} 之平均进行差分

ϕij+1ϕijΔt=D2(ϕi+1j2ϕij+ϕi1jh2+ϕi+1j+12ϕij+1+ϕi1j+1h2)+Sij\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

DΔt2h2ϕi1j+1+(1+DΔth2)ϕij+1DΔt2h2ϕi+1j+1=DΔt2h2ϕi+1j+(1DΔth2)ϕij+DΔt2h2ϕi1j+SijΔt-\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=DΔth2r = \frac{D \Delta t}{h^2},由数值分析可得隐式算法无条件稳定。方程可写为

r2ϕi1j+1+(1+r)ϕij+1r2ϕi+1j+1=r2ϕi+1j+(1r)ϕij+r2ϕi1j+SijΔt-\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

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

[1+rr/2r/21+rr/2r/21+rr/2r/21+r][ϕ2j+1ϕ3j+1ϕN2j+1ϕN1j+1]=[r2ϕ3j+(1r)ϕ2j+r2ϕ1j+S2jΔt+r2ϕ1j+1r2ϕ4j+(1r)ϕ3j+r2ϕ2j+S3jΔtr2ϕN1j+(1r)ϕN2j+r2ϕN3j+SN2jΔtr2ϕNj+(1r)ϕN1j+r2ϕN2j+SN1jΔt+r2ϕNj+1]\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}