This is machine translation

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

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Inhomogeneous Heat Equation on a Square Domain

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

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 to specify initial conditions.

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

Create a PDE Model with a single dependent variable

numberOfPDE = 1;
pdem = createpde(numberOfPDE);

Create a geometry and append it to the PDE Model


Apply Boundary Conditions

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'

% Solution is zero at all four outer edges of the square

Specify PDE Coefficients


Specify Initial Conditions

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

% Specify non-zero initial condition inside the circle, face ID 2.

Generate Mesh

msh = generateMesh(pdem);
axis equal

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 Solution with SolverOptions.ReportStatistics of pdem set to 'on'.

pdem.SolverOptions.ReportStatistics ='on';
result = solvepde(pdem,tlist);
u1 = result.NodalSolution;
75 successful steps
0 failed attempts
152 function evaluations
1 partial derivatives
16 LU decompositions
151 solutions of linear systems

Plot Solution

umax = max(max(u1));
umin = min(min(u1));
for j = 1:nframes,
    caxis([umin umax]);
    axis([-1 1 -1 1 0 1]);
    Mv(j) = getframe;

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

Was this topic helpful?