test10.m

matlab中PDE工具箱的帮助文件pde.pdf的例子全部用代码实现,花了我一个月的时间,有限元应用程序

%       We solve the standard heat equation with a source term
%        du/dt-div(grad(u))=1
%       on a square with a discontinuous initial condition.
pause % Strike any key to continue.
%       Problem definition
g='squareg'; % The unit square
b='squareb1'; % 0 on the boundary
c=1;
a=0;
f=1;
d=1;

%       Mesh
[p,e,t]=initmesh(g);

%       Initial condition: 1 inside the circle with radius 0.4.
%       0 otherwise.
u0=zeros(size(p,2),1);
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);
u0(ix)=ones(size(ix));
pause % Strike any key to continue.
%       Solve parabolic problem
u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);
% Time: 0.00526316
% Time: 0.0105263
% Time: 0.0157895
% Time: 0.0210526
% Time: 0.0263158
% Time: 0.0315789
% Time: 0.0368421
% Time: 0.0421053
% Time: 0.0473684
% Time: 0.0526316
% Time: 0.0578947
% Time: 0.0631579
% Time: 0.0684211
% Time: 0.0736842
% Time: 0.0789474
% Time: 0.0842105
% Time: 0.0894737
% Time: 0.0947368
% Time: 0.1
pause % Strike any key to continue.
%       To speed up the plotting, we interpolate to a rectangular grid.
x=linspace(-1,1,31);y=x;
[unused,tn,a2,a3]=tri2grid(p,t,u0,x,y);

%       Make the animation
newplot;
Mv = moviein(nframes);
umax=max(max(u1));
umin=min(min(u1));
for j=1:nframes,...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
  u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));...
  surf(x,y,u);caxis([umin umax]);colormap(cool),...
  axis([-1 1 -1 1 0 1]);...
  Mv(:,j) = getframe;...
end
% Show movie
movie(Mv,10)