Lab 1 -- Solving a heat equation in Matlab


  • To introduce students to programming, and Matlab programming in particular.
  • To learn how pseudocode is used as an intermediate step in converting an algorithm stated in English and Math into a computer program.
  • To operationalize calculus concepts covered in lecture.
  • To introduce finite difference methods for the numerical solution of differential equations.


Periodic Heat Diffusion in Subsurface Rocks

Consider the problem of determining the subsurface temperature fluctuations of rock as the result of daily or seasonal temperature variations. As the surface is heated or cooled, the heat will diffuse through the soil and rock. Mathematically, diffusion may be represented by the equation

where T(t,z) is the temperature at time t and depth z, and K is a constant which measures the "diffusivity" of the rock, i.e. how quickly heat moves through the rock. A typical value for K is 10^-6 m^2/s.

The above equation may be solved "analytically", i.e. T(t,z) may be written exactly as an equation. Matlab is more appropriate, however, for solving problems numerically, i.e. approximately as (a matrix or a map) of numbers.

In order to solve the system numerically, we must first choose an appropriate "space" to solve the problem in. We will have two axes, the vertical axes representing depth, from 0m to 100 m below the surface, and the horizontal axis will represent one year's worth of time.

Boundary Conditions

The heat equation is an example of what is known as a "partial differential equation." A differential equation is any equation in which a function (Temperature in time and space in this instance) is not represented directly, but via it's derivative. Partial indicates that there are at least two variables (time and space) in the derivatives. An ordinary DE only has one derivative.

In order to solve a PDE numerically, we need to specify boundary conditions. The four boundaries of our "space" are the surface (0m), the bottom depth (100m), an arbitrary beginnning time, and one year later.

The top boundary is the simplest, we will simply specify a temperature that varies seasonally, i.e.

This will create an average annual temperature of 15, an annual low of 5 and an annual high of 25 (time is measured in months, and we are not starting on Jan 1).

*Bottom boundary

The bottom boundary condition is easy as well... we will assume that no heat reaches the bottom from the top, i.e.

This specifies that there is 0 heat flow at the bottom. We do not yet know what this bottom temperature will be.

The left and right boundary conditions are a bit more problematic, as they require that we already know the temperature at every depth, which is exactly what we are trying to find out! But we do know that the temperature now and exactly one year from now should be the same.

Finite difference approximations

The final step is to determine how to use finite differences to evaluate the derivatives. For the first derivative of temperature with respect to time we will use the "forward derivative"

For the first derivative of temperature with respect to depth, we will use the "central derivative"

% For the second derivative of temperature with respect to depth, we will use

See your class notes for the logic use here.

The following little program implements this procedure

Matlab script

%%% Script to solve the heat equation in soil with
%%% seasonal temperature forcing

%%% Initialize variables
dz = .25; %each depth step is 1/4 meter
Nz = 400; % Choose the number of depth steps (should go to at least 100 m)
Nt = 5000; % Choose the number of time steps
dt = (365*24*60*60)/Nt; %Length of each time step in seconds (~ 6.3*10^3 seconds, or ~105 minutes)
K =  2*10^-6; % "canonical" K is 10e-6 m^2/s

T = 15*ones(Nz+1,Nt+1);  %Create temperature matrix with Nz+1 rows, and Nt+1 columns
                         % Initial guess is that T is 15 everywhere.
time = [0:12/Nt:12];
T(1,:) = 15-10*sin(2*pi*time/12);  %Set surface temperature

maxiter = 500
for iter = 1:maxiter
    Tlast = T; %Save the last guess
    T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp
    for i=2:Nt+1,
        depth_2D = (T(1:end-2,i-1)-2*T(2:end-1,i-1)+T(3:end,i-1))/dz^2;
        time_1D = K*depth_2D;
        T(2:end-1,i) = time_1D*dt + T(2:end-1,i-1);
        T(end,i) = T(end-1,i); % Enforce bottom BC
    err(iter) = max(abs(T(:)-Tlast(:)));  %Find difference between last two solutions
    if err(iter)<1E-4
        break; % Stop if solutions very similar, we have convergence
if iter==maxiter;
    warning('Convergence not reached')
maxiter =



Plot the solution

plot(log(err)), title('Convergence plot')
imagesc([0 12],[0 100],T); title('Temperature plot (imagesc)')
depth = [0:dz:Nz*dz];
contourf(time,-depth,T); title('Temperature plot (contourf)')

Plot the heat at depths of 0, 5 10, 15 and 20 m.

xlabel('Time (months)'); ylabel('Temperature (C)');
legend('0m','5m','10m','15m', '20m')

Plot the solution using contourf. Adjust your plot so that only the "interesting" depths are shown. Plot the solution at 0, 1/2, 1 and 1.5 meters below surface. About how many hours after the surface reaches its hottest, does it reach its hottest at 1/2 meters below the surface?