This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Inhomogeneous Heat Equation on Square Domain

This example shows how to solve the heat equation with a source term.

The basic heat equation with a unit source term is


This equation is solved on a square domain with a discontinuous initial condition and zero temperatures on the boundaries.

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 of the material.


Specify internal heat source.


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'

Set zero temperatures on all four outer edges of the square.


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 a mesh.

msh = generateMesh(thermalmodel);
axis equal

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);
99 successful steps
0 failed attempts
200 function evaluations
1 partial derivatives
20 LU decompositions
199 solutions of linear systems
T = result.Temperature;

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.