matlab程序(解泊松方程)matlab程序(解泊松方程)matlab程序(解泊松方程)matlab程序(解泊松方程)編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:求解泊松方程的functionFinite_element_tri(Imax)%用有限元法求解三角形形區(qū)域上的Possion方程Jmax=2*Imax;%其中ImaxJmax分別表示x軸和y軸方向的網(wǎng)格數(shù)，其中Jmax等于Imax的兩倍%定義一些全局變量globalndmnelna%ndm總節(jié)點數(shù)%nel基元數(shù)%na表示活動節(jié)點數(shù)V=0;J=0;X0=1/Imax;Y0=X0;%V=0為邊界條件domain_tri%調(diào)用函數(shù)畫求解區(qū)域[X,Y,NN,NE]=setelm_tri(Imax,Jmax);%給節(jié)點和三角形元素編號，并設(shè)定節(jié)點坐標(biāo)%以下求解有限元方程的求系數(shù)矩陣T=zeros(ndm,ndm);forn=1:neln1=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;%三角形面積fork=1:3ifn1<=na|n2<=naT(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則邊界積分為零，非零時積分編程類似，再加邊界積分。endk=n1;n1=n2;n2=n3;n3=k;%輪換坐標(biāo)將值賦入3階主子矩陣中endendM=T(1:na,1:na);%求有限元方程的右端項f=X;%場源函數(shù)G=zeros(na,1);forn=1:neln1=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;fork=1:3ifn1<=naG(n1)=G(n1)+(2*f(n1)+f(n2)+f(n3))*s/12;%f在單元上為線性差值時場域單元的積分公式endn4=n1;n1=n2;n2=n3;n3=n4;%輪換坐標(biāo)標(biāo)endendF=M\G;%求解方程得結(jié)果NNV=zeros(Imax+1,Jmax+1);fi=zeros(ndm,1);fi(1:na)=F(1:na);fi(na+1:ndm)=V;forj=0:Jmaxfori=0:Imaxn=NN(i+1,j+1);ifn<=0n=na+1;endNNV(i+1,j+1)=fi(n);endendfigure(2)imagesc(NNV);%畫解函數(shù)的平面圖X1=zeros(1,Imax+1);Y1=zeros(1,Jmax+1);fori=1:Imax+1X1(i)=(i-1)*X0;endfori=1:Jmax+1Y1(i)=(i-1)*Y0;endfigure(3)surf(X1,Y1,NNV');%畫解函數(shù)的曲面圖%以下是結(jié)果的輸出fid=fopen('','w');fprintf(fid,'\n*********有限元法求解三角形區(qū)域上Possion方程的結(jié)果**********\n\n');L=[1:ndm]';fprintf(fid,'\n\n節(jié)點編號坐標(biāo)分量x坐標(biāo)分量yu（x,y）的值\n\n');fori=1:ndmfprintf(fid,'%8d%%%\n',L(i),X(i),Y(i),fi(i));endfclose(fid);functiondomain_tri%畫求解區(qū)域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)%給節(jié)點和三角形單元編號，并設(shè)定節(jié)點坐標(biāo)%定義一些全局變量globalndmnelna%I1I2J1J2ImaxJmax分別描述網(wǎng)線縱向和橫向數(shù)目的變量%XY表示節(jié)點坐標(biāo)%NN描述節(jié)點編號%NE描述各單元局部節(jié)點編號與總體編號對應(yīng)的矩陣%ndm總節(jié)點數(shù)%nel單元數(shù)%na表示不含邊界的節(jié)點數(shù)nlm=Imax*Jmax;dx=1/Imax;dy=1/Jmax;X=nlm:1;Y=nlm:1;NN=zeros(Imax+1,Jmax+1);n1=0;forj=3:Jmax/2fori=2:j-1n1=n1+1;NN(i,j)=n1;%X=i列，Y=j行處節(jié)點X(n1)=(i-1)*dx;Y(n1)=-1+(j-1)*dy;endendk=Jmax/2+1;forj=Jmax/2+1:Jmax-1%三角形區(qū)域上下兩部分節(jié)點坐標(biāo)分別求k=k-1;fori=2:kn1=n1+1;NN(i,j)=n1;X(n1)=(i-1)*dx;Y(n1)=1+(j-Jmax-1)*dy;endendna=n1;%不含邊界節(jié)點數(shù)forj=Jmax+1:-1:Jmax/2+1%降序n1=n1+1;NN(1,j)=n1;X(n1)=0;Y(n1)=1+(j-Jmax-1)*dy;endforj=Jmax/2:-1:1n1=n1+1;NN(1,j)=n1;X(n1)=0;Y(n1)=-1+(j-1)*dy;end%fori=2:Imax+1n1=n1+1;NN(i,i)=n1;X(n1)=(i-1)*dx;Y(n1)=-1+(i-1)*dy;endK=0;fori=Imax:-1:2K=K+2;n1=n1+1;NN(i,i+K)=n1;X(n1)=(i-1)*dx;Y(n1)=1+(i+K-Jmax-1)*dy;end%以上四個循環(huán)為對邊界節(jié)點進(jìn)行編號ndm=n1;NE=zeros(3,2*ndm);n1=0;forj=3:Jmax/2fori=2:j-1n1=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;forj=Jmax/2+1:Jmax-1k=k-1;fori=2:kn1=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%內(nèi)部節(jié)點對應(yīng)左上角正方形的兩個三角形單元，上左，左上fori=2:Imaxn1=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%下斜邊界節(jié)點要對應(yīng)三個單元，上左，左上，左下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;fori=Imax:-1:2K=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值為總的單元個數(shù)求解偏微分方程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;0];qb=[0;1];clx=0::1;t=0::2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pde