www.gusucode.com > pde 案例源码 matlab代码程序 > pde/ParabolicEquationExample.m

    %% Parabolic Equation
% Solve the parabolic equation
%
% $$\frac{\partial u}{\partial t} = \Delta u$$
%
% on the square domain specified by |squareg|.
%%
% Create a PDE model and import the geometry.
model = createpde;
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on')
ylim([-1.1,1.1])
axis equal
%%
% Set Dirichlet boundary conditions $u = 0$ on all edges.
applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geometry.NumEdges,'u',0);
%%
% Generate a relatively fine mesh.
generateMesh(model,'Hmax',0.02);
%%
% Set the initial condition to have $u(0) = 1$ on the disk $x^2 + y^2 \le
% 0.4^2$ and $u(0) = 0$ elsewhere.
p = model.Mesh.Nodes;
u0 = zeros(size(p,2),1); 
ix = find(sqrt(p(1,:).^2 + p(2,:).^2) <= 0.4); 
u0(ix) = ones(size(ix));
%%
% Set solution times to be from 0 to 0.1 with step size 0.005.
tlist = linspace(0,0.1,21);
%%
% Create the PDE coefficients.
c = 1;
a = 0;
f = 0;
d = 1;
%%
% Solve the PDE.
u = parabolic(u0,tlist,model,c,a,f,d);
%%
% Plot the initial condition, the solution at the final time, and two
% intermediate solutions.
figure
subplot(2,2,1)
pdeplot(model,'XYData',u(:,1));
axis equal
title('t = 0')
subplot(2,2,2)
pdeplot(model,'XYData',u(:,5))
axis equal
title('t = 0.02')
subplot(2,2,3)
pdeplot(model,'XYData',u(:,11))
axis equal
title('t = 0.05')
subplot(2,2,4)
pdeplot(model,'XYData',u(:,end))
axis equal
title('t = 0.1')