用有限差分法求解矩形域上的Poisson方程的Matlab实现

物理过程,都可用椭圆型方程来描述。其中最典型的方程是泊松(Poisson)方程。 本文件给出了使用有限差分法求解泊松方程的两个算例及MATLAB源程序。

function poisson5spot（h）
if nargin<1 %默认步长h=0.02
    h=0.02;
end
x=0:h:3/2;
y=0:h:1/2;
nx=length（x）-1;%取区域的内点
ny=length（y）-1;
%===========构造矩阵B==========
B=eye（nxnx）*4;
for i=1:nx-1
    B（ii+1）=-1;
end
for i=2:nx
    B（ii-1）=-1;
end
I=-eye（nxnx）;
%===========构造系数矩阵A========
A=zeros（nx*nynx*ny）;
A（1:nx1:nx）=B;
%由区域的对称性，正方形网格最下面一行的差分形式
%改为4u（ij）-u（i-1j）-u（i+1j）-2*u（ij+1）
%因为u（ij+1）=u（ij-1）
A（1:nxnx+1:2*nx）=2*I;  
A（nx+1:2*nxnx+1:2*nx）=B;
%矩形网格左下角第一个点的差分形式改为：
%4u（ij）-u（i-1j）-u（i+1j）-u（ij+1）-u（ij-1）
%=4u（ij）-2u（ij+1）-2u（i+1j）
A（12）=-2;
%为了方便，本程序往梯形中增加了一个三角形，以方便编程
A（nx+1:2*nx1:nx）=I;
for i=2:ny-1
    A（i*nx+1:（i+1）*nxi*nx+1:（i+1）*nx）=B;
    A（（i-1）*nx+1:i*nxi*nx+1:（i+1）*nx）=I;
    A（i*nx+1:（i+1）*nx（i-1）*nx+1:i*nx）=I;
end
%==========================
%bi=h*h*f;右端向量
b=h*h*ones（nx*ny1）;
%由于对称性，左侧三角形中有u（ij）-u（ji）=0     i>j
%因此令A（ij）=1 A（ji）=-1
%所以在本程序中多计算了（ny^2-my）/2 个点的函数值
%但对程序的影响并不是很大
for i=2:ny
    A（（i-1）*nx+1:（i-1）*nx+i-1:）=0;
    for j=1:i-1
        A（（i-1）*nx+j（i-1）*nx+j）=1;
        A（（i-1）*nx+j（j-1）*nx+i）=-1;
        b（（i-1）*nx+j）=0;
    end
end
x=A\b;%为了方便采用左除求解网格点上的函数值
%x=gmres（Ab100）;
%按顺序将x赋值给u
u=zeros（ny+1nx+1）;
for i=1:ny
    for j=1:nx
        u（ij）=x（（i-1）*nx+j）;
    end
end
%根据对称性，给网格最上面一行的点赋值
u（ny+11:ny）=u（1:nyny+1）;
 
%=============作Poisson方程在区域上的图形===========
[xy]=meshgrid（0:h:3/20:h:1/2）;
hold on
surf（xyu）%11第一象限的第一块区域，下面的以此类推
surf（yxu）%12
surf（-yxu）%21
surf（-xyu）;%22
surf（-x-yu）%31
surf（-y-xu）%32
surf（y-xu）%41
surf（x-yu）;%42

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件     1796928  2019-03-19 13:34  泊松方程的有限差分法的MATLAB实现.pdf
     文件        3945  2019-03-19 13:15  用有限差分法求解矩形域上的Poisson方程.m
     文件      600064  2019-03-19 13:42  Possion方程的差分法求解.doc
     文件        1807  2019-03-19 13:40  可分割为矩形的复杂区域的Possion 方程的差分方法求解.m