一维差分格式

在偏微分方程的数值解法中，有限差分法数学概念直观，推导自然，是发展较早且比较成熟的数值方法。由于计算机只能存储有限个数据和做有限次运算，所以任何一种用计算机解题的方法，都必须把连续问题（微分方程的边值问題、初值问题等）离散化，最终化成有限形式的线性代数方程组。

学习有限差分方法，差分解的存在唯一性、收敛性及稳定性等理论知识请参考书籍《微分方程数值解法》（李荣华）。

这里重点介绍问题差分格式的 MATLAB 编程实现。

两点边值问题（常系数）

考虑二阶常微分方程边值问题（常系数）：

$$\begin{aligned} & L u=-\frac{\mathrm{d}^{2} u}{\mathrm{d} x^{2}}+\frac{\mathrm{d} u}{\mathrm{d} x}+q u=f, \quad a < x < b, \\ &u(a)=\alpha, \quad u(b)=\beta. \end{aligned}$$

其中 为 上的连续函数，； 为给定常数。这是最简单的椭圆方程第一边值问题。

将区间 分成 等分，分点为

$$x_{i}=a+i h \quad i=0,1, \cdots, N,$$

。于是我们得到区间 的一个网格剖分。 称为网格的节点， 称为步长。

差分方程：

$$L_{h} u_{i}=-\frac{u_{i+1}-2 u_{i}+u_{i-1}}{h^{2}}+\frac{u_{i+1}-u_{i-1}}{2h}+q_{i} u_{i}=f_{i},\quad 1 \leqslant j \leqslant N-1.$$

其中 。

以上差分方程对于 都成立，加上边值条件 ，就得到关于 的差分格式：

$$\begin{aligned} & L_{h} u_{i}=-\frac{u_{i+1}-2 u_{i}+u_{i-1}}{h^{2}}+\frac{u_{i+1}-u_{i-1}}{2h}+q_{i} u_{i}=f_{i}, \quad i=1,2, \cdots, N-1, \\ & u_{0}=\alpha, \quad u_{N}=\beta. \end{aligned}$$

它的解 是 在 处的差分解。

先定义向量 ：

$$\boldsymbol{u}=(u_{1}, u_{2}, \cdots, u_{N-1})^{\mathrm{T}}.$$

差分格式可以写为矩阵形式：

$$\boldsymbol{A}\boldsymbol{u}=\boldsymbol{f}.$$

其中矩阵 、向量 的定义如下，注意向量 的首尾元素已包含了 和 处的边界条件。

$$\boldsymbol{A}=\left[\begin{array}{cccccc} \displaystyle \frac{2}{h^{2}}+q_{1} & \displaystyle\frac{1}{2h}-\frac{1}{h^{2}} & & & \\[8pt] \displaystyle-\frac{1}{2h}-\frac{1}{h^{2}} & \displaystyle \frac{2}{h^{2}}+q_{2} & \displaystyle\frac{1}{2h}-\frac{1}{h^{2}} & & \\[8pt] & & & & \\ & \ddots & \ddots & \ddots & \\[8pt] & & & & \\ & & \displaystyle-\frac{1}{2h}-\frac{1}{h^{2}} & \displaystyle\frac{2}{h^{2}}+q_{N-2}& \displaystyle\frac{1}{2h}-\frac{1}{h^{2}} \\[8pt] & & & \displaystyle-\frac{1}{2h}-\frac{1}{h^{2}} & \displaystyle \frac{2}{h^{2}}+q_{N-1} \end{array}\right],$$ $$\boldsymbol{f}=(f_1+\frac{\alpha}{h^2}+\frac{\alpha}{2h}, f_2, \cdots, f_{N-1}+\frac{\beta}{h^2}-\frac{\beta}{2h})^{\mathrm{T}}.$$

数值例子

$$\left\{\begin{aligned} &-\frac{\mathrm{d}^{2} u}{\mathrm{d} x^{2}}+\frac{\mathrm{d} u}{\mathrm{d} x}=\pi^{2} \sin (\pi x)+\pi \cos (\pi x),\quad 0 < x < 1, \\ &u(0)=0,\quad u(1)=0. \end{aligned} \right.$$

方程的真解：.

计算数值解

% fdm1d1.m
% finite difference method for 1D problem
% -u''+u'=pi^2*sin(pi*x)+pi*cos(pi*x) in [0,1]
% u(0)=0, u(1)=0 ;
% exact solution: u=sin(pi*x)
clear all; close all;
h=0.05;
x=0:h:1;
N=length(x)-1;
A=diag((2/h^2)*ones(N-1,1))...
    +diag((1/(2*h)-1/h^2)*ones(N-2,1),1)...
    +diag((-1/(2*h)-1/h^2)*ones(N-2,1),-1);
b=pi^2*sin(pi*x(2:N))+pi*cos(pi*x(2:N));
u=A\b';
u=[0;u;0];
ue=sin(pi*x)';
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
Error=max(abs(u-ue))
legend('Exact ','Numerical','location','NorthEast')
%title('Finite Difference Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize', 16), ylabel('u','fontsize',16)

% print -dpng -r600 fdm1d1.png
% print -depsc2 fdm1d1.eps

计算收敛阶

% fdm1d1_error.m
% finite difference method for 1D problem
% -u''+u'=pi^2*sin(pi*x)+pi*cos(pi*x)  in [0,1]
% u(0)=0, u(1)=0 ;
% exact solution : u=sin(pi*x)
clear all; close all;
Nvec=[10 50 100 500 1000];    % Number of partitions
Error=[];
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    x=0:h:1;
    N=length(x)-1;
    A=diag((2/h^2)*ones(N-1,1))...
        +diag((1/(2*h)-1/h^2)*ones(N-2,1),1)...
        +diag((-1/(2*h)-1/h^2)*ones(N-2,1),-1);
    b=pi^2*sin(pi*x(2:N))+pi*cos(pi*x(2:N));
    u=A\b';
    u=[0;u;0];
    ue=sin(pi*x');
    error=max(abs(u-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec), log10(Nvec.^(-2)), '--')
grid on
%title('Convergence of Finite Difference Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize', 14), ylabel('log_{10}Error','fontsize',14),

% add annotation of slope
ax = [0.55 0.52];
ay = [0.67 0.62];
annotation('textarrow',ax,ay,'String','slope = -2 ','fontsize',14)

% computating convergence order
for i=1:length(Nvec)-1
    order(i)=-log(Error(i)/Error(i+1))/(log(Nvec(i)/Nvec(i+1)));
end
Error   
order

% print -dpng -r600 fdm1d1_error.png
% print -depsc2 fdm1d1_error.eps

输出结果

Error =
    0.0084    0.0003    0.0001    0.0000    0.0000

order =
    1.9990    2.0001    2.0000    2.0000

两点边值问题（变系数）

考虑二阶常微分方程边值问题（变系数）：

$$\begin{aligned} & L u=-\frac{\mathrm{d}}{\mathrm{d} x}\left(p \frac{\mathrm{d} u}{\mathrm{d} x}\right)+r \frac{\mathrm{d} u}{\mathrm{d} x}+q u=f, \quad a < x < b, \\ &u(a)=\alpha, \quad u(b)=\beta. \end{aligned}$$

假定 ，，， 是给定的常数。

首先取 个节点：

$$a=x_{0}<x_{1}<\cdots<x_{i}<\cdots<x_{N}=b$$

将区间 $I=[a,b]$ 分成 $N$ 个小区间：

$$I_{i}: x_{i-1} \leqslant x \leqslant x_{i}, \quad i=1,2, \cdots, N.$$

记 ，称 为最大网格步长。

取相邻节点 的中点 称为半整数点，则由节点

$$a=x_{0}<x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{i-\frac{1}{2}}<\cdots<x_{N-\frac{1}{2}}<x_{N}=b.$$

又构成 $[a,b]$ 的一个剖分，称为对偶剖分。

差分方程：

$$\begin{aligned} L_{h} u_{i} \equiv & -\frac{2}{h_{i}+h_{i+1}}\left[p_{i+\frac{1}{2}} \frac{u_{i+1}-u_{i}}{h_{i+1}}-p_{i-\frac{1}{2}} \frac{u_{i}-u_{i-1}}{h_{i}}\right] +\\ & \frac{r_{i}}{h_{i}+h_{i+1}}\left(u_{i+1}-u_{i-1}\right)+q_{i} u_{i}=f_{i}, \quad i=1, \cdots, N-1 \\ & u_{0}=\alpha, \quad u_{N}=\beta. \end{aligned}$$

当网格均匀，即 时，差分格式：

$$\begin{aligned} L_{h} u_{i}=&-\frac{1}{h^{2}}\left[p_{i+\frac{1}{2}} u_{i+1}-\left(p_{i+\frac{1}{2}}+p_{i-\frac{1}{2}}\right) u_{i}+p_{i-\frac{1}{2}} u_{i-1}\right]+\\ & r_{i} \frac{u_{i+1}-u_{i-1}}{2 h}+q_{i} u_{i}=f_{i},\quad i=1, \cdots, N-1 \\ & u_{0}=\alpha, \quad u_{N}=\beta. \end{aligned}$$

数值例子

$$\left\{ \begin{aligned} &-\frac{\mathrm{d}}{\mathrm{d} x}\left(x \frac{\mathrm{d} u}{\mathrm{d} x}\right)+x \frac{\mathrm{d} u}{\mathrm{d} x}=\pi^{2}x \sin(\pi x)+\pi(x-1) \cos(\pi x),\quad 0 < x < 1, \\ &u(0)=0,\quad u(1)=0. \end{aligned}\right.$$

方程的真解：.

计算数值解

% fdm1d2.m
% finite difference method for 1D problem
% -(xu')'+x*u'=pi^2*x*sin(pi*x)-pi*cos(pi*x)+pi*x*cos(pi*x) in [0,1]
% u(0)=0, u(1)=0 ;
% exact solution : u=sin(pi*x)
clear all; close all;
h=0.05;
x=0:h:1;
N=length(x)-1;
A=diag(2*x(2:N)./h^2)+diag(x(2:N-1)./(2*h)-(x(2:N-1)+0.5*h)./h^2,1)...
    +diag(-x(3:N)./(2*h)-(x(3:N)-0.5*h)./h^2,-1);
b=pi^2*x(2:N).*sin(pi*x(2:N))+pi*(x(2:N)-1).*cos(pi*x(2:N));
u=A\b';
u=[0;u;0];
ue=sin(pi*x');
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
Error=max(abs(u-ue))
legend('Exact ','Numerical','location','NorthEast')
%title('Finite Difference Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize', 16), ylabel('u','fontsize',16)

% print -dpng -r600 fdm1d2.png
% print -depsc2 fdm1d2.eps

计算收敛阶

% fdm1d2_error.m
% finite difference method for 1D problem
% -(xu')'+x*u'=pi^2*x*sin(pi*x)-pi*cos(pi*x)+pi*x*cos(pi*x) in [0,1]
% u(0)=0, u(1)=0 ;
% exact solution : u=sin(pi*x)
clear all; close all;
Nvec=[10 20 50 100 200 500 1000];    % Number of partitions
Error=[];
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    x=0:h:1;
    N=length(x)-1;
    A=diag(2*x(2:N)./h^2)+diag(x(2:N-1)./(2*h)-(x(2:N-1)+0.5*h)./h^2,1)...
        +diag(-x(3:N)./(2*h)-(x(3:N)-0.5*h)./h^2,-1);
    b=pi^2*x(2:N).*sin(pi*x(2:N))+pi*(x(2:N)-1).*cos(pi*x(2:N));
    u=A\b';
    u=[0;u;0];
    ue=sin(pi*x');
    error=max(abs(u-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec), log10(Nvec.^(-2)), '--')
grid on
%title('Convergence of Finite Difference Method','fontsize',14)
set(gca,'fontsize',14)
xlabel('log_{10}N','fontsize', 14), ylabel('log_{10}Error','fontsize',14),

% add annotation of slope
ax = [0.57 0.53];
ay = [0.68 0.63];
annotation('textarrow',ax,ay,'String','slope = -2 ','fontsize',14)

% computating convergence order
for i=1:length(Nvec)-1
    order(i)=-log(Error(i)/Error(i+1))/(log(Nvec(i)/Nvec(i+1)));
end
Error
order

% print -dpng -r600 fdm1d2_error.png
% print -depsc2 fdm1d2_error.eps

输出结果

Error =
    0.0059    0.0013    0.0002    0.0001    0.0000    0.0000    0.0000

order =
    2.1277    2.0216    1.7479    1.8159    1.8623    1.8930