矩形网的差分格式

二维 Possion 方程五点差分格式的 MATLAB 编程实现。

二维 Possion 方程

考虑二维 Possion 方程:

$$\left\{ \begin{aligned} &-\Delta u=f(x,y),\quad (x,y)\in \Omega, \\ &u|_{\partial \Omega}=\phi(x,y),\quad (x, y) \in \partial \Omega , \end{aligned}\right.$$

其中 为区域 的边界， 和 为已知函数， 为边界条件。

五点差分格式

考虑 为矩形的情况，即 ，. 取定沿 轴和 轴的步长 和 ，. 则

$$x_{i}=a+i h_{1}, \quad 0 \leqslant i \leqslant N,$$ $$y_{j}=c+j h_{2}, \quad 0 \leqslant j \leqslant M.$$

称为网格节点。

二维 Possion 方程的五点差分格式：

$$-\frac{1}{h_{2}^{2}} u_{i, j-1}-\frac{1}{h_{1}^{2}} u_{i-1,j} +2\left(\frac{1}{h_{1}^{2}}+\frac{1}{h_{2}^{2}}\right) u_{i, j} -\frac{1}{h_{1}^{2}} u_{i+1, j}-\frac{1}{h_{2}^{2}} u_{i, j+1}=f\left(x_{i}, y_{j}\right),$$ $$1 \leqslant i \leqslant N-1,\quad 1 \leqslant j \leqslant M-1.$$

先定义向量 ：

$$\boldsymbol{u}_{j}=\left(u_{1 j}, u_{2 j}, \ldots, u_{N-1, j}\right)^{\mathrm{T}}, \quad 0 \leqslant j \leqslant M.$$

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

$$\boldsymbol{D} \boldsymbol{u}_{j-1}+\boldsymbol{C u}_{j}+\boldsymbol{D} \boldsymbol{u}_{j+1}=\boldsymbol{f}_{j}, \quad 1 \leqslant j \leqslant M-1.$$

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

$$\boldsymbol{C} = \left( \begin{array}{ccccc} \displaystyle 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right) & \displaystyle -\frac{1}{h_1^2} & & & \\ \displaystyle -\frac{1}{h_1^2} & \displaystyle 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right) & \displaystyle -\frac{1}{h_1^2} & & \\ & & & & \\ & \ddots & \ddots & \ddots & \\ & & & & \\ & &\displaystyle -\frac{1}{h_1^2} &\displaystyle 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right) & \displaystyle-\frac{1}{h_1^2} \\ & & & \displaystyle -\frac{1}{h_1^2} & \displaystyle 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right) \\ \end{array} \right),$$ $$\boldsymbol{D} = \left( \begin{array}{ccccc} \displaystyle-\frac{1}{h_2^2} & & & & \\ & \displaystyle-\frac{1}{h_2^2}& & & \\ & & \ddots & & \\ & & & \displaystyle-\frac{1}{h_2^2}& \\ & & & & \displaystyle-\frac{1}{h_2^2} \\ \end{array} \right),\quad \boldsymbol{f}_{j} = \left( \begin{array}{c} \displaystyle f(x_{1},y_{j})+\frac{1}{h_1^2}\phi(x_{0},y_{j}) \\ \displaystyle f(x_{2},y_{j}) \\ \vdots \\ \displaystyle f(x_{N-2},y_{j}) \\ \displaystyle f(x_{N-1},y_{j})+\frac{1}{h_1^2}\phi(x_{N},y_{j}) \\ \end{array} \right).$$

以上矩阵形式的差分格式还可以进一步写为如下的矩阵形式，注意等号右边向量的首尾元素加入了 和 处的边界条件。

$$\left(\begin{array}{ccccc} \boldsymbol{C} & \boldsymbol{D} & & & \\ \boldsymbol{D} & \boldsymbol{C}& \boldsymbol{D} & & \\ & \ddots & \ddots & \ddots & \\ & & \boldsymbol{D} & \boldsymbol{C}& \boldsymbol{D} \\ & & & \boldsymbol{D} & \boldsymbol{C} \\ \end{array} \right) \left(\begin{array}{c} \boldsymbol{u}_{1} \\ \boldsymbol{u}_{2} \\ \vdots \\ \boldsymbol{u}_{M-2} \\ \boldsymbol{u}_{M-1} \\ \end{array} \right) =\left(\begin{array}{c} \boldsymbol{f}_{1}-D\boldsymbol{u}_{0} \\ \boldsymbol{f}_{2} \\ \vdots \\ \boldsymbol{f}_{M-2} \\ \boldsymbol{f}_{M-1}-D\boldsymbol{u}_{M} \\ \end{array} \right).$$

数值例子 1

$$\left\{ \begin{aligned} &-\Delta u=f(x,y),\quad (x,y)\in \Omega=(0,1)\times (0,1), \\ &u=0,\quad (x, y) \in \partial \Omega . \end{aligned}\right.$$

其中

方程的真解：

$$u(x, y)=e^{\pi(x+y)} \sin (\pi x) \sin(\pi y), \quad (x,y)\in \Omega=(0,1)\times (0,1).$$

计算数值解

% fdm2d1.m
% finite difference method for 2D problem
% -d^2u/dx^2-d^2u/dy^2=f(x,y)
% f(x,y)=-2*pi^2*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y))
% exact solution: ue=exp(pi*x+pi*y)*sin(pi*x)*sin(pi*y)
clear all; close all;
% generate coordinates on the grid
h=0.02;
x=[0:h:1]';
y=[0:h:1]';
N=length(x)-1;
M=length(y)-1;
[X,Y]=meshgrid(x,y);
X=X(2:M,2:N);
Y=Y(2:M,2:N);
% generate the matrix of RHS
f=-2*pi^2*exp(pi*X+pi*Y).*(sin(pi*X).*cos(pi*Y)+cos(pi*X).*sin(pi*Y));
% constructing the coefficient matrix
C=4/h^2*eye(N-1)-1/h^2*diag(ones(N-2,1),1)-1/h^2*diag(ones(N-2,1),-1);
D=-1/h^2*eye(N-1);
A=kron(eye(M-1),C)+kron(diag(ones(M-2,1),1)+diag(ones(M-2,1),-1),D);
% solving the linear system
f=f';
u=zeros(M+1,N+1);
u(2:M,2:N)=reshape(A\f(:),N-1,M-1)';
u(:,1)=0;
u(:,end)=0;
ue=zeros(M+1,N+1);
ue(2:M,2:N)=exp(pi*X+pi*Y).*sin(pi*X).*sin(pi*Y);
% compute maximum error
Error=max(max(abs(u-ue)))
mesh(x,y,u)
%title('Finite Difference Method','fontsize',14)
set(gca,'fontsize',12)
xlabel('x','fontsize', 16)
ylabel('y','fontsize',16)
zlabel('u','fontsize',16)

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

计算收敛阶

% fdm2d1_error.m 
% finite difference method for 2D problem
% -d^2u/dx^2-d^2u/dy^2=f(x,y)
% f(x,y)=-2*pi^2*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y))
% exact solution: ue=exp(pi*x+pi*y)*sin(pi*x)*sin(pi*y)
clear all; close all;
Nvec=2.^[4:10];
Error=[];
for n=Nvec
   h=1/n;
   x=[0:h:1]';
   y=[0:h:1]';
   N=length(x)-1;
   M=length(y)-1;
   [X,Y]=meshgrid(x,y);
   X=X(2:M,2:N);
   Y=Y(2:M,2:N);
   % generate the matrix of RHS
   f=-2*pi^2*exp(pi*X+pi*Y).*(sin(pi*X).*cos(pi*Y)+cos(pi*X).*sin(pi*Y));
   % constructing the coefficient matrix
   e=ones(N-1,1);
   C=1/h^2*spdiags([-e 4*e -e],[-1 0 1],N-1,N-1);
   D=-1/h^2*eye(N-1);
   e=ones(M-1,1);
   A=kron(eye(M-1),C)+kron(spdiags([e e],[-1 1],M-1,M-1),D);
   % solving the linear system
   f=f';
   u=zeros(M+1,N+1);
   u(2:M,2:N)=reshape(A\f(:),N-1,M-1)';
   u(:,1)=0;
   u(:,end)=0;
   ue=zeros(M+1,N+1);
   % numerical solution
   ue(2:M,2:N)=exp(pi*X+pi*Y).*sin(pi*X).*sin(pi*Y);
   error=max(max(abs(u-ue)));     % maximum error
   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.46 0.50];
ay = [0.41 0.46];
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 fdm2d1_error.png
% print -depsc2 fdm2d1_error.eps

输出结果

Error =
    0.4587    0.1145    0.0286    0.0072    0.0018    0.0004    0.0001

order =
    2.0026    2.0002    1.9994    2.0000    2.0000    2.0000

数值例子 2

$$\left\{ \begin{aligned} &-\Delta u=\cos 3 x \sin \pi y, && (x, y) \in G=(0, \pi) \times(0,1), \\ &u(x, 0)=u(x, 1)=0, && 0 \leqslant x \leqslant 1, \\ &u_{x}(0, y)=u_{x}(\pi, y)=0, && 0 \leqslant y \leqslant 1. \end{aligned}\right.$$

方程的真解：

对应《微分方程数值解法》104-105 页数值例子。

差分格式:

$$-\left(\frac{u_{i+1, j}-2 u_{i j}+u_{i-1, j}}{h_{1}^{2}}+\frac{u_{i, j+1}-2 u_{i j}+u_{i, j-1}}{h_{2} 2}\right)=\cos(3 x_{i}) \sin(\pi y_{j}),$$ $$i, j=1,2, \cdots, N-1.$$

边界条件为

$$\begin{aligned} &u_{i 0}=u_{i N}=0,\quad i=0, \cdots, N, \\ &u_{0 j}=u_{1 j},\quad j=1, \cdots, N-1, \\ &u_{N j}=u_{N-1, j},\quad j=1, \cdots, N-1. \end{aligned}$$

计算数值解

% fdm2d2.m
% finite difference method for 2D problem
% -\Delta u=cos(3*x)*sin(pi*y)  in (0,pi)x(0,1)
% u(x,0)=u(x,1)=0  in [0,pi]
% u_x(0,y)=u_x(pi,y)=0 in [0,1]
% exact solution: ue=(9+pi^2)^(-1)*cos(3*x)*sin(pi*y)
clear all; close all;
N=4    % N=4 8 16 32
% coordinates on the grid
h1=pi/N;  h2=1/N;
x=[0:h1:pi]';
y=[0:h2:1]';
[X,Y]=meshgrid(x,y);
X1=X(2:N,2:N);
Y1=Y(2:N,2:N);
% generate the matrix of RHS
f=cos(3*X1).*sin(pi*Y1);
% constructing the coefficient matrix
e=ones(N-1,1);
C=diag([1/h1^2+2/h2^2, (2/h1^2+2/h2^2)*ones(1,N-3), 1/h1^2+2/h2^2])...
    -1/h1^2*diag(ones(N-2,1),1)-1/h1^2*diag(ones(N-2,1),-1);
D=-1/h2^2*eye(N-1);
% A=kron(eye(N-1),C)+kron(diag(ones(N-2,1),1)+diag(ones(N-2,1),-1),D);
A=kron(eye(N-1),C)+kron(spdiags([e e],[-1 1],N-1,N-1),D);
% solving the linear system
f=f';
u=zeros(N+1,N+1);
u(2:N,2:N)=reshape(A\f(:),N-1,N-1)';
% Neumann boundary condition
u(:,1)=u(:,2);
u(:,end)=u(:,end-1);
ue=1/(9+pi^2)*(cos(3*X)).*(sin(pi*Y));

format long
% value of u and ue in the selected points (i*pi/4,j/4), i,j=1,2,3.
u_select=u(N/4+1:N/4:3*N/4+1,N/4+1:N/4:3*N/4+1)
ue_select=u(N/4+1:N/4:3*N/4+1,N/4+1:N/4:3*N/4+1)

输出结果

对应书上结果

N =
     4

u_select =
  -0.045480502664596  -0.000000000000000   0.045480502664596
  -0.064319143691817  -0.000000000000000   0.064319143691817
  -0.045480502664596  -0.000000000000000   0.045480502664596

ue_select =
  -0.045480502664596  -0.000000000000000   0.045480502664596
  -0.064319143691817  -0.000000000000000   0.064319143691817
  -0.045480502664596  -0.000000000000000   0.045480502664596

% N = 8, 16, 32
  ......

计算收敛阶

% fdm2d2_error.m 
% finite difference method for 2D problem
% -\Delta u=cos(3*x)*sin(pi*y)  in (0,pi)x(0,1)
% u(x,0)=u(x,1)=0  in [0,pi]
% u_x(0,y)=u_x(pi,y)=0 in [0,1]
% exact solution: ue=(9+pi^2)^(-1)*cos(3*x)*sin(pi*y)
clear all; close all;
Nvec=2.^[2:7];
Error=[];
for N=Nvec
    % coordinates on the grid
    h1=pi/N;  h2=1/N;
    x=[0:h1:pi]'; 
    y=[0:h2:1]';
    [X,Y]=meshgrid(x,y);
    X1=X(2:N,2:N);
    Y1=Y(2:N,2:N);
    % generate the matrix of RHS
    f=cos(3*X1).*sin(pi*Y1);
    % constructing the coefficient matrix
    e=ones(N-1,1);
    C=diag([1/h1^2+2/h2^2, (2/h1^2+2/h2^2)*ones(1,N-3), 1/h1^2+2/h2^2])...
        -1/h1^2*diag(ones(N-2,1),1)-1/h1^2*diag(ones(N-2,1),-1);
    D=-1/h2^2*eye(N-1);
    %A=kron(eye(N-1),C)+kron(diag(ones(N-2,1),1)+diag(ones(N-2,1),-1),D);
    A=kron(eye(N-1),C)+kron(spdiags([e e],[-1 1],N-1,N-1),D);
    % solving the linear system
    f=f'; 
    u=zeros(N+1,N+1);
    u(2:N,2:N)=reshape(A\f(:),N-1,N-1)';
    % Neumann boundary condition
    u(:,1)=u(:,2); 
    u(:,end)=u(:,end-1);
    ue=1/(9+pi^2)*(cos(3*X)).*(sin(pi*Y));
    error=max(max(abs(u-ue)));   % maximum error
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec), log10(Nvec.^(-1)), '--')
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.64 0.60];
ay = [0.69 0.64];
annotation('textarrow',ax,ay,'String','slope = -1 ','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 fdm2d2_error.png
% print -depsc2 fdm2d2_error.eps

输出结果

Error =
    0.1173    0.0473    0.0191    0.0085    0.0040    0.0019

order =
   -1.3104   -1.3070   -1.1767   -1.0911   -1.0458

思考：五点差分格式的收敛阶应该 2 阶，为什么结果是 1 阶，思考之后发现是因为边界条件。