FD1D_HEAT_EXPLICIT - Time Dependent 1D Heat Equation, Finite Difference, Explicit Time Stepping

A reference on a finite difference solution of the time dependent 1D heat equation using explicit time stepping in MATLAB (also available are codes in C, C++, and Fortran 77 and 90. There are also other links to similar programs and methods. I try to repost many things as they do not last on the internet forever and this will prove to be an invaluable collection.

FD1D_HEAT_EXPLICIT - TIme Dependent 1D Heat Equation, Finite Difference, Explicit Time Stepping

FD1D_HEAT_EXPLICIT is a MATLAB program which solves the time-dependent 1D heat equation, using the finite difference method in space, and an explicit version of the method of lines to handle integration in time.

FD1D_HEAT_EXPLICIT is available in a C version, a C++ version, a FORTRAN77 version, a FORTRAN90 version and a MATLAB version.

Related Data and Programs:

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 is a MATLAB program which applies the finite difference method to a two point boundary value problem in one spatial dimension.

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

FD1D_HEAT_STEADY is a MATLAB program which uses the finite difference method to solve the steady (time independent) heat equation in 1D.

FD1D_PREDATOR_PREY is a MATLAB program which uses finite differences to solve a 1D predator prey problem.

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 is MATLAB program which applies the finite element method to solve the 2D heat equation.

FEM1D is a MATLAB program which applies the finite element method, with piecewise linear basis functions, to a linear two point boundary value problem;

FEM1D_ADAPTIVE is a MATLAB program which applies the finite element method to a linear two point boundary value problem in a 1D region, using adaptive refinement to improve the solution.

FEM1D_NONLINEAR is a MATLAB program which applies the finite element method to a nonlinear two point boundary value problem in a 1D region.

FEM1D_PMETHOD is a MATLAB program which applies the p-method version of the finite element method to a linear two point boundary value problem in a 1D region.

FEM2D_HEAT is a MATLAB program which applies the finite element method to solve the 2D heat equation.

FREE_FEM_HEAT is a MATLAB program which applies the finite element method to solve the time dependent heat equation in an arbitrary triangulated 2D region.

HOT_PIPE is a MATLAB program which uses FEM_50_HEAT to solve a heat problem in a pipe.

HOT_POINT is a MATLAB program which uses FEM_50_HEAT to solve a heat problem with a point source.

Reference:

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_explicit.m, the MATLAB function which carries out the calculation;

function [ H, x, t ] = fd1d_heat_explicit ( xn, tn )

%*****************************************************************************80
%
%% FD1D_HEAT_EXPLICIT: Finite difference solution of 1D heat equation.
%
% Discussion:
%
% This program solves
%
% dUdT - k * d2UdX2 = F(X,T)
%
% over the interval [A,B] with boundary conditions
%
% U(A,T) = UA(T),
% U(B,T) = UB(T),
%
% over the time interval [T0,T1] with initial conditions
%
% U(X,T0) = U0(X)
%
% The code uses the finite difference method to approximate the
% second derivative in space, and an explicit forward Euler approximation
% to the first derivative in time.
%
% The finite difference form can be written as
%
% U(X,T+dt) - U(X,T) ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) )
% ------------------ = F(X,T) + k * ------------------------------------
% dt dx * dx
%
% or, assuming we have solved for all values of U at time T, we have
%
% U(X,T+dt) = U(X,T) + cfl * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) + dt * F(X,T)
%
% Here "cfl" is the Courant-Friedrichs-Loewy coefficient:
%
% cfl = k * dt / dx / dx
%
% In order for accurate results to be computed by this explicit method,
% the cfl coefficient must be less than 0.5!
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 April 2010
%
% Author:
%
% John Burkardt
%
% Reference:
%
% George Lindfield, John Penny,
% Numerical Methods Using MATLAB,
% Second Edition,
% Prentice Hall, 1999,
% ISBN: 0-13-012641-1,
% LC: QA297.P45.
%
% Parameters:
%
% Input, integer XN, the number of points to use in the spatial dimension.
%
% Input, integer TN, the number of equally spaced points in the time dimension.
%
% Output, real H(XN,TN), the computed solution.
%
% Output, real X(XN), the spatial points.
%
% Output, real T(TN), the time points.
%
timestamp ( );
fprintf ( 1, '\n' );
fprintf ( 1, 'FD1D_HEAT_EXPLICIT:\n' );
fprintf ( 1, ' MATLAB version.\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' Compute an approximate solution to the time-dependent\n' );
fprintf ( 1, ' one dimensional heat equation:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' dH/dt - K * d2H/dx2 = f(x,t)\n' );
%
% Heat coefficient.
%
k = 0.002;
%
% X values;
%
xmin = 0.0;
xmax = 1.0;
x = linspace ( xmin, xmax, xn );
dx = ( xmax - xmin ) / ( xn - 1 );
%
% T values;
%
tmin = 0.0;
tmax = 80.0;
dt = ( tmax - tmin ) / ( tn - 1 );
t = linspace ( tmin, tmax, tn );
%
% Check the CFL condition, have processor 0 print out its value,
% and quit if it is too large.
%
cfl = k * dt / dx / dx;

fprintf ( 1, '\n' );
fprintf ( 1, ' CFL stability criterion value = %f\n', cfl );

if ( 0.5 <= cfl ) fprintf ( 1, '\n' ); fprintf ( 1, 'FD1D_HEAT_EXPLICIT - Fatal error!\n' ); fprintf ( 1, ' CFL condition failed.\n' ); fprintf ( 1, ' 0.5 <= K * dT / dX / dX = %f\n', cfl ); error ( 'FD1D_HEAT_EXPLICIT' ); end % % Compute and save initial values. % h = initial_condition ( xn, x, t(1) ); h = boundary_conditions ( xn, x, t(1), h ); H = zeros ( tn, xn ); H(1,1:xn) = h(1:xn); % % Compute the values of H at the next time, based on current data. % L = 1 : xn - 2; C = 2 : xn - 1; R = 3 : xn; for i = 1 : tn - 1 h(C) = h(C) + cfl * ( h(L) - 2.0 * h(C) + h(R) ) + dt * rhs ( x(C), t(i) ); h = boundary_conditions ( xn, x, t(i+1), h ); H(i+1,1:xn) = h(1:xn); end % % Terminate. % fprintf ( 1, '\n' ); fprintf ( 1, 'FD1D_HEAT_EXPLICIT:\n' ); fprintf ( 1, ' Normal end of execution.\n' ); fprintf ( 1, '\n' ); timestamp ( ); return end function h = boundary_conditions ( xn, x, t, h ) %*****************************************************************************80 % %% BOUNDARY_CONDITIONS evaluates the boundary conditions. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 02 April 2010 % % Author: % % John Burkardt % % Parameters: % % Input, integer XN, the number of nodes. % % Input, real X(XN), the node coordinates. % % Input, real T, the current time. % % Input, real H(XN), the current heat values. % % Output, real H(XN), the current heat values, after boundary % conditions have been imposed. % h(1) = 90.0; h(xn) = 70.0; return end function h = initial_condition ( xn, x, t ) %*****************************************************************************80 % %% INITIAL_CONDITION evaluates the initial condition. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 05 April 2010 % % Author: % % John Burkardt % % Parameters: % % Input, integer XN, the number of nodes. % % Input, real X(XN), the node coordinates. % % Input, real T, the initial time. % % Output, real H(XN), the heat values at the initial time. % h(1:xn) = 50.0; return end function value = rhs ( x, t ) %*****************************************************************************80 % %% RHS evaluates the right hand side of the differential equation. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 05 April 2010 % % Author: % % John Burkardt % % Parameters: % % Input, real X(*), the node coordinates. % % Input, real T, the current time. % % Output, real VALUE(*), the source term. % value = x; value = 0.0; return end function timestamp ( ) %*****************************************************************************80 % %% TIMESTAMP prints the current YMDHMS date as a timestamp. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 14 February 2003 % % Author: % % John Burkardt % t = now; c = datevec ( t ); s = datestr ( c, 0 ); fprintf ( 1, '%s\n', s ); return end

run_and_plot.m, a MATLAB script that supplies input to the function and shows how to plot the results.

%% RUN_AND_PLOT runs the code and plots the results.
%
% Discussion:
%
% This script suggests how to call the code that solves the
% heat equation, and uses some data that is reasonable.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'RUN_AND_PLOT:\n' );
fprintf ( 1, ' Run FD1D_HEAT_EXPLICIT.\n' );
fprintf ( 1, ' Plot the results.\n' );
%
% XN is the number of equally spaced nodes to use between 0 and 1.
%
xn = 21;
%
% TN is the number of equally spaced time points between 0 and 10.0.
%
tn = 201;
%
% Running the code produces an array H of temperatures H(t,x),
% and vectors x and t.
%
[ H, x, t ] = fd1d_heat_explicit ( xn, tn );
%
% Now make a plot of the data.
%
[ X, T ] = meshgrid ( x, t );

mesh ( X, T, H );
title ( 'H(X,T) by FD1D\_HEAT\_EXPLICIT' );
xlabel ( '-- X --' );
ylabel ( '-- Time --' );
zlabel ( '-- H(X,T) --' );

fprintf ( 1, '\n' );
fprintf ( 1, 'RUN_AND_PLOT:\n' );
fprintf ( 1, ' Normal end of execution.\n' );

plot.png, a PNG image of the solution, using the MESH command to emphasize the method of lines approach underlying the solution.