matlab程序解泊松方程Word格式.docx
《matlab程序解泊松方程Word格式.docx》由会员分享,可在线阅读,更多相关《matlab程序解泊松方程Word格式.docx(9页珍藏版)》请在冰豆网上搜索。
![matlab程序解泊松方程Word格式.docx](https://file1.bdocx.com/fileroot1/2022-11/29/4ea0835a-9858-4874-84bc-85be4583c989/4ea0835a-9858-4874-84bc-85be4583c9891.gif)
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
M=T(1:
na,1:
na);
%求有限元方程的右端项
f=X;
%场源函数
G=zeros(na,1);
—
=na
G(n1)=G(n1)+(2*f(n1)+f(n2)+f(n3))*s/12;
%f在单元上为线性差值时场域单元的积分公式
n4=n1;
n1=n2;
n2=n3;
n3=n4;
%轮换坐标标
》
F=M\G;
%求解方程得结果
NNV=zeros(Imax+1,Jmax+1);
fi=zeros(ndm,1);
fi(1:
na)=F(1:
fi(na+1:
ndm)=V;
forj=0:
Jmax
fori=0:
Imax
(
n=NN(i+1,j+1);
ifn<
=0
n=na+1;
NNV(i+1,j+1)=fi(n);
figure
(2)
imagesc(NNV);
%画解函数的平面图
X1=zeros(1,Imax+1);
Y1=zeros(1,Jmax+1);
fori=1:
Imax+1
X1(i)=(i-1)*X0;
Jmax+1
Y1(i)=(i-1)*Y0;
:
figure(3)
surf(X1,Y1,NNV'
);
%画解函数的曲面图
%以下是结果的输出
fid=fopen('
'
'
w'
fprintf(fid,'
\n*********有限元法求解三角形区域上Possion方程的结果**********\n\n'
L=[1:
ndm]'
;
\n\n节点编号坐标分量x坐标分量yu(x,y)的值\n\n'
^
ndm
fprintf(fid,'
%8d%%%\n'
L(i),X(i),Y(i),fi(i));
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描述各单元局部节点编号与总体编号对应的矩阵
%nel单元数
%na表示不含边界的节点数
nlm=Imax*Jmax;
dx=1/Imax;
dy=1/Jmax;
X=nlm:
1;
Y=nlm:
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;
.
k=Jmax/2+1;
forj=Jmax/2+1:
Jmax-1%三角形区域上下两部分节点坐标分别求
k=k-1;
k
[
Y(n1)=1+(j-Jmax-1)*dy;
na=n1;
%不含边界节点数
forj=Jmax+1:
-1:
Jmax/2+1%降序
NN(1,j)=n1;
X(n1)=0;
forj=Jmax/2:
1
end%
fori=2:
NN(i,i)=n1;
Y(n1)=-1+(i-1)*dy;
K=0;
fori=Imax:
2
K=K+2;
NN(i,i+K)=n1;
Y(n1)=1+(i+K-Jmax-1)*dy;
%以上四个循环为对边界节点进行编号
ndm=n1;
《
NE=zeros(3,2*ndm);
n1=0;
Jmax/2
NE(1,n1)=NN(i,j);
NE(2,n1)=NN(i-1,j+1);
NE(3,n1)=NN(i-1,j);
NE(2,n1)=NN(i,j+1);
NE(3,n1)=NN(i-1,j+1);
Jmax-1
end%内部节点对应左上角正方形的两个三角形单元,上左,左上
NE(1,n1)=NN(i,i);
NE(2,n1)=NN(i-1,i);
NE(3,n1)=NN(i-1,i-1);
NE(2,n1)=NN(i-1,i+1);
NE(3,n1)=NN(i-1,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);
%右顶点的左下
NE(2,n1)=NN(Imax,Imax+2);
NE(3,n1)=NN(Imax,Imax+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;
qb=[0;
cl
x=0:
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:
+,..;
)
TheSolutionofu_2"
xlabel(X'
zlabel("
U'