fd1d_heat_steady, a MATLAB code which applies the finite difference method to estimate the solution of the steady state heat equation over a one dimensional region, which can be thought of as a thin metal rod.

We will assume the rod extends over the range A <= X <= B.

The quantity of interest is the temperature U(X) at each point in the rod.

We will assume that the temperature of the rod is fixed at known values at the two endpoints. Symbolically, we write the boundary conditions:

U(A) = UA
U(B) = UB

Inside the rod, we assume that a version of the steady (time independent) heat equation applies. This assumes that the situation in the rod has "settled down", so that the temperature configuration has no further tendency to change. The equation we will consider is

- d/dx ( K(X) * d/dx U(X) ) = F(X)

Here, the right hand side term F(X) allows us to consider internal heat sources in the metal - perhaps a portion of the rod is sitting above a blow torch. The coefficient K(X) is a measure of heat conductivity. It measures the rate at which the heat from a local hot spot will spread out.

If the heat source function F(X) is zero everywhere, and if K(X) is constant, then the solution U(X) will be the straight line function that matches the two endpoint values. Making F(X) positive over a small interval will "heat up" that portion. We simulate a rod that is divided into regions of different materials by setting the function K(X) to have a given value K1 over some subinteral of [A,B] and value K2 over the rest of the region.

To estimate the value of the solution, we will pick a uniform mesh of N points X(1) through X(N), from A to B. At the I-th point, we will compute an estimated temperature U(I). To do this, we will need to use the boundary conditions and the differential equation.

Since X(1) = A and X(N) = B, we can use the boundary conditions to set U(1) = UA and U(N) = UB.

For the points in the interior, we need to approximate the differential equation in a way that allows us to determine the solution values. We will do this using a finite difference approximation.

Suppose we are working at node X(I), which is associated with U(I). Then a centered difference approximation to

- d/dx ( K(X) * d/dx U(X) )

is
- ( + K(X(I)+DX/2) * ( U(X(I+1)) - U(X(I))   ) / DX )
- K(X(I)-DX/2) * ( U(X(I))   - U(X(I-1)) ) / DX ) / DX

If we rearrange the terms in this approximation, and set it equal to F(X(I)), we get the finite difference approximation to the differential equation at X(I):

- K(X(I)-DX/2)                   * U(X(I-1)
+ (   K(X(I)-DX/2) + K(X(I)+DX(2)) ) * U(X(I))
- K(X(I)+DX(2))   * U(X(I+1)) = DX * DX * F(X(I))

This means that we have N-2 equations, each of which involves an unknown U(I) and its left and right neighbors, plus the two boundary conditions, which give us a total of N equations in N unknowns.

We can set up and solve these linear equations using a matrix A for the coefficients, and a vector RHS for the right hand side terms, and the simple MATLAB command u = A \ rhs will give us the solution!

Because finite differences are only an approximation to derivatives, this process only produces estimates of the solution. Usually, we can reduce this error by decreasing DX.

This program assumes that the user will provide a calling program, and functions to evaluate K(X) and F(X).

Usage:

[ x, u ] = fd1d_steady_heat ( n, a, b, ta, tb, @k, @f )
where
• n, the number of spatial points.
• a, b, the left and right endpoints.
• ta, tb, the temperature values at the left and right endpoints.
• @k, the name of the function which evaluates K(X);
• @f, the name of the function which evaluates the right hand side F(X).
• x, output, the location of the nodes;
• u, output, the computed value of the temperature at the nodes.

Languages:

fd1d_heat_steady is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version

Related Data and Programs:

fd1d, a data directory which contains examples of 1d fd files, two text files that can be used to describe many finite difference models with one space variable, and either no time dependence or a snapshot at a given time;

fd1d_advection_ftcs, a MATLAB code which applies the finite difference method to solve the time-dependent advection equation ut = - c * ux in one spatial dimension, with a constant velocity, using the ftcs method, forward time difference, centered space difference.

fd1d_burgers_lax, a MATLAB code which applies the finite difference method and the lax-wendroff method to solve the non-viscous time-dependent burgers equation in one spatial dimension.

fd1d_burgers_leap, a MATLAB code which applies the finite difference method and the leapfrog approach to solve the non-viscous time-dependent burgers equation in one spatial dimension.

fd1d_bvp, a MATLAB code which applies the finite difference method to a two point boundary value problem in one spatial dimension.

fd1d_display, a MATLAB code which reads a pair of files defining a 1d finite difference model, and plots the data.

fd1d_heat_explicit, a MATLAB code which uses the finite difference method and explicit time stepping to solve the time dependent heat equation in 1d.

fd1d_heat_implicit, a MATLAB code which uses the finite difference method (fdm) and implicit time stepping to solve the time dependent heat equation in 1d.

fd1d_predator_prey, a MATLAB code which implements a finite difference algorithm for predator-prey system with spatial variation in 1d.

fd1d_predator_prey_plot, a MATLAB code which plots the output from fd1d_predator_prey();

fd1d_wave, a MATLAB code which applies the finite difference method to solve the time-dependent wave equation in one spatial dimension.

fem2d_heat, a MATLAB code which solves the 2d time dependent heat equation on the unit square.

Reference:

1. George Lindfield, John Penny,
Numerical Methods Using MATLAB,
Second Edition,
Prentice Hall, 1999,
ISBN: 0-13-012641-1,
LC: QA297.P45.

Source Code:

Last revised on 14 January 2019.