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