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.

## No comments:

## Post a Comment