基于matlab的有限元法求解泊松方程

1、% 用有限元法求解三角形形区域上的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)

2、; 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; endendm=t(1:na,1:na);% 求有限元方程的右端项g=z

3、eros(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; endend% 求解方程得结果f=mg; nnv=zeros(imax+1,jmax+1);fi=zeros(ndm,1);fi(1:na)=f(1:na);fi(na+1:ndm)

4、=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); endend% 画等电势线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;end% 画解函数的曲面图figure(2)surf(x1,y1,nnv);fid=fopen(finite_element_tri.txt,w);fprintf(fid,n *有限元法求解三角形区域上po

5、ssion方程的结果* n n);fprintf(fid,n 节点编号 n n);nna=fliplr(nn);fprintf(fid,%4d%4d%4d%4d%4d%4d%4d%4d%4dn,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.6fn,nnv);l=1:ndm;fprintf(fid,nn 节点编号 坐标分量x 坐标分量y u（x,y）的值nn);for i=1:ndm fprintf(fid,%8d%14.

6、5f%14.5f%14.5fn,l(i),x(i),y(i),fi(i);endfclose(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

7、(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; endendk=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

8、)=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)*dy;endfor 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;

9、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); 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; n

10、e(1,n1)=nn(i,j); ne(2,n1)=nn(i,j+1); ne(3,n1)=nn(i-1,j+1); endendfor 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);endn1=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);endnel=n1;end function