Numerical Math - The Jacobi method

A simple explanation giving by Landau et al [1] defines the Jacobi method as an initial sweep of a values without updating if available. The Jacobi method is very basic. For example, from the square wire problem in Landau et al the numerical algorithm to be solved is

\[ U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \]

The initialization and BCs symmetry are preserved in this way.

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

Numerical Math - The Gauss-Seidel (GS) method

The Gauss-Seidel method is an improvement upon the basic Jacobi method [1]. The GS method utilizes known values and current updates them in the algorithm as opposed to the Jacobi method which sweeps across the domain without utilizing known values to accelerate convergence. The accelerated convergence also produces less round-off error and utilizes less memory since only one generation of guesses is needed.

If the sweep begins at the top left corner then the GS algorithm looks like

\begin{equation} U_{i, j} = \dfrac{1}{4} \left[ U^{(old)}_{i + 1, j} + U^{(new)}_{i - 1, j} + U^{(old)}_{i, j + 1} + U^{(new)}_{i, j - 1} \right] \end{equation}

Compare to the Jacobi method.  Both algorithms shown are for the square wire finite difference technique problem in Landau et al [1].

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \end{equation}

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

Numerical Math - The Finite Element Method (FEM)

From Landau [1], the finite element method (FEM) is explained as solving a PDE where the whole region or domain is subsectioned into smaller areas known as elements. Next an initial solution is formulated for the PDE in each of the elements.  The next step comes in the form of modifying the parameters of the initial formalization known as a best fit.

Pepper and Heinrich [2] open their FEM book by explaining that the method is a numerical formulation that solves physics and engineering models described by differential equations  (DEs). Similar to finite difference methods (FDM), FEM models a geometric region which is then broken up into a number of smaller subregions which results in a network known as a mesh. One difference between FDM and FEM includes mesh type (geometrically). The FDM requires that the network be made of orthogonal rows and columns (squares and rectangles) while the FEM does not require this limitation and can be, in fact, any shape such as triangles and/or quadrilaterals in two dimensions and tetrahedrons and/or hexahedrons in three dimensions.  Next, each FE is initialized with an approximate functions of the unknown variables for which to be solved. Expansions determine the variables approximations and appear as a linear or higher-order polynomials functions. These expansions, in turn, depend upon the geometric shape of the elements and location known as nodes. A second difference between FDM and FEM noted by Pepper and Heinrich [3] entails the solution method. The FEM integrates over each subregion then adds or connects them to make up the whole. This integration and summation engenders a set of finite linear equations per section which can then be solved utilizing linear algebra methods. Jiang [6] also discusses that FEM does not...

The network of elements and nodes or where these elements connect make up discrete systems such as a trusses, circuits, and fluid transport pipes [5]. In order to solve for the system variables such as displacements, electric potentials, and pressures, one can begin with a known and simple parameter such as

Hooke’s Law:

\begin{equation} F = \dfrac{\Delta L}{R} = E A \dfrac{\Delta L}{L} \end{equation}

Ohm’s Law:

\begin{equation} i = \dfrac{\Delta V}{R} = \dfrac{A}{\rho} \dfrac{\Delta V}{L} \end{equation}

Poiseuille’s Law:

\begin{equation} \dot{m} = \dfrac{\Delta p}{R} = \dfrac{\rho \pi D^4}{128 \mu} \dfrac{\Delta p}{L} \end{equation}

Jiang [6] notes that FEM has been utilized as one of the most general numerical techniques to solve DEs and done so with great success. Jiang even quotes Oden as

“Perhaps no other family of approximation methods has had a greater impact on the theory and application of numerical methods during the twentieth century”

In progress...to be continued.

---

References:

[1] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. The Finite Element Method - It’s Basis and Fundamentals, 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

[2] D. W. Pepper and J. C. Heinrich. The Finite Element Method: Basic Concepts and Applications, Taylor & Francis Hemisphere Publishing Corporation,. 1992

[4] J. C. Heinrich and D. W. Pepper. The Intermediate Finite Element Method: Fluid Flow and Heat Transfer Applications, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1999

[5] G. Comini, S. D. Giudice and C. Nonino. Finite Element Analysis in Heat Transfer: Basic Formulation & Linear Problems, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1994

[6] B.-N. Jiang. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Springer-Verlag, Berlin, Germany. 1998

[7] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. The Finite Element Method for Fluid Dynamics, (Volume 3) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

Numerical Math - Square wire problem - finite difference method

The square-wire problem is described by the 2-D Laplace and Poisson equations [1].

\begin{equation} \dfrac{\partial^2 U \left( x, y \right)}{\partial x^2} + \dfrac{\partial^2 U \left( x, y \right)}{\partial y^2} = \left\{ \begin{array}{cr} 0 & \text{Laplace's Equation}\\ -4 \pi \rho \left( \mathbf{x} \right) & \text{Poisson's Equation} \end{array} \right. \end{equation}

Next, in order to obtain the second order derivatives a Taylor series expansion with steps of positive and negative values in both \( x \) and \( y \) are executed for \( U \).

\begin{equation} U \left( x + \Delta x, y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x - \Delta x, y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 - \cdots, \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x, y - \Delta y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 - \cdots, \end{equation}

Now, adding both equations each for \( x \) and \( y \) the odd terms cancel, and truncating at \( O \left( \Delta^4 \right) \) obtains

\begin{equation} U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + O \left( \Delta x^4 \right) \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + O \left( \Delta y^4 \right) \end{equation}

Rearranging to segregate the second derivative engenders

\begin{equation}  \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 = U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right) \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 = U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right) \end{equation}

Finally, dividing produces

\begin{equation} \dfrac{\partial^2 U}{\partial x^2} = \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} = \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} \end{equation}

Next, substituting the second order derivatives into the Poisson equation produces

\begin{equation} \begin{array}{l} \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \\ \hspace{1.15 in} + \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} = -4 \pi \rho \left( x, y \right) \end{array} \end{equation}

If the grid spacing is the same in both the \( x \) and \( y \) coordinates, that is \( \Delta x = \Delta y \), then the equation reduces to

\begin{equation} \begin{array}{l} U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) + U \left( x, y + \Delta \right) \\ \hspace{2 in} + \hspace{1 mm} U \left( x, y - \Delta \right) - 4 U \left( x, y \right) = -4 \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Algebraically, the solution \( U \left( x, y \right) \) can be shown as

\begin{equation} \begin{array}{l} U \left( x, y \right) \simeq \dfrac{1}{4} \left[ U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) \right. \\ \hspace{2 in} \left. + \hspace{1 mm} U \left( x, y + \Delta \right) + U \left( x, y - \Delta \right) \right] + \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Lattice spacings can be identified as \( x = x_0 \pm \Delta \) and \( y = y_0 \pm \Delta \) which correspond to \(i \pm 1 \) and \( j \pm 1 \) where \( i, j = 0, 1, \ldots, N_{max - 1} \).  Thus, the finite difference algorithm appears as

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) + \pi \rho_{i, j}\Delta^2 \end{equation}

However, for simplicity Landau et al just ask to solve for the homogeneous equation

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \end{equation}

A few notes here from Landau et al.  For one a direct solution cannot be obtain from this algorithm.  That is, a process must be repeated several times in order for a solution to converge!  Thus, this is where efficiency and proficiency becomes very important in numerical methods.

First, a guess is imposed then swept over the domain which finds values for all variables at each point or node.  The solution has "converged" when the current generation does not change compared to the previous generation of values within a "tolerance" or error/precision (absolute value of the difference between old values and new) specified by the user or technique.  Another outcome could be that the values emerge as not expected which at that point the user must go back and figure out from where the error is coming.  Landau et al give the definition of relaxed as when the first guess has converged into a solution.

Landau et al pose "Does it always converge, and if so, does it converge fast enough to be useful?"

Another important note is the coordinate system used.  This problem was done in rectangular Cartesian coordinates.  Other problems may work better in other coordinate systems such as cylindrical or spherical depending upon the geometry of the problem in which case the Laplacian expansion \( \Delta^2 \) will be different.

One way to solve the algorithm is to utilize the Jacobi method.  This most basic method simply evaluates each point, sweeping across the domain.  The initialization and BCs symmetry are preserved in this way.  However, in order to speed up this method, a simple modification can be implemented where values that are known from the current run compared to other nodes can be updated immediately.  Thus, reducing the amount of needed storage by half since two nodes will be know for one in this case and the previous sweep will not be necessary to store since information is updated instantaneously when available.  This technique is called the Gauss-Seidel method.  The GS method leads to quicker convergence, less round-off error, and less needed storage.  However, symmetry is not conserved which the user hopes does not affected the solution outcome.

If the sweep begins at the top left corner then the GS algorithm looks like

\begin{equation} U_{i, j} = \dfrac{1}{4} \left[ U^{(old)}_{i + 1, j} + U^{(new)}_{i - 1, j} + U^{(old)}_{i, j + 1} + U^{(new)}_{i, j - 1} \right] \end{equation}

Octave code:

%This code was written by Tim Barber for the numerical PDE class taught by Dr. Parigger in the Summer of 2011

%This file corresponds specifically to HW set #1

%Problem 1.3 Square Wire Problem: Chapter 17 of the Landau-course book

discusses the square-wire problem, viz. elliptic PDE. Section 17.4.2

elaborates on the solution using the finite difference method. Your task

here is to (i) derive the finite difference algorithm Eq. 17.27, (ii)

implement the code using source programing (see the listed source

code), and (iii) display the results using gnuplot (along the lines of

Section 17.4.2).

clf

clc

clear

%input data manually

N_max = 100;

iter = 1000;

%initialize matrix

for i = 1:N_max

for j = 1:N_max

if j = 100

U(i, 100) = 99;

else

U(i, j) = 0;

j = j + 1;

end

end

i = i +1;

end

%iterate matrix using the Gauss-Seidel technique

for n = 1:iter + 1

for i = 2:N_max - 1

for j = 2:N_max - 1

U(i, j) = (1/4)*(U(i + 1, j) + U(i - 1, j)

+ U(i, j + 1) + U(i, j - 1));

j = j + 1;

end

i = i +1;

end

n = n + 1;

end

%plot of U(x, y)

mesh (U)

%label title and axes

title('HW 1.3', "fontsize", 20)

ylabel('y', "fontsize", 16)

xlabel('x', "fontsize", 16)

zlabel('U(x, y)', "fontsize", 16)

%commands to save plot as a .png file

print("HW_1_3_mesh.png")

replot

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics -

a Introductory Computational Science, Princeton University Press, Princeton, New Jersey.

2008

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.