matlab有限元解二维抛物方程

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

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

本文来源：https://www.bwwdw.com/article/um1p.html