MATLAB Examples

Inhomogeneous Heat Equation on a Square Domain

This example shows how to solve the heat equation with a source term using 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 temperatures on the boundaries.


Create Thermal Model and Geometry

Create a transient thermal model.

thermalmodel = createpde('thermal','transient');

Create a square geometry centered at x = 0 and y = 0 with sides of length 2. Include a circle of radius 0.4 concentric with the square.

R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
C1 = [1;0;0;0.4];
C1 = [C1;zeros(length(R1) - length(C1),1)];
gd = [R1,C1];
sf = 'R1+C1';
ns = char('R1','C1')';
g = decsg(gd,sf,ns);

Append the geometry to the model.


Specify Thermal Properties and Internal Heat Source

Specify thermal properties of the material.


Specify internal heat source.


Set Zero Temperatures on Boundaries

Plot the geometry and display the edge labels for use in the boundary condition definition.

axis([-1.1 1.1 -1.1 1.1]);
axis equal
title 'Geometry With Edge and Subdomain Labels'

Temperature is zero at all four outer edges of the square.


Specify Initial Temperature

The discontinuous initial value is 1 inside the circle and zero outside. Specify zero initial temperature everywhere.

% Specify non-zero initial temperature inside the circle (Face 2).

Generate and Plot Mesh

msh = generateMesh(thermalmodel);
axis equal

Find and Plot Solution

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

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

thermalmodel.SolverOptions.ReportStatistics ='on';
result = solve(thermalmodel,tlist);
T = result.Temperature;
98 successful steps
0 failed attempts
198 function evaluations
1 partial derivatives
21 LU decompositions
197 solutions of linear systems

Plot the solution.

Tmax = max(max(T));
Tmin = min(min(T));
for j = 1:nframes,
    caxis([Tmin Tmax]);
    axis([-1 1 -1 1 0 1]);
    Mv(j) = getframe;

To play the animation, use the movie(Mv,1) command.