Inhomogeneous Heat Equation on a Square Domain

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

$$\frac{\partial u}{\partial t} - \Delta u = 1$$

This is solved on a square domain with a discontinuous initial condition and zero Dirichlet boundary conditions.

Problem Definition

  • 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;

Boundary Conditions

% 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);

Generate Mesh

[p,e,t] = initmesh(g);
figure;
pdemesh(p,e,t);
axis equal

Generate Initial Conditions

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));

Generate Time Discretization

We want the solution at 20 points in time between 0 and 0.1.

nframes = 20;
tlist = linspace(0,0.1,nframes);

Find FEM Solution

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

Plot FEM Solution

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);

Was this topic helpful?