This example shows how to solve the heat equation with a source term using the parabolic function in the Partial Differential Equation Toolbox™.
The basic heat equation with a unit source term is
This is solved on a square domain with a discontinuous initial condition and zero Dirichlet boundary conditions.
|On this page…|
g: A specification function that is used by initmesh. For more information, please see the documentation page for squareg and pdegeom.
c, a, f, d: The coefficients of the PDE.
g = @squareg; c = 1; a = 0; f = 1; d = 1;
% Plot the geometry and display the edge labels for use in the boundary % condition definition. figure pdegplot(g, 'edgeLabels', 'on'); axis([-1.1 1.1 -1.1 1.1]); title 'Geometry With Edge Labels Displayed' % Create a pde entity for a PDE with a single dependent variable numberOfPDE = 1; pb = pde(numberOfPDE); % Create a geometry entity pg = pdeGeometryFromEdges(g); % Solution is zero at all four outer edges of the square pb.BoundaryConditions = pdeBoundaryConditions(pg.Edges(1:4), 'u', 0);
[p,e,t] = initmesh(g); figure; pdemesh(p,e,t); axis equal
The discontinuous initial value is 1 inside a circle of radius 0.4 and zero outside.
u0 = zeros(size(p,2),1); ix = find(sqrt(p(1,:).^2+p(2,:).^2)<0.4); u0(ix) = ones(size(ix));
We want the solution at 20 points in time between 0 and 0.1.
nframes = 20; tlist = linspace(0,0.1,nframes);
u1 = parabolic(u0,tlist,pb,p,e,t,c,a,f,d);
75 successful steps 1 failed attempts 154 function evaluations 1 partial derivatives 17 LU decompositions 153 solutions of linear systems
To speed up the plotting, we interpolate to a rectangular grid.
figure colormap(cool); x = linspace(-1,1,31); y = x; [~,tn,a2,a3] = tri2grid(p,t,u0,x,y); 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]); axis([-1 1 -1 1 0 1]); shading interp; Mv(j) = getframe; end movie(Mv,1);