0%

极坐标形式的差分格式

极坐标形式的 Poisson 方程差分方法的 MATLAB 编程实现

极坐标形式的 Poisson 方程

其中 , 平面半带形域 .

方程的系数于 处奇异, 因此只当 时有意义。在 需补充 为光滑的条件 (原点是可去奇点). 由 , , 则知

差分格式

对变量 分别取等步长 . 令

则半带形域的网格节点为 , 它们在 平面上的分布如图

fdm_polar_mesh.png

差分方程:

下面是关于点 的差分方程

这样就得到 个方程的方程组.

注意: 实际计算时, 取 为边界比较合适.

数值例子 1

单位圆上的 Poisson 方程边值问题:

方程的真解:

差分格式推导:

fem_polar_derive1.jpg

fem_polar_derive2.jpg

计算数值解

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
% fdm_polar0.m
% finite difference method for polar coordinate problem
% -{Delta}_{r,theta}u(r,theta)=f(r,theta), [0,1]x[0,2*pi]
% u(1,theta)=0, theta in [0,2*pi]
% exact solution: ue=(1-r^2)/4.
clear all; close all;
N=50;
M=100;
dr=1/(N+1/2);
dthe=2*pi/M;
r=(dr/2:dr:1)';
the=(dthe:dthe:2*pi)';
% generate coordinates on the grid
[The,R]=meshgrid(the,r);
% generate the matrix of RHS
f=ones(M,N);
% constructing the coefficient matrix
e=[2/dr^2+8/(dr^2*dthe^2);(2./(dthe^2*r(2:N).^2)+2/dr^2)];
e1=[-2/dr^2;-(r(2:N-1)+dr/2)./(dr^2*r(2:N-1))];
e2=-(r(2:N)-dr/2)./(dr^2*r(2:N));
C=diag(e)+diag(e1,1)+diag(e2,-1);
D=diag([-4/(dr^2*dthe^2);-1./(dthe^2*r(2:N).^2)]);
E=diag([-4/(dr^2*dthe^2);-1./(dthe^2*r(2:N).^2)]);
A=kron(eye(M),C)+kron(diag(ones(M-1,1),1)+diag(1,1-M),D)...
+kron(diag(ones(M-1,1),-1)+diag(1,M-1),E);
% solving the linear system
f=f';
uh=reshape(A\f(:),N,M);
un=[uh;zeros(1,M)];
ue=(1-R.^2)/4;
% error on the node of the mesh
Err=abs(un-ue);
% compute maximum error
MaxErr=max(max(abs(un-ue)))
% plot the figure
[X,Y] = pol2cart(The,R);
mesh(X,Y,un)
set(gca,'fontsize',12)
xlabel('x','fontsize', 16)
ylabel('y','fontsize',16)
zlabel('u','fontsize',16,'Rotation',0)

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

输出结果

1
2
MaxErr =
1.9429e-14

当我取不同的步长 , 误差总是这么小, 是解的性质太好还是其他方面有问题, 还值得思考.

fdm_polar.png

数值例子 2

单位圆上的 Poisson 方程边值问题:

方程的真解:

计算数值解

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
% fdm_polar.m
% finite difference method for polar coordinate problem
% -{Delta}_{r,theta}u(r,theta)=f(r,theta), [0,1]x[0,2*pi]
% u(1,theta)=0, theta in [0,2*pi]
% exact solution: ue=(1-r^2)*r*(sin(theta)+cos(theta))/4.
clear all; close all;
N=50;
M=100;
dr=1/(N+1/2);
% dr=0.02;
dthe=2*pi/M;
r=(dr/2:dr:1)';
the=(0:dthe:2*pi)';
% N=length(r)-1;
% M=length(the);
% generate coordinates on the grid
[The,R]=meshgrid(the,r);
% generate the matrix of RHS
The1=The(1:N,2:end); R1=R(1:N,2:end);
f=2*R1.*(sin(The1)+cos(The1));
f=f';
% constructing the coefficient matrix
e=[2/dr^2+8/(dr^2*dthe^2);(2./(dthe^2*r(2:N).^2)+2/dr^2)];
e1=[-2/dr^2;-(r(2:N-1)+dr/2)./(dr^2*r(2:N-1))];
e2=-(r(2:N)-dr/2)./(dr^2*r(2:N));
C=diag(e)+diag(e1,1)+diag(e2,-1);
D=diag([-4/(dr^2*dthe^2);-1./(dthe^2*r(2:N).^2)]);
E=diag([-4/(dr^2*dthe^2);-1./(dthe^2*r(2:N).^2)]);
A=kron(eye(M),C)+kron(diag(ones(M-1,1),1)+diag(1,1-M),D)...
+kron(diag(ones(M-1,1),-1)+diag(1,M-1),E);
% solving the linear system
f=f';
% u=zeros(M+1,N+1);
uh=reshape(A\f(:),N,M);
un=[uh;zeros(1,M)];
un=[un(:,end),un];
ue=(1-R.^2).*R.*(sin(The)+cos(The))/4;
% error on the node of mesh
Err=abs(un(1:N,2:end)-ue(1:N,2:end));
% compute maximum error
MaxErr=max(max(abs(un-ue)))
% plot the figure
[X,Y] = pol2cart(The,R);
mesh(X,Y,ue)
set(gca,'fontsize',12)
xlabel('x','fontsize', 16)
ylabel('y','fontsize',16)
zlabel('u','fontsize',16)
view(36,24)

% print -dpng -r600 fdm_polar.png
% print -depsc2 fdm_polar.eps
1
2
MaxErr =
9.3183e-06

fdm_polarx.png

参考书籍

  • 李荣华, 微分方程数值解法, 第四版.