极坐标形式的差分格式

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

极坐标形式的 Poisson 方程

$$-\Delta_{r, \theta} u=-\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial u}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} u}{\partial \theta^{2}}\right]=f(r, \theta),$$

其中 $r=\sqrt{x^{2}+y^{2}},~\tan \theta =y / x$, $r \theta$ 平面半带形域 $\{0 \leqslant r < \infty, 0 \leqslant \theta \leqslant 2\pi \}$.

方程的系数于 $r=0$ 处奇异, 因此只当 $r>0$ 时有意义。在 $r=0$ 需补充 $u$ 为光滑的条件 (原点是可去奇点). 由 $u(0, \theta)=u(r, \theta)-r u_{r}(r, \theta)+O(r^{2})$, $r \rightarrow 0$, 则知

$$\lim _{r \rightarrow 0^{+}} r \frac{\partial u}{\partial r}=0.$$

差分格式

对变量 $r, \theta$ 分别取等步长 $h_r$ 和 $h_{\theta}$. 令

$$r_{i}=(i+0.5) h_{r},\ i=0,1,2, \cdots,N-1,$$ $$\theta_{j}=(j+1) h_{\theta},\ j=0,1, \cdots, M-1,\ h_{\theta}=2 \pi / M.$$

则半带形域的网格节点为 $(r_i,\theta_j)$, 它们在 $r\theta$ 平面上的分布如图

差分方程：

$$\begin{aligned} &-\frac{1}{r_{i}} \frac{r_{i+\frac{1}{2}} u_{i+1, j}-\left(r_{i+\frac{1}{2}}+r_{i-\frac{1}{2}}\right) u_{i j}+r_{i-\frac{1}{2}} u_{i-1, j}}{h_{r}^{2}} \\ &-\frac{1}{r_{i}^{2}} \frac{u_{i, j+1}-2 u_{i j}+u_{i, j-1}}{h_{\theta}^{2}}=f\left(r_{i}, \theta_{j}\right),\quad 1 \leqslant i \leqslant N-1. \end{aligned}$$

下面是关于点 $(r_0, \theta_{j})$ 的差分方程

$$-\frac{2}{h_{r}} \frac{u_{1 j}-u_{0 j}}{h_{r}}-\frac{4}{h_{r}^{2}} \frac{u_{0, j+1}-2 u_{0 j}+u_{0, j-1}}{h_{\theta}^{2}}=f_{0 j}.$$

这样就得到 $N~(i=0,1,\cdots,N-1)$ 个方程的方程组.

注意： 实际计算时, 取 $(N+h_r/2)$ 为边界比较合适.

数值例子 1

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

$$\left\{\begin{aligned} &-\Delta u=1, \quad \Omega=\left\{(x,y) \mid x^{2}+y^{2}<1\right\}, \\ &\left.u\right|_{\partial \Omega}=0. \end{aligned}\right.$$

方程的真解：

$$u(x, y)=\frac{1-x^{2}-y^{2}}{4}.$$

差分格式推导：

计算数值解

% 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

输出结果

MaxErr =
   1.9429e-14

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

数值例子 2

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

$$\left\{\begin{aligned} &-\Delta u=1, \quad \Omega=\left\{(x,y) \mid x^{2}+y^{2}<1\right\}, \\ &\left.u\right|_{\partial \Omega}=0. \end{aligned}\right.$$

方程的真解：

$$u(x, y)=\frac{(1-x^{2}-y^{2})(x+y)}{4}.$$

计算数值解

% 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

MaxErr =
    9.3183e-06

参考书籍

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