有限元法求解泊松方程，用matlab编写的计算程序

% 用有限元法求解三角形形区域上的Possion方程 function Finite_element_tri(Imax,Jmax) global ndm nel na % ndm 总节点数 % nel 基元数 % na 活动节点数 Imax=30;Jmax=60;%设定网格数 V=0; J=0;X0=1/Imax;Y0=X0; domain_tri [X,Y,NN,NE]=setelm_tri(Imax,Jmax); % 给节点和三角形元素编号，并设定节点坐标 % 求解有限元方程的求系数矩阵 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(n2)-Y(n3))^2+(X(n3)-X(n2))^2)/(4*s); end k=n1;n1=n2;n2=n3;n3=k; end end M=T(1:na,1:na); % 求有限元方程的右端项 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)+(2*X(n1)+X(n2)+X(n3))*s/12; 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; 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); end end % 画等电势线 X1=zeros(1,Imax+1); Y1=zeros(1,Jmax+1); for i=1:Imax+1 X1(i)=(i-1)*X0; end for i=1:Jmax+1 Y1(i)=(i-1)*Y0; end % 画解函数的曲面图 figure(2) surf(X1,Y1,NNV'); fid=fopen('Finite_element_tri.txt','w'); fprintf(fid,'\n *********有限元法求解三角形区域上Possion方程的结果********** \n \n'); fprintf(fid,'\n 节点编号 \n \n'); Nna=fliplr(NN); fprintf(fid,'%4d%4d%4d%4d%4d%4d%4d%4d%4d\n',Nna);fprintf(fid,'\n 各节点的电势 \n \n'); NNV=fliplr(NNV); fprintf(fid,'%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f\n',NNV); L=[1:ndm]'; fprintf(fid,'\n\n 节点编号 坐标分量x 坐标分量y u（x,y）的值\n\n'); for i=1:ndm fprintf(fid,'%8d%14.5f%14.5f%14.5f\n',L(i),X(i),Y(i),fi(i)); end fclose(fid); end function [X,Y,NN,NE]=setelm_tri(Imax,Jmax) % 给节点和三角形元素编号，并设定节点坐标 global ndm nel na % 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(n1)=(i-1)*dx; Y(n1)=-1+(j-1)*dy; end end k=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; end end na=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; end for j=Jmax/2:-1:1 n1=n1+1; NN(1,j)=n1; X(n1)=0; Y(n1)=-1+(j-1)*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; end K=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); 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; 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); end end 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,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; 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; end function domain_tri % 定义求解区域 xy=[0 -1;-1 1;1 1;]; 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);