二阶椭圆型方程的五点差分格式-MATLAB实现

---

前言

数学来源人类的社会生产活动。现代数学的产生和发展与力学、物理学、天文学等应用学科的发展相辅相成的：它们为数学提出问题，而数学在解决这些问题的过程中所获得的更广泛、更深刻的结果反过来推动这些学科的发展。 现实世界中绝大多数事物的内外联系是及其复杂的，其状态随着时间、地点、条件的不同而不同，我们只能通过对问题进行简化和作某些假定，从中找出其状态和状态的变化规律之间的关系，也即一个或一些函数与它们的导数之间的关系，这种关系的数学表达就是微分方程。其中偏微分方程是指从物理问题中导出的反映客观物理量在各个地点、时刻之间相互制约关系的一些偏微分方程。本文主要介绍一些有关于求解二阶椭圆型方程中的poisson方程五点差分格式及其数值算法的实现。

一、问题介绍

此处我们主要研究熟知poisson方程：$$Lu=-\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=f(x,y),(x,y)\in \Omega$$其中$△u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$为laplace算子
且边界条件第一类边界条件：$$u(x,y)=\alpha (x,y),\forall (x,y)\in \Gamma$$其中$\Gamma$为边界曲线。

二、五点差分格式

1.区域剖分

对矩形区域$\Omega$进行矩形网格剖分：
取沿x轴和y轴方向上的步长分别为$h_1$和$h_2$,对区域$\Omega$进行剖分，则有$x_i$=$a+i{h}_1$, $y_i=b+j{h}_2$,其中 i ,j=0,1,2,… ,同时记$M=（b-a）/h_1$，$N=(d-c)/h_2，h=max(h_1,h_2)$.

2.差商代替导数

利用Taylor展式可以得到在内点$(x_i,y_j)$处的差分方程为$$L_h{u}_{ij}=-\frac{u_{i+1,j}-2u_{ij}+u_{i-1,j}}{h_1^2}-\frac{u_{i,j+1}-2u_{ij}+u_{i,j-1}}{h_2^2}=f_{ij}$$
$$i=1,...,M-1,j=1,...,N-1$$
特别的取$h_1=h_2=h$则上述差分格式可以简化为：$$4u_{ij}-(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})=h^2{f}_{ij}$$
如果$f\equiv 0(laplace方程)$，则有$$u_{ij}=\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})$$

3.解的存在唯一性

对于取$h_1=h_2=h$的情况，我们按自然顺序，即从下到上，从左到右的次序排列各内点，定义向量$uh=[u_1,u_2,...,u_{(M-1)(N-1)}]$,再结合边界条件，可以将上述各内点$(x_i,y_j)$差分格式改写成如下形势：$$H{u}_h=K$$
其中$$H=\left[\begin{array}{ccccc}B& -I& & & \\ -I& B& -I& & \\ & \ddots & \ddots & \ddots & \\ & & -I& B& -I\\ & & & -I& B\end{array}\right]$$为$(M-1)(N-1)$阶矩阵。$B为(M-1)阶矩阵$，$I为M-1阶单位矩阵$
由于H为块三对角矩阵且是非奇异的，所以上述线性方程组有唯一解。可以采用追赶法等数值算法求解。

三、实例及代码实现

$$f(x,y)=-4,\Omega =(0,1)\times (0,2)$$
$$\{\begin{array}{cc}& u(x,0)=x^2,u(x,2)=(x-2{)}^2,0<x<1\\ & u(0,y)=y^2,u(1,y)=(y-1{)}^2,0<y<2\end{array}$$
利用MATLAB可以求出精确解为$u(x,y)=(x-y{)}^2$
下面介绍求数值解的具体代码：
首先需要对区域$\Omega$进行剖分，取$h_1=h_2=h$,相应的代码实现如下：

function [yc,yd,xa,xb,x,y]=subdivision2(a,b,c,d,h)
%yc,yd,xa,xb分别为对应边界点的取值
%h为x轴和y轴方向的步长，即取h1=h2=h
%对应的区域为[a,b]×[c,d]
M=(b-a)/h;%剖分
N=(d-c)/h;
yc=zeros(M-1,1);%初始量
yd=zeros(M-1,1);
xa=zeros(N-1,1);
xb=zeros(N-1,1);
x=zeros(M-1,1);
y=zeros(N-1,1);
for i=1:M-1%上下边界
    x(i)=a+i*h;
    yc(i)=x(i)^2;%y=c边界上的值
    yd(i)=(x(i)-2)^2;%y=d
end
for i=1:N-1%左右边界
    y(i)=c+i*h;
    xa(i)=y(i)^2;%x=a
    xb(i)=(y(i)-1)^2;%x=b
end

其次需要生成系数矩阵$H$,由于其是块三对角矩阵，考虑使用cell数组产生块三对角矩阵，然后再利用cell2mat函数转换成所需要的矩阵。
用于生成系数矩阵$H$代码如下：

function [H] = Matrix2(M,N)%(M-1)*(N-1)为方程组的阶数
H=cell(N-1,N-1);
for i=1:N-1
    for j=1:N-1
        if i==j%主对角线
            H(i,j)={diag(ones(1,M-1)*4)-diag(ones(1,M-2),1)-diag(ones(1,M-2),-1)};
        else
if i==j-1||i==j+1%次对角线
   H(i,j)={-eye(M-1)};
else%其他地方
H(i,j)={zeros(M-1,M-1)};
end
        end
    end
end
H=cell2mat(H);%cell类型转换为mat
end

给出步长$h=h_1=h_2=0.25$时的系数矩阵$H$如下：

最后对于常数项的生成以及方程的求解和误差计算的相应代码如下：

function[error,Uh,t]=ellipticequations1(a,b,c,d,h)%矩形区域
%yc,yd,xa,xb分别为对应边界点的取值
%h为x轴和y轴方向的步长，即取h1=h2=h
%对应的区域为[a,b]×[c,d]
%error为误差，Uh为数值解，t为精确解
M=(b-a)/h;
N=(d-c)/h;
f=-4;%函数f(x,y)的值
[yc,yd,xa,xb,x,y]=subdivision2(a,b,c,d,h);%区域剖分，并求边界值
t=zeros((M-1)*(N-1),1);%精确值
K=(h^2)*f*ones((M-1)*(N-1),1);%处理边界
if M>=3&&N>=3 %特殊情形的处理
K(1)=yc(1)+xa(1)+f*(h^2);
K(M-1)=yc(M-1)+xb(1)+f*(h^2);
K((N-1)*(M-1))=xb(N-1)+yd(M-1)+f*(h^2);
K(1+(M-1)*(N-2))=xa(N-1)+yd(1)+f*(h^2);
else
    K(1)=K(1)+xa(1)+yc(1)+xb(1);
    K((M-1)*(N-1))=K((M-1)*(N-1))+xa(N-1)+xb(N-1)+yd(1);
end
A=zeros(M-1,N-1);
uh=zeros(N+1,M+1);%数值解矩阵
[H] = Matrix2(M,N);%产生方程组的系数矩阵
%特殊情形处理
if M>=4
for i=2:M-2
    K(i)=K(i)+yc(i);
end
for i=2+(N-2)*(M-1):(N-1)*(M-1)-1
    K(i)=K(i)+yd(i-(N-2)*(M-1));
end
end
if N>=4
for j=1:N-3
    K(M+(j-1)*(M-1))=K(M+(j-1)*(M-1))+xa(j+1);
    K(2*M-2+(j-1)*(M-1))=K(2*M-2+(j-1)*(M-1))+xb(j+1);
end
end
Uh=H\K;%求解
for j=1:N-1%精确值的产生
    for i=1:M-1
        A(i,j)=(x(i)-y(j))^2;
    end
end
for i=1:N-1
    t(1+(i-1)*(M-1):i*(M-1))=A(:,i);
end
%error=norm(Uh-t);%误差
error=max(abs(Uh-t));%误差

对于$h=h_1=h_2=0.25$的情况下，给出相应的误差如下：

error=1.06581410364015e-14;

补充：还可以在求解代码后面加上如下代码，从而作出数值解和精确解的图像：

%精确解画图
x0=a:h:b;
y0=c:h:d;
[X,Y]=meshgrid(x0,y0);%生成xoy网格点
U=(X-Y).^2;
mesh(X,Y,U);
colormap cool
xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis');
%legend('-','ln(x^2+y^2)');
hold on;
%数值解画图
%构造数值解矩阵
for j=1:M+1
    uh(1,j)=x0(j)^2;%yc
    uh(N+1,j)=(x0(j)-2)^2;%yd
end
for i=2:N
    uh(i,1)=y0(i)^2;%xa
    uh(i,M+1)=(y0(i)-1)^2;%xb
end
for i=2:N
    for j=2:M
        uh(i,j)=Uh((i-2)*(M-1)+j-1);
    end
end
plot3(X,Y,uh,'+');
legend('+',"数值解");
hold off;
end

对于取步长$h=h_1=h_2=0.025$的情况，作出数值解和精确解的图像如下：

注：以上代码是针对本例题进行编写的，如果需要用其求解其他类似问题，需要对其中的边界函数和精确解函数进行相应的修改，才可以使用。

由于本人水平有限，如有错误之处，敬请批评指正！