Andy123t
0%

矩形网的差分格式

发表于 更新于 分类于

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

二维 Possion 方程

考虑二维 Possion 方程:

其中 $\partial \Omega$ 为区域 $\Omega$ 的边界，$f(x,y)$ 和 $\phi(x)$ 为已知函数，$u|_{\partial \Omega}=\phi(x,y)$ 为边界条件。

五点差分格式

考虑 $\Omega$ 为矩形的情况，即 $a<x<b$，$c<x<d$. 取定沿 $x$ 轴和 $y$ 轴的步长 $h_1$ 和 $h_2$，$h_{1}=(b-a) / N, h_{2}=(d-c) / M$. 则

$(x_{i}, y_{j})$ 称为网格节点。

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

先定义向量 $\boldsymbol{u}_{j}$：

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

其中矩阵 $\boldsymbol{C}$、$\boldsymbol{D}$、向量 $\boldsymbol{f}_j$ 的定义如下，注意向量 $\boldsymbol{f}_j$ 的首尾元素已包含了 $x=a$ 和 $x=b$ 处的边界条件。

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

数值例子 1

其中 $f(x,y)=-2 \pi^{2} e^{\pi(x+y)}(\sin \pi x \cos \pi y+\cos \pi x \sin \pi y).$

方程的真解：

计算数值解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
% 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.png

计算收敛阶

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
% 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

输出结果

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

fdm2d1_error.png

数值例子 2

方程的真解：$u=(9+\pi^{2})^{-1} \cos(3 x) \sin(\pi y).$

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

差分格式:

边界条件为

计算数值解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
% 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)

输出结果

对应书上结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
  ......

计算收敛阶

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
% 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

输出结果

1
2
3
4
5
Error =
    0.1173    0.0473    0.0191    0.0085    0.0040    0.0019

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

fdm2d2_error.png

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