有限元解二维变系数抛物方程

%%%%% 真解u=sin(pi*x)*sin(pi*y)*sin(t)
%%%%% 方程diff(u,t)-div(a(x,y)*grad(u))=f
%%%%% a(x,y)=x^2+y^2+1
%%%%%
f=sin(pi*x)*sin(pi*y)*cos(t)-2*pi^2*sin(pi*x)*sin(pi*y)*sin(t)*(x^2+y^2+1)+2*pi*sin(t)*(x*cos( pi*x)*sin(pi*y)+y*sin(pi*x)*cos(pi*y))
%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)*(x x(1,1)^2+xx(1,2)^2+1)+2*pi*sin(t)*(xx(1,1)*cos(pi*xx(1,1))*sin(pi*xx(1,2))+y*sin(pi*xx(1,1))* cos(pi*xx(1,2));
a=inline(‘xx(1,1)^2+xx(1,2)^2+1’,’xx’);
[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
z=0;
for k=1:7
z=z+S*weight(k)*lambda(k)*a(xx(k,:));
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*z*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)\H(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);
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*px y(:,2)).*sin((j-1)*r)]’,’pxy’);
H1_err=getH1error((node,elem,Dut(:,j),u(:,j),5);。