matlab有限元解二维抛物方程.docVIP

%%%%% 真解 u=sin(pi*x)*sin(pi*y)*sin(t) %%%%% 方程 diff(u,t)-Laplace(u)=f %%%%% f=sin(pi*x)*sin(pi*y)*cos(t)+2*pi^2*sin(pi*x)*sin(pi*y)*sin(t) %clear all % clc %%%%finite element code for parabolic equation with constant coefficient %%%mesh%% node=[0,0;1,0;1,1;0,1]; elem=[2,3,1;4,1,3]; T=1; bdEdge=setboundary(node,elem,’Dirichlet’); n=input(‘Please input initial mesh:’); M=input(‘M=’); for i=1:n [node,elem,bdEdge]=uniformrefine(node,elem,bdEdge); end N=size(node,1); NT=size(elem,1); S=1/NT; r=1/M; A=zeros(N,N); u=zeros(N,M+1); u1=zeros(N,1); f=inline(‘sin(pi*xx(1,1))*sin(pi*xx(1,2))*cos(t)+2*pi^2*sin(pi*xx(1,1))*sin(pi*xx(1,2))*sin(t)’,’xx’,’t’); [lambda,weight]=quadpts(5); p=node’; T=elem’; totalEdge=[elem(:,[2,3]);elem(:,[3,1]);elem(:,[1,2])]; isBdEdge=reshape(bdEdge,3*NT,1); Dirichlet=totalEdge(isBdEdge==1),:); isBdNode=false(N,1); isBdNode(Dirichlet)=true; bdNode=find(isBdNode); freeNode=find(~isBdNode); for j=2:M+1 for i=1:NT F=zeros(N,1); F_ele=zeros(1,3); T_ele=elem(i,:); for m=1:7 xx(m,1)=(p(1,t(1,i))-p(1,t(3,i)))*lambda(m,1)+(p(1,t(2,i))-p(1,t(3,i))*lambda(m,2)+p(1,t(3,i)); xx(m,1)=(p(2,t(1,i))-p(2,t(3,i)))*lambda(m,1)+(p(2,t(2,i))-p(2,t(3,i))*lambda(m,2)+p(2,t(3,i)); end for i=1:3 for k=1:7 F_ele(i)=F_ele(i)+S*weight(k)*lambda(k,i)*f(xx(k,:),(j-1)*r); end end x=node(T_ele,:); [a,b]=Basis_coeff(x); A_ele=[a,b]’*[a,b]/(4*s); B_ele=zeros(3,3); for i=1:3 for j=1:3 if i==j B_ele(i,j)=1/12; else B_ele(i,j)=1/24; end end end A(T_ele,T_ele)=A(T_ele,T_ele)+S*B_ele+r*A_ele; F(T_ele,1)=F(T_ele,1)+r*F_ele’+S*B_ele; end uj=zeros(N,1); uj(freeNode)=A(freeNode,freeNode)\F(freeNode); u(:,j)=uj; end showresult(node,elem,uj) u_exact=zeros(N,M+1); for j=1:M+1 u_exact(：，j)=inline(‘sin(pi*pxy(:,1)).*sin(pi*pxy(:,2)).*sin((j-1)*r)’,’pxy’); end L2_err=getL2error(node,elem,u_exact(:,j),u(:,j),5); %t=(j-1)*r L’2_err=getL2error(node,elem,u_exact,u,5); Du(:,j)=inline(‘[pi*cos(pi*pxy(:,1)).*sin(pi*pxy(:,2)).*sin((j-1)*r),pi*sin(pi*pxy(:,1)