This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the 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?