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

1、% 用有限元法求解三角形形區(qū)域上的possion方程function finite_element_tri(imax,jmax)global ndm nel na% ndm 總節(jié)點(diǎn)數(shù)% nel 基元數(shù)% na 活動(dòng)節(jié)點(diǎn)數(shù)imax=30;jmax=60;%設(shè)定網(wǎng)格數(shù)v=0; j=0;x0=1/imax;y0=x0;domain_tri x,y,nn,ne=setelm_tri(imax,jmax); % 給節(jié)點(diǎn)和三角形元素編號(hào)，并設(shè)定節(jié)點(diǎn)坐標(biāo)% 求解有限元方程的求系數(shù)矩陣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);% 求有限元方程的右端項(xiàng)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% 求解方程得結(jié)果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% 畫(huà)等電勢(shì)線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% 畫(huà)解函數(shù)的曲面圖figure(2)surf(x1,y1,nnv);fid=fopen(finite_element_tri.txt,w);fprintf(fid,n *有限元法求解三角形區(qū)域上po

5、ssion方程的結(jié)果* n n);fprintf(fid,n 節(jié)點(diǎn)編號(hào) n n);nna=fliplr(nn);fprintf(fid,%4d%4d%4d%4d%4d%4d%4d%4d%4dn,nna);fprintf(fid,n 各節(jié)點(diǎn)的電勢(shì) 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 節(jié)點(diǎn)編號(hào) 坐標(biāo)分量x 坐標(biāo)分量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)% 給節(jié)點(diǎn)和三角形元素編號(hào)，并設(shè)定節(jié)點(diǎn)坐標(biāo)global ndm nel na% i1 i2 j1 j2 imax jmax分別描述網(wǎng)線縱向和橫向數(shù)目的變量% x y表示節(jié)點(diǎn)坐標(biāo)% nn描述節(jié)點(diǎn)編號(hào)% ne 描述各基點(diǎn)局域節(jié)點(diǎn)的矩陣% ndm 總節(jié)點(diǎn)數(shù)% nel 基元數(shù)% na 表示活動(dòng)節(jié)點(diǎn)數(shù)nlm=imax*jmax;dx=1/imax;dy=1/jmax;x=nlm:1;y=nlm:1;nn=zeros

7、(imax+1,jmax+1);n1=0; % 活動(dòng)節(jié)點(diǎn)編號(hào)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% 以上為對(duì)節(jié)點(diǎn)進(jìn)行編號(hào)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