FD1D_HEAT_STEADY
Finite Difference Solution of a 1D Steady State Heat Equation

FD1D_HEAT_STEADY is a MATLAB program 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)
      

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. You can 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) )
      

       - ( + 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 )

• 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.

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

Languages:

FD1D_HEAT_STEADY is available in a C version and a C++ version and a FORTRAN77 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_BURGERS_LAX, a MATLAB program 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 program 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 program which applies the finite difference method to a two point boundary value problem in one spatial dimension.

FD1D_DISPLAY, a MATLAB program which reads a pair of files defining a 1D finite difference model, and plots the data.

FD1D_HEAT_EXPLICIT, a MATLAB program which uses the finite difference method and explicit time stepping to solve the time dependent heat equation in 1D.

FD1D_HEAT_IMPLICIT, a MATLAB program which uses the finite difference method and implicit time stepping to solve the time dependent heat equation in 1D.

FD1D_PREDATOR_PREY, a MATLAB program which implements a finite difference algorithm for predator-prey system with spatial variation in 1D.

FD1D_PREDATOR_PREY_PLOT, a MATLAB program which plots the output from the FD1D_PREDATOR_PREY program;

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

FEM_50_HEAT, a MATLAB program which implements a finite element calculation for the heat equation.

FEM2D_HEAT, a MATLAB program 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:

• fd1d_heat_steady.m, the program.
• r8mat_write.m, writes an R8MAT to a file.
• timestamp.m, prints the YMDHMS date as a timestamp.

Examples and Tests:

• problem1.m, uses K(X) = 1, F(X) = 0, so the solution should be the straight line that connects the boundary values.
• problem1_nodes.txt, the coordinates of the nodes.
• problem1_values.txt, the computed temperatures at the nodes.
• problem1.png, a PNG image of the solution.

• problem2.m, uses K(X) which is set to different constants over three subregions, and F(X) = 0.0, so the solution will be a piecewise linear function that connects the boundary values.
• problem2_nodes.txt, the coordinates of the nodes.
• problem2_values.txt, the computed temperatures at the nodes.
• problem2.png, a PNG image of the solution.

• problem3.m, uses K(X) = 1, F(X) defines a heat source, so the solution can rise above the boundary values.
• problem3_nodes.txt, the coordinates of the nodes.
• problem3_values.txt, the computed temperatures at the nodes.
• problem3.png, a PNG image of the solution.

• problem4.m, uses K(X) = 1, F(X) defines a heat source and a heat sink, so the solution can go above and below the boundary values.
• problem4_nodes.txt, the coordinates of the nodes.
• problem4_values.txt, the computed temperatures at the nodes.
• problem4.png, a PNG image of the solution.

You can go up one level to the MATLAB source codes.