📢该推文部分内容是作者个人猜测，谨慎阅读。欢迎在评论区纠正！
⛏ 作者正在摸鱼，该推文尚未完成

微分方程往往是连续的，可计算机更偏好以离散的视角观察世界。如何用计算机求解微分方程？FDM给出了一种近似计算方法。

理论

$\begin{align} Taylor's\,Polynomial\quad f(x_0+h)&= f(x_0)+\frac{f'(x_0)}{1!}h+\frac{f'(x_0)}{2!}h^2...+\frac{f'(x_0)}{n!}h^n+R_n(x)\\ \Rightarrow f'(x_0)&= \frac{f(x_0+h)-f(x_0)}h-\frac{R_1(x_0+h)}h\\ f'(x_0)&\approx \frac{f(x_0+h)-f(x_0)}h \end{align}$

✋Q：为什么要近似？知道解析式不能直接求出导数吗？

🧐A（个人猜测）：MATLAB的符号函数可以实现求部分导函数解析式的功能了（并不知道怎么实现的）。相比DFM是数值的，可以求解没有导函数，以及函数只有数据点的情况，更加普适。在运算量上，说不定也有优势。

误差分析

Error

round-off error
truncation error

truncation error

$\begin{align} y&\equiv \lim_{h\rightarrow0}\frac{f(x+h)-f(x)}{h}\\ y^\ast&\equiv \frac{f(x+h)-f(x)}h\\ |\Delta y|&=|y^\ast-y|=|\frac{R_1(x+h)}h|\\ Lagrange\;Form\quad R_n(x+h)&= \frac{f^{(n+1)}(\xi)}{(n+1)!}h^{n+1}\\ \Rightarrow|\Delta y|&=|f''(\xi)h|\\ \end{align}$

$h>0,x<\xi<x+h;h<0,x+h<\xi<x$

$x\rightarrow y$ 唯一，而 $x\rightarrow y^\ast$ 和 $h$ 的取值相关。FDM求解的是 $y^\ast$ ，选择的映射方式不同，解的值 $y^\ast$ 也不同。

对比求导表达式，当 $h\rightarrow 0$ 时，收敛到理想值。步长 $h$ 绝对值越小，精度越大，运算量越大。

不同的映射方式分为forward difference, backward difference和central difference。他们的误差都是 $f’’(\xi)h$ 。

$\begin{align} Define\;&h>0\\ Forward\;Difference\quad \bigtriangleup_h [f](x)&\equiv f(x+h)-f(x)\\ Backword\;Difference\quad \bigtriangledown_h [f](x)&\equiv f(x)-f(x-h)\\ Central\;Difference\quad \delta_h [f](x)&\equiv f(x+h/2)-f(x-h/2)\\ &=\bigtriangleup_{h/2} [f](x)+\bigtriangledown_{h/2} [f](x)\\ y^\ast&=\frac{Difference|h}{h} \end{align}$

应用

求解1-D热传导方程

根据表达$\frac{\delta U}{\delta t}$，$\frac{\delta ^2 U}{\delta t^2}$的差分方式，可分为三种方法。

	时间一阶偏微	空间二阶偏微	时间基准
Explicit Method	forward	center	$t_n$
Implicit Method	backward	center	$t_{n+1}$
Crank-Nicolson Method	center	center	$t_{n+1/2}$

MATLAB仿真

$\begin{align} &U_t=\alpha U_{xx},\quad\alpha=\frac1{\pi^2}\\ &U(0,t)=U(1,t)=0\\ &U(x,0)=\frac{1}{\pi^2}\sin(\pi x) \end{align}$

本题的理论解

$U(x,t)=\frac1{\pi^2}e^{-t}\sin(\pi x)$

Explicit Method

Explicit Method最容易实现，计算量最小。

$\frac{u_j^{n+1}-u_j^{n}}{k} =\alpha\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}$ $u_j^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^n+ru_{j+1}^n\quad r=\frac{\alpha k}{h^2}$

循环实现

Alpha = 1/pi^2;
k = 0.1; % step of time
h = 0.1; % step of space
r = Alpha*k/h^2;

t = 0:k:2;
x = 0:h:1;
nt = length(t); % the length of t
nx = length(x); % the length of x

u = zeros(nx,nt); % 初始化的同时设置了边界
u(:,1) = Alpha * sin(pi*x); % 初值条件

for n = 1:nt-1 % loop of time
    for j = 1+1:nx-1 % loop of space
        u(j,n+1)=(1-2*r)*u(j,n)+r*u(j-1,n)+r*u(j+1,n);
    end
end

clf;
mesh(t,x,u);
xlabel('t','FontSize',20);
ylabel('x','FontSize',20);
title({'Explicit Method by Loop'; ...
    [' k = ',num2str(k)];...
    [' h = ',num2str(h)];...
    [' r = ',num2str(r)]})

矩阵实现

$\begin{aligned} u_j^{n+1} &= \begin{bmatrix} r&1-2r&r \end{bmatrix} \begin{bmatrix} u_{j+1}^{n}\\ u_{j}^{n}\\ u_{j-1}^{n} \end{bmatrix} \\ \begin{bmatrix} u_{J-1}^{n+1}\\u_{J-2}^{n+1}\\\vdots\\u_{2}^{n+1} \end{bmatrix} &= \begin{bmatrix} r & 1-2r & r & 0 & 0 & \cdots & 0 \\ 0 & r & 1-2r & r & 0 & \cdots & 0 \\ 0 & 0 & r & 1-2r & r & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \cdots & r \\ \end{bmatrix} \begin{bmatrix} u_{J}^{n}\\u_{J-1}^{n}\\\vdots\\u_{1}^{n} \end{bmatrix}\\ \implies U_{n+1}&=LU_{n} \end{aligned}$

$U_{n+1}$和$U_n$的维度并不相同。删减$L$的最左最右列变成方阵（还是三对角矩阵）。

$\begin{bmatrix} u_{J-1}^{n+1}-ru_{J}^{n}\\u_{J-2}^{n+1}\\\vdots\\u_{2}^{n+1}-ru_{1}^{n} \end{bmatrix} = \begin{bmatrix} 1-2r & r & 0 & 0 & \cdots & 0 \\ r & 1-2r & r & 0 & \cdots & 0 \\ 0 & r & 1-2r & r & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 1-2r \\ \end{bmatrix} \begin{bmatrix} u_{J-1}^{n}\\u_{J-2}^{n}\\\vdots\\u_{2}^{n} \end{bmatrix}\\$

Alpha = 1/pi^2;
k = 0.1; % step of time
h = 0.1; % step of space
r = Alpha*k/h^2;

t = 0:k:2;
x = 0:h:1;
nt = length(t); % the length of t
nx = length(x); % the length of x

u = zeros(nx,nt); % 初始化的同时设置了边界
u(:,1) = Alpha * sin(pi*x); % 初值条件

d = ones(nx-2,1)*(1-2*r);
e = ones(nx-3,1)*r;
L = full(gallery('tridiag',e,d,e));

for n = 1:nt-1 % loop of time
        u(2:nx-1,n+1)=L*u(2:nx-1,n);
end

clf;
mesh(t,x,u);
xlabel('t','FontSize',20);
ylabel('x','FontSize',20);
title('Explicit Method by Matrix')

Implicit Method

$\frac{u_j^{n+1}-u_j^{n}}{k} =\alpha\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}$ $u_j^{n}=(1+2r)u_j^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}\quad r=\frac{\alpha k}{h^2}$ $\begin{aligned} \begin{bmatrix} u_{J-1}^{n}\\u_{J-2}^{n}\\\vdots\\u_{2}^{n} \end{bmatrix} &= \begin{bmatrix} -r & 1+2r & -r & 0 & 0 & \cdots & 0 \\ 0 & -r & 1+2r & -r & 0 & \cdots & 0 \\ 0 & 0 & -r & 1+2r & -r & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \cdots & -r \\ \end{bmatrix} \begin{bmatrix} u_{J}^{n+1}\\u_{J-1}^{n+1}\\\vdots\\u_{1}^{n+1} \end{bmatrix}\\ \begin{bmatrix} u_{J-1}^{n}+ru_{J}^{n+1}\\u_{J-2}^{n}\\\vdots\\u_{2}^{n}+ru_{1}^{n+1} \end{bmatrix} &= \begin{bmatrix} 1+2r & -r & 0 & 0 & \cdots & 0 \\ -r & 1+2r & -r & 0 & \cdots & 0 \\ 0 & -r & 1+2r & -r & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 1+2r \\ \end{bmatrix} \begin{bmatrix} u_{J-1}^{n+1}\\u_{J-2}^{n+1}\\\vdots\\u_{2}^{n+1} \end{bmatrix} \end{aligned}$

Alpha = 1/pi^2;
k = 0.1; % step of time
h = 0.1; % step of space
r = Alpha*k/h^2;

t = 0:k:2;
x = 0:h:1;
nt = length(t); % the length of t
nx = length(x); % the length of x

u = zeros(nx,nt); % 初始化的同时设置了边界
u(:,1) = Alpha * sin(pi*x); % 初值条件

d = ones(nx-2,1)*(1+2*r);
e = ones(nx-3,1)*(-r);
L = full(gallery('tridiag',e,d,e));

for n = 1:nt-1 % loop of time
        u(2:nx-1,n+1)=L\u(2:nx-1,n);
end

clf;
mesh(t,x,u);
xlabel('t','FontSize',20);
ylabel('x','FontSize',20);
title({'Implicit Method by Matrix'; ...
    [' k = ',num2str(k)];...
    [' h = ',num2str(h)];...
    [' r = ',num2str(r)]})

Crank-Nicolson Method

$\frac{u_j^{n+1}-u_j^{n}}{k} =\alpha\frac12 ( \frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2} + \frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2} )$ $(2+2r)u_j^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1} = (2-2r)u_j^{n}+ru_{j-1}^{n}+ru_{j+1}^{n} \quad r=\frac{\alpha k}{h^2}$