matlab程序(解泊松方程)

求解泊松方程的functionFinite_element_tri(Imax)%用有限元法求解三角形形区域上的Possion方程Jmax=2*Imax;%其中ImaxJmax分别表示x轴和y轴方向的网格数，其中Jmax等于Imax的两倍%定义一些全局变量globalndmnelna%ndm总节点数%nel基元数%na表示活动节点数V=0;J=0;X0=1/Imax;Y0=X0;%V=0为边界条件domain_tri%调用函数画求解区域[X,Y,NN,NE]=setelm_tri(Imax,Jmax);%给节点和三角形元素编号，并设定节点坐标%以下求解有限元方程的求系数矩阵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;%轮换坐标将值赋入3阶主子矩阵中endendM=T(1:na,1:na);%求有限元方程的右端项f=X;%场源函数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;%轮换坐标标endendF=M\G;%求解方程得结果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);%画解函数的平面图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');%画解函数的曲面图%以下是结果的输出fid=fopen('Finite_element_tri.txt','w');fprintf(fid,'\n*********有限元法求解三角形区域上Possion方程的结果**********\n\n');L=[1:ndm]';fprintf(fid,'\n\n节点编号坐标分量x坐标分量yu（x,y）的值\n\n');fori=1:ndmfprintf(fid,'%8d%14.5f%14.5f%14.5f\n',L(i),X(i),Y(i),fi(i));endfclose(fid);functiondomain_tri%画求解区域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)%给节点和三角形单元编号，并设定节点坐标%定义一些全局变量globalndmnelna%I1I2J1J2ImaxJmax分别描述网线纵向和横向数目的变量%XY表示节点坐标%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;forj=3:Jmax/2fori=2:j-1n1=n1+1;NN(i,j)=n1;%X=i列，Y=j行处节点X(n1)=(i-1)*dx;Y(n1)=-1+(j-1)*dy;endendk=Jmax/2+1;forj=Jmax/2+1:Jmax-1%三角形区域上下两部分节点坐标分别求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;%不含边界节点数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%以上四个循环为对边界节点进行编号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%内部节点对应左上角正方形的两个三角形单元，上左，左上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%下斜边界节点要对应三个单元，上左，左上，左下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值为总的单元个数求解偏微分方程function[c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*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:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);figure('numberite'off'name';'PDEDemo--byMatlabsky')%创建个窗口，窗口名字是name后边的名字'NumberTitle'off'是关掉默认显示名字title('TheSolutionofu_1')xlabel('X')ylabel('T)subplot(212)st+f:+,..;,)title('TheSolutionofu_2)xlabel(X')ylabel('T)zlabel(U')