fd1d_heat_steady, a FORTRAN90 code which applies the finite difference method to estimate the solution of the steady state (time independent) 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. 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) )
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 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 code assumes that the user will provide a calling code, and functions to evaluate K(X) and F(X).
call fd1d_steady_heat ( n, a, b, ua, ub, k, f, x, u )
where
The computer code and data files described and made available on this web page are distributed under the MIT license
fd1d_heat_steady is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version
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 FORTRAN90 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 FORTRAN90 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 FORTRAN90 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 FORTRAN90 code which uses the finite difference method and explicit time stepping to solve the time dependent heat equation in 1D.
FD1D_HEAT_IMPLICIT, a FORTRAN90 code which uses the finite difference method and implicit time stepping to solve the time dependent heat equation in 1D.
FD1D_PREDATOR_PREY, a FORTRAN90 code which implements a finite difference algorithm for a predator-prey system with spatial variation in 1D.
FD1D_WAVE, a FORTRAN90 code which applies the finite difference method to solve the time-dependent wave equation utt = c * uxx in one spatial dimension.
FD2D_HEAT_STEADY, a FORTRAN90 code which uses the finite difference method (FDM) to solve the steady (time independent) heat equation in 2D.
FEM1D_HEAT_STEADY, a FORTRAN90 code which uses the finite element method to solve the steady (time independent) heat equation in 1D.
FEM2D_HEAT, a FORTRAN90 code which solves the 2D time dependent heat equation on the unit square.