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

FD1D_HEAT_STEADY is a FORTRAN77 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 calling a function to solve the system A*U=RHS will give us the solutioni.

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:

        call fd1d_steady_heat ( n, a, b, ua, ub, k, f, x, u )
      

n,

the number of spatial points.

a, b,

the left and right endpoints.

ua, ub,

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 X coordinates of 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.

Related Data and Programs:

FD1D_BVP is a FORTRAN77 program which applies the finite difference method to a two point boundary value problem in one spatial dimension.

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

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

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

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

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

FEM1D_HEAT is a FORTRAN90 program which uses the finite element method to solve the 1D Time Dependent Heat Equation.

FEM2D_HEAT is a FORTRAN90 program which solves the 2D time dependent heat equation on the unit square.

FFH_SPARSE is a MATLAB program which solves the time dependent heat equation in an arbitrary triangulated 2D region, using MATLAB's sparse matrix storage format and solver.

FREE_FEM_HEAT is a FORTRAN90 program which uses the finite element method and the backward Euler method to solve the 2D time-dependent heat equation on an arbitrary triangulated region.

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.f, the library.
• fd1d_heat_steady.csh, commands to compile the library.

Examples and Tests:

• problem1.f, uses K(X) = 1, F(X) = 0, so the solution should be the straight line that connects the boundary values.
• problem1.csh, commands to compile the problem and run it with the library.
• problem1_x.txt, the coordinates of the nodes.
• problem1_u.txt, the computed temperatures at the nodes.
• problem2.f, 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.csh, commands to compile the problem and run it with the library.
• problem2_x.txt, the coordinates of the nodes.
• problem2_u.txt, the computed temperatures at the nodes.
• problem3.f, uses K(X) = 1, F(X) defines a heat source, so the solution can rise above the boundary values.
• problem3.csh, commands to compile the problem and run it with the library.
• problem3_x.txt, the coordinates of the nodes.
• problem3_u.txt, the computed temperatures at the nodes.
• problem4.f, 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.csh, commands to compile the problem and run it with the library.
• problem4_x.txt, the coordinates of the nodes.
• problem4_u.txt, the computed temperatures at the nodes.

List of Routines:

• DTABLE_DATA_WRITE writes DTABLE data to a file.
• DTABLE_WRITE writes a DTABLE to a file.
• FD1D_HEAT_STEADY solves the steady 1D heat equation.
• GET_UNIT returns a free FORTRAN unit number.
• R83_NP_FS factors and solves a R83 system.
• S_BLANK_DELETE removes blanks from a string, left justifying the remainder.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• TIMESTRING writes the current YMDHMS date into a string.

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