matlab有限元解二维抛物方程

%%%%% 方程 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

  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

  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

 ：，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)).*cos(pi*pxy(:,2)).*sin((j-1)*r)]’,’pxy’);

H1_err=getH1error((node,elem,Dut(:,j),u(:,j),5);