matlab程序解泊松方程.docx

1、matlab程序解泊松方程求解泊松方程的function Finite_element_tri(Imax)% 用有限元法求解三角形形区域上的Possion方程Jmax=2*Imax;% 其中Imax Jmax分别表示x轴和y轴方向的网格数，其中Jmax等于Imax的两倍% 定义一些全局变量global ndm nel na % ndm 总节点数% nel 基元数% na 表示活动节点数V=0; J=0;X0=1/Imax;Y0=X0;%V=0为边界条件domain_tri % 调用函数画求解区域X,Y,NN,NE=setelm_tri(Imax,Jmax); % 给节点和三角形元素编号，并设定

2、节点坐标% 以下求解有限元方程的求系数矩阵T=zeros(ndm,ndm);for n=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;%三角形面积 for k=1:3 if n1=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(n

3、2)-Y(n3)2+(X(n3)-X(n2)2)/(4*s);%V=0则边界积分为零，非零时积分编程类似，再加边界积分。 end| k=n1;n1=n2;n2=n3;n3=k; % 轮换坐标将值赋入3阶主子矩阵中 endendM=T(1:na,1:na);% 求有限元方程的右端项f=X;%场源函数G=zeros(na,1);for n=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; for k=1:3 if n1=na G(n1)=G(n1

4、)+(2*f(n1)+f(n2)+f(n3)*s/12;%f在单元上为线性差值时场域单元的积分公式 end n4=n1; n1=n2; n2=n3; n3=n4; % 轮换坐标标 endendF=MG; % 求解方程得结果NNV=zeros(Imax+1,Jmax+1);fi=zeros(ndm,1);fi(1:na)=F(1:na);fi(na+1:ndm)=V;for j=0:Jmax for i=0:Imax( n=NN(i+1,j+1); if n=0 n=na+1; end NNV(i+1,j+1)=fi(n); endendfigure(2)imagesc(NNV);%画解函数的平

5、面图X1=zeros(1,Imax+1);Y1=zeros(1,Jmax+1);for i=1:Imax+1 X1(i)=(i-1)*X0;endfor i=1:Jmax+1 Y1(i)=(i-1)*Y0;:endfigure(3)surf(X1,Y1,NNV);% 画解函数的曲面图% 以下是结果的输出fid=fopen(,w);fprintf(fid,n *有限元法求解三角形区域上Possion方程的结果* n n);L=1:ndm;fprintf(fid,nn 节点编号 坐标分量x 坐标分量y u（x,y）的值nn);for i=1:ndm fprintf(fid,%8d%n,L(i),X

6、(i),Y(i),fi(i);endfclose(fid);function domain_tri% 画求解区域xy=0 1;0 -1;1 0;-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) % 给节点和三角形单元编号，并设定节点坐标 % 定义一些全局变量global ndm nel na

7、% I1 I2 J1 J2 Imax Jmax分别描述网线纵向和横向数目的变量% X Y表示节点坐标% 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; for j=3:Jmax/2 for i=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; en

8、d.endk=Jmax/2+1;for j=Jmax/2+1:Jmax-1 %三角形区域上下两部分节点坐标分别求 k=k-1; for i=2:k n1=n1+1; NN(i,j)=n1; X(n1)=(i-1)*dx; Y(n1)=1+(j-Jmax-1)*dy; endendna=n1;%不含边界节点数for j=Jmax+1:-1:Jmax/2+1 %降序 n1=n1+1; NN(1,j)=n1; X(n1)=0; Y(n1)=1+(j-Jmax-1)*dy;endfor j=Jmax/2:-1:1 n1=n1+1; NN(1,j)=n1; X(n1)=0; Y(n1)=-1+(j-1)

9、*dy;end %-for i=2:Imax+1 n1=n1+1; NN(i,i)=n1; X(n1)=(i-1)*dx; Y(n1)=-1+(i-1)*dy;endK=0;for i=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;for j=3:Jmax/2 for i=2:j-1 n1=n1+1; NE(1,n1)=NN(i,j); NE(2,n1)=NN(i-1,j+1)

10、; 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); endendk=Jmax/2+1;for j=Jmax/2+1:Jmax-1 k=k-1;（ for i=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); endend %内部节

11、点对应左上角正方形的两个三角形单元，上左，左上for i=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,Im

12、ax+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;for i=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值为总的单元个数求解偏微分方程funct

13、ion 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);function u0=pdeic(x)functionpa,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;clx=0:1;t=0:2;m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t);figure(numberiteoff name;PDE Demo- -by Matlabsky)%创建个窗口，窗口名字是name后边的名字NumberTitleoff是关掉默认显示名字title(The Solution ofu_ 1)xlabel(X)ylabel(T)subplot(212)st+f:+,.;,)title(The Solutionofu_ 2)xlabel(X)ylabel(T)zlabel(U)