matlab程序解泊松方程.docx
《matlab程序解泊松方程.docx》由会员分享，可在线阅读，更多相关《matlab程序解泊松方程.docx（9页珍藏版）》请在冰豆网上搜索。
matlab程序解泊松方程
求解泊松方程的
functionFinite_element_tri（Imax）
%用有限元法求解三角形形区域上的Possion方程
Jmax=2*Imax;
%其中ImaxJmax分别表示x轴和y轴方向的网格数，其中Jmax等于Imax的两倍
%定义一些全局变量
globalndmnelna
%ndm总节点数
@
%nel基元数
%na表示活动节点数
V=0;J=0;X0=1/Imax;Y0=X0;%V=0为边界条件
domain_tri%调用函数画求解区域
[X,Y,NN,NE]=setelm_tri（Imax,Jmax）;%给节点和三角形元素编号，并设定节点坐标
%以下求解有限元方程的求系数矩阵
T=zeros（ndm,ndm）;
forn=1:
nel
…
n1=NE（1,n）;n2=NE（2,n）;n3=NE（3,n）;%整体编号
s=abs（（X（n2）-X（n1））*（Y（n3）-Y（n1））-（X（n3）-X（n1））*（Y（n2）-Y（n1）））/2;%三角形面积
fork=1:
3
ifn1<=na|n2<=na
T（n1,n2）=T（n1,n2）+（（Y（n2）-Y（n3））*（Y（n3）-Y（n1））+（X（n3）-X（n2））*（X（n1）-X（n3）））/（4*s）;
T（n2,n1）=T（n1,n2）;
T（n1,n1）=T（n1,n1）+（（Y（n2）-Y（n3））^2+（X（n3）-X（n2））^2）/（4*s）;%V=0则边界积分为零，非零时积分编程类似，再加边界积分。
end
|
k=n1;n1=n2;n2=n3;n3=k;%轮换坐标将值赋入3阶主子矩阵中
end
end
M=T（1:
na,1:
na）;
%求有限元方程的右端项
f=X;%场源函数
G=zeros（na,1）;
forn=1:
nel
—
n1=NE（1,n）;n2=NE（2,n）;n3=NE（3,n）;
s=abs（（X（n2）-X（n1））*（Y（n3）-Y（n1））-（X（n3）-X（n1））*（Y（n2）-Y（n1）））/2;
fork=1:
3
ifn1<=na
G（n1）=G（n1）+（2*f（n1）+f（n2）+f（n3））*s/12;%f在单元上为线性差值时场域单元的积分公式
end
n4=n1;n1=n2;n2=n3;n3=n4;%轮换坐标标
end
》
end
F=M\G;%求解方程得结果
NNV=zeros（Imax+1,Jmax+1）;
fi=zeros（ndm,1）;
fi（1:
na）=F（1:
na）;
fi（na+1:
ndm）=V;
forj=0:
Jmax
fori=0:
Imax
（
n=NN（i+1,j+1）;
ifn<=0
n=na+1;
end
NNV（i+1,j+1）=fi（n）;
end
end
figure
（2）
—
imagesc（NNV）;%画解函数的平面图
X1=zeros（1,Imax+1）;
Y1=zeros（1,Jmax+1）;
fori=1:
Imax+1
X1（i）=（i-1）*X0;
end
fori=1:
Jmax+1
Y1（i）=（i-1）*Y0;
:
end
figure（3）
surf（X1,Y1,NNV'）;%画解函数的曲面图
%以下是结果的输出
fid=fopen（'','w'）;
fprintf（fid,'\n*********有限元法求解三角形区域上Possion方程的结果**********\n\n'）;
L=[1:
ndm]';
fprintf（fid,'\n\n节点编号坐标分量x坐标分量yu（x,y）的值\n\n'）;
^
fori=1:
ndm
fprintf（fid,'%8d%%%\n',L（i）,X（i）,Y（i）,fi（i））;
end
fclose（fid）;
functiondomain_tri
%画求解区域
xy=[01;0-1;10];
-
A=zeros（3,3）;
A（1,1）=2;A（1,2）=-1;A（1,3）=-1;
A（2,2）=2;A（2,1）=-1;A（2,3）=-1;
A（3,3）=2;A（3,2）=-1;A（3,1）=-1;
A=sparse（A）;
figure
（1）;
gplot（A,xy）;
function[X,Y,NN,NE]=setelm_tri（Imax,Jmax）
^
%给节点和三角形单元编号，并设定节点坐标
%定义一些全局变量
globalndmnelna
%I1I2J1J2ImaxJmax分别描述网线纵向和横向数目的变量
%XY表示节点坐标
%NN描述节点编号
%NE描述各单元局部节点编号与总体编号对应的矩阵
%ndm总节点数
（
%nel单元数
%na表示不含边界的节点数
nlm=Imax*Jmax;
dx=1/Imax;
dy=1/Jmax;
X=nlm:
1;
Y=nlm:
1;
NN=zeros（Imax+1,Jmax+1）;
@
n1=0;
forj=3:
Jmax/2
fori=2:
j-1
n1=n1+1;
NN（i,j）=n1;%X=i列，Y=j行处节点
X（n1）=（i-1）*dx;
Y（n1）=-1+（j-1）*dy;
end
.
end
k=Jmax/2+1;
forj=Jmax/2+1:
Jmax-1%三角形区域上下两部分节点坐标分别求
k=k-1;
fori=2:
k
n1=n1+1;
NN（i,j）=n1;
X（n1）=（i-1）*dx;
[
Y（n1）=1+（j-Jmax-1）*dy;
end
end
na=n1;%不含边界节点数
forj=Jmax+1:
-1:
Jmax/2+1%降序
n1=n1+1;
NN（1,j）=n1;
X（n1）=0;
…
Y（n1）=1+（j-Jmax-1）*dy;
end
forj=Jmax/2:
-1:
1
n1=n1+1;
NN（1,j）=n1;
X（n1）=0;
Y（n1）=-1+（j-1）*dy;
end%
-
fori=2:
Imax+1
n1=n1+1;
NN（i,i）=n1;
X（n1）=（i-1）*dx;
Y（n1）=-1+（i-1）*dy;
end
K=0;
fori=Imax:
-1:
2
^
K=K+2;
n1=n1+1;
NN（i,i+K）=n1;
X（n1）=（i-1）*dx;
Y（n1）=1+（i+K-Jmax-1）*dy;
end
%以上四个循环为对边界节点进行编号
ndm=n1;
《
NE=zeros（3,2*ndm）;n1=0;
forj=3:
Jmax/2
fori=2:
j-1
n1=n1+1;
NE（1,n1）=NN（i,j）;
NE（2,n1）=NN（i-1,j+1）;
NE（3,n1）=NN（i-1,j）;
n1=n1+1;
（
NE（1,n1）=NN（i,j）;
NE（2,n1）=NN（i,j+1）;
NE（3,n1）=NN（i-1,j+1）;
end
end
k=Jmax/2+1;
forj=Jmax/2+1:
Jmax-1
k=k-1;
（
fori=2:
k
n1=n1+1;
NE（1,n1）=NN（i,j）;
NE（2,n1）=NN（i-1,j+1）;
NE（3,n1）=NN（i-1,j）;
n1=n1+1;
NE（1,n1）=NN（i,j）;
NE（2,n1）=NN（i,j+1）;
NE（3,n1）=NN（i-1,j+1）;
end
end%内部节点对应左上角正方形的两个三角形单元，上左，左上
fori=2:
Imax
n1=n1+1;
NE（1,n1）=NN（i,i）;
NE（2,n1）=NN（i-1,i）;
NE（3,n1）=NN（i-1,i-1）;
-
n1=n1+1;
NE（1,n1）=NN（i,i）;
NE（2,n1）=NN（i-1,i+1）;
NE（3,n1）=NN（i-1,i）;
n1=n1+1;
NE（1,n1）=NN（i,i）;
NE（2,n1）=NN（i,i+1）;
NE（3,n1）=NN（i-1,i+1）;
）
end%下斜边界节点要对应三个单元，上左，左上，左下
n1=n1+1;
NE（1,n1）=NN（Imax+1,Imax+1）;
NE（2,n1）=NN（Imax,Imax+1）;
NE（3,n1）=NN（Imax,Imax）;%右顶点的左下
n1=n1+1;
NE（1,n1）=NN（Imax+1,Imax+1）;
NE（2,n1）=NN（Imax,Imax+2）;
:
NE（3,n1）=NN（Imax,Imax+1）;%右顶点的左上
K=0;
fori=Imax:
-1:
2
K=K+2;
n1=n1+1;
NE（1,n1）=NN（i,i+K）;
NE（2,n1）=NN（i-1,i+K+1）;
NE（3,n1）=NN（i-1,i+K）;
）
end%右上边界的左上
nel=n1;%此时n1值为总的单元个数
求解偏微分方程
function[c,f,s]=pdefun（x,t,u,du）
c=[1;1];
f=[*du
（1）;*du
（2）];
temp=u
（1）-u
（2）;
s=[-1;1].*（exp*temp）-exp*temp））;
functionu0=pdeic（x）
function
[pa,qa,pb,qb]=pdebc（xa,ua,xb,ub,t）
%a表示左边界，b表示右边界
pa=[0;ua
（2）];qa=[1;0];
pb=[ub
（1）-1;0];qb=[0;1];
cl
x=0:
:
1;
t=0:
:
2;
m=0;
sol=pdepe（m,@pdefun,@pdeic,@pdebc,x,t）;
figure（'numberite'off'name';'PDEDemo--byMatlabsky'）%创建个窗口，窗口名字
是name后边的名字'NumberTitle'off'是关掉默认显示名字
title（'TheSolutionofu_1'）
xlabel（'X'）
ylabel（'T）
subplot（212）
st+f:
+,..;,）
title（'TheSolutionofu_2"）
xlabel（X'）
ylabel（'T）
zlabel（"U'）
升级会员 