0%

一维差分格式

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

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

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

两点边值问题(常系数)

考虑二阶常微分方程边值问题(常系数):

其中 $q, f$ 为 $[a,b]$ 上的连续函数,$q \geqslant 0$;$\alpha, \beta$ 为给定常数。这是最简单的椭圆方程第一边值问题。

将区间 $[a,b]$ 分成 $N$ 等分,分点为

$h=(b-a) / N$。于是我们得到区间 $I=[a,b]$ 的一个网格剖分。$x_i$ 称为网格的节点,$h$ 称为步长

差分方程:

其中 $q_{i}=q(x_{i}), f_{i}=f(x_{i})$。

以上差分方程对于 $i=1,2, \cdots, N-1$ 都成立,加上边值条件 $u_{0}=\alpha, u_{N}=\beta$,就得到关于 $u_i$ 的差分格式:

它的解 $u_i$ 是 $u(x)$ 在 $x=x_i$ 处的差分解。

先定义向量 $\boldsymbol{u}$:

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

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

数值例子

方程的真解:$u(x)=\sin (\pi x)$.

计算数值解

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
% 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.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
% 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

输出结果

1
2
3
4
5
Error =
0.0084 0.0003 0.0001 0.0000 0.0000

order =
1.9990 2.0001 2.0000 2.0000

fdm1d1_error.png

两点边值问题(变系数)

考虑二阶常微分方程边值问题(变系数):

假定 $p \in C^{1}[a, b]$,$p(x) \geqslant p_{\min }>0$,$r, q, f \in C[a, b]$,$\alpha, \beta$ 是给定的常数。

首先取 $N+1$ 个节点:

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

记 $h_{i}=x_{i}-x_{i-1}$,称 $h=\max\limits _{i} h_{i}$ 为最大网格步长。

取相邻节点 $x_{i-1},x_i$ 的中点 $x_{i-\frac{1}{2}}=\frac{1}{2}\left(x_{i-1}+x_{i}\right)$ $(i=1,2, \cdots, N)$ 称为半整数点,则由节点

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

差分方程:

当网格均匀,即 $h=h_i (i=1,2, \cdots, N)$ 时,差分格式:

数值例子

方程的真解:$u(x)=\sin (\pi x)$.

计算数值解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
% 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.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
% 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

输出结果

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

fdm1d2_error.png