http://interval.louisiana.edu/Classical-and-Modern-NA/index.html

# Classical and Modern Numerical Analysis

## Theory, Methods, and Practice

### Azmy S. Ackleh, Edward J. Allen, R. Baker Kearfott, and Padmanabhan Seshaiyer

Instructors:  I am continuing to correct small errors in the booklet with answers to the exercises.  Please contact me (Ralph Baker Kearfott) if you wish to have an updated copy.

We list here, by chapter and context, Matlab codes that support topics in our book.  We will continue to add to this collection.  To the extent practical, we have put these programs in a uniform style, and have provided in-line help.  For example, to use bisect.m, put bisect.m in a directory within Matlab's search path;  you can then see an explanation by typing "help bisect" from Matlab's command window.

Functions associated with interval arithmetic may be obtained from the web site for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, that is, at:  http://www.siam.org/books/ot110.  These Matlab functions require INTLAB, a high-quality Matlab toolbox that is available free of charge for non-commercial use.  INTLAB has many demos integrated with Matlab's help system.  Jiri Rohn has provided a short INTLAB primer.

Another resource for Matlab-related illustrative functions is the on-line book Numerical Computing with Matlab (SIAM, Philadelphia, 2004) and associated set of  Matlab "m" files, available at  http://www.mathworks.com/moler/index_ncm.html.

/ Chapter 1 / Chapter 2 / Chapter 3 / Chapter 4 / Chapter 5Chapter 6 / Chapter 7 / Chapter 8 / Chapter 9 / Chapter 10 /

#### Chapter 1

difference_table.m
This Matlab function (also introduced in Chapter 6) can be used to produce a table like that on page 9 of the text, illustrating roundoff error.
mean_value_form.m
Function that computes the mean value form of an interval extension of a function, as explained in problem 26, page 33 of our text.  This function is available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.

#### Chapter 2

bisect_method.m
Function for the univariate method of bisection, in Section 2.2, and possibly in problems 2, 3, and 4 of Chapter 2 (Section 2.10).  We also include the sample function xsqm2.m and  the Matlab script run_bisect_method.m.
steffensen.m
Function for Steffensen's method, in Section 2.7.2 and problems 16 and 34 of Chapter 2 (Section 2.10).  We also include the script run_steffensen.m.
newton.m
Function for the univariate Newton method, as described in Section 2.4 and in  problems 20, 21, and (with modification) problem 31 of Chapter 2 (Section 2.10).
i_newton_step.m and i_newton_step_no_fp.m
Functions that perform a step of  the univariate interval Newton method, as explained in Section 2.5, starting on page 54 of our text, and may be useful in  problem 28 of Chapter 2 (in Section 2.10).  i_newton_step requires programming not only the function f, but also the derivative of f, whereas i_newton_step_no_fp uses INTLAB's automatic differentiation toolbox, and thus does not require programming of the derivative of f.  These functions are available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.
Krawczyk_step.m
Function that performs a step of the Krawczyk method, as explained in problem 29 of Chapter 2 (page 81, Section 2.10) of our text.  This function can be used both for the univariate and multivariate Krawczyk methods.  This function is available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.

#### Chapter 3

Numerical linear algebra lies at the foundation of the Matlab system, so many of the concepts in this chapter on numerical linear algebra are available directly as Matlab functions.  Browse and search within Matlab "help" and explore the demos.  (From the Matlab command window, select "help" on the menu bar at the top, then select "demos".)

For  more detailed explanations, as well as additional examples and support for elementary interval linear algebra routines for concepts in our text, see the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009.  We list some Matlab / INTLAB functions from the collection below.

gauss_seidel_step.m
Function that performs a  step of the interval Gauss--Seidel method for a single specified coordinate.  This function is available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.  It is called by the function Gauss_Seidel_image.m, also available with the collection on http://www.siam.org/books/ot110.  The interval Gauss--Seidel method is covered in Section 3.4.5, page 153, of our text.

#### Chapter 4

Matlab has much built-in support for interpolation, using state-of-the art implementations.  This includes splines, Fourier transforms, etc.  Browse and search within Matlab's help facility.  We have supplied the following functions, meant to illustrate and serve as a tutorial for basic concepts in the text.

plot_Runge_and_Lagrange_interpolant.m
Matlab script that first computes the coefficients of the degree 10 Lagrange interpolating polynomial (in terms of powers of x) with equally spaced points over the interval -5 to 5 to Runge's function, then plots both the values of this interpolating polynomial and the values of Runge's function with 200 equally spaced points on the interval [-5,5].  This script calls the function Lagrange_interp_poly_coeffs.m and Lagrange_interp_poly_val.m.  The script and Lagrange_interp_poly_coeffs.m use the supplied function runge.m.   Lagrange_interp_poly_coeffs.m and Lagrange_interp_poly_val.m are general, and can be used for Lagrange interpolation of any degree for any function over any interval.  The Lagrange interpolating polynomial is explained in Section 4.3.3, starting on page 209 of the text.  Runge's function is introduced in problem 5 of Chapter 4 (Section 4.9) as an example of a function for which interpolation with equally spaced points goes awry.  The m-files listed here can be used to do that problem.
piecewise linear interpolation
We cover piecewise linear interpolation in Section 4.4.1 of the text (starting on page 238). Students can try piecewise linear interpolation with the Matlab function interp1, with the call
yi = interp1(x,y,xi,'linear'), where x is the vector of abscissas, y is the vector of ordinates, xi is the vector of points at which values should be computed, and the returned vector yi contains the values of the linear interpolating function at corresponding components of xi.  The m-file for interp1 can be found at present in the Matlab subdirectory
[Matlab root directory]/toolbox/matlab/polyfun/interp1.m, and can be examined by students, copied, and modified.  In addition, we have supplied the function plot_piecewise_linear_interpolant.m, that are somewhat simpler than interp1.m, and that students can modify and study.  The function plot_piecewise_linear_interpolant.m computes a piecewise linear interpolant at m points, then plots it;  it uses either a straightforward computation or hat function basis.  The straightforward computation is done with either piecewise_linear_interpolant_value.m or with interpolant_with_hat_function_basis.m, depending on the user's choice.  The routine interpolant_with_hat_function_basis.m uses the auxiliary routine hat_function_value.m.
cubic spline interpolation
This can be done with an option in the Matlab function interp1 or directly through the Matlab function spline.  Spline interpolation is analyzed in Section 4.4.2 of the text, and is featured in problems 25 and 26 of Chapter 4 (Section 4.9).
piecewise cubic Hermite interpolation
This can also be done with an option in the Matlab function interp1.  Hermite interpolation is discussed in Section 4.3.5 of the text, and piecewise cubic Hermite interpolation is featured in Problem 19 of Chapter 4 (Section 4.9).
plot_trigonometric_interpolant.m
Function that plots a trigonometric interpolant to a function over [0,2π].  It calls trigonometric_interpolant_value.m.  We also supply the example function trig_interp_example.m.

#### Chapter 5

As with the linear algebra concepts in Chapter 3, eigenvalue / eigenvector computations lie at the foundation of the Matlab system.  Browse and search within Matlab's help facility.  For example, the Matlab function eig  computes eigenvalues and eigenvectors with the QZ algorithm, a form of the QR method, where the QR method is treated  in Section 5.5 of the text.  As with most Matlab functions, the Matlab source code is distributed with Matlab, in this case in the file
[Matlab root directory]/toolbox/eml/lib/matlab/matfun/eig.m (as of Version R2008a).

We supply the following simple routines for illustrative and didactic purposes.
inverse_power_method.m
Function that computes an approximation to an eigenvalue and corresponding eigenvector using the inverse power method as described in Section 5.3 of the text.  We also supply a simple Matlab script run_inverse_power_method.m that tries inverse_power_method.m with some simple examples.

#### Chapter 6

difference_table.m
Function that produces a table of finite differences, using the function difference_quotient.m.  The function difference_quotient.m programs the forward difference as in formula (6.1) (page 323) of the text.  A table similar to that of Example 6.2 (page 326) may be produced with this function.  Also, with modification of difference_quotient.m, problem 3 of Chapter 6 (Section 6.4) may be done.
composite_Newton_Cotes.m
Function for composite Newton-Cotes quadrature, as explained in Section 6.3.1 and Section 6.3.2 of the text;  type "help composite_Newton_Cotes" from the Matlab command window after downloading.  We also supply the example function f_sin_x_over_xsqp1.m and the example Matlab script run_composite_Newton_Cotes.m.
n_d_Monte_Carlo_integral.m
Function that computes an approximation to an n-dimensional integral using Monte Carlo integration, as explained in Section 6.3.6 of the text.  We also include the example function sum_of_squares.m.  Use of n_d_Monte_Carlo_integral.m is explained in the in-line documentation (in the m-file itself).

#### Chapter 7

Matlab provides a suite of  programs for the numerical solution of systems of ordinary differential equations,  of various orders (as explained in Sections 7.3 and 7.4) of the text, and for both stiff and non-stiff systems (where stiff systems are explained in Section 7.7 of the text.)  An example of this is
[Matlab root directory]/toolbox/eml/lib/matlab/ode45.m
which implements a classic Runge-Kutta fourth order / fifth order pair of methods.  (Here,  the different orders are used to estimate the error in a step.  The fourth-order method is given on page 398 of the text.)  Further examples can be obtained by entering "differential equations" into the search box in Matlab's help window.  For professional solution to initial value problems for systems of ordinary differential equations, we recommend using these functions provided with Matlab.

For illustration and teaching, we provide the following.
Runge_Kutta_2_for_systems.m
Function that implements the modified Euler method (a second-order Runge--Kutta method) as explained at the bottom of page 397 of the text.  We also include Lotka_Volterra.m as an example system for  Runge_Kutta_2_for_systems to integrate, and we include the script  run_Runge_Kutta_2_for_systems.m to run the Lotka--Volterra model and plot the results.
Function that implements an order 3 Adams Bashforth multistep method as explained in Section 7.5 of the text.  We also include Lotka_Volterra.m as an example system  to integrate, and we include the script  run_Adams_Bashforth_order_3.m to run the Lotka--Volterra model and plot the results.  (The script run_Adams_Bashforth_order_3.m differs little from run_Runge_Kutta_2_for_systems.m.)

#### Chapter 8

For  more detailed explanations, as well as additional examples and support for the interval Newton methods described below, see the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009.
newton_sys.m
Function that implements the multi-dimensional Newton method, as explained in Section 8.3 of the text (and, especially, Algorithm 8.1 on page 458, section 8.3.7).  (Only the range tolerance is used, and not the domain tolerance.)  We also include the functions simple_2d.m and simple_2d_prime.m, that can be used for Problem 9 on page 483 (Section 8.7).  The routine is a slight modification of newton.m, used for Chapter 2;  studying newton.m and newton_sys.m may be helpful in understanding Matlab syntax.  Finally, we include the Matlab diary newton_sys_example.txt, to illustrate how newton_sys.m can be used.
nonlinear_Gauss_Seidel_image.m
Function that performs a sweep of the nonlinear interval Gauss--Seidel method, as explained in Section 8.4.1, starting on page 465 of our text.  This function is available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.
Krawczyk_step.m
Function that performs an iteration of the nonlinear Krawczyk method, as explained in Section 8.4.2, starting on page 470 of our text. This function is available with the collection for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009, at the web page  http://www.siam.org/books/ot110.
homotopy_step.m
Function that implements the homotopy method for finding the zeros of a function g, given a related function f whose zeros are known, according to formulas (8.92) through (8.96) in Section 8.6.1 (on page 480).  We also include the functions f_prob_8_27.mg_prob_8_27.m, fp_prob_8_27.m, and gp_prob_8_27.m, that can be used to do problem 27 of Chapter 8 (Section 8.7).  Finally, we include the Matlab diary problem_8_27_example.txt, illustrating use of homotopy_step.m to do problem 8(a) in Section 8.7.

#### Chapter 9

fminsearch.m
Matlab function that implements the Nelder--Mead simplex method, as described in Section 9.1.6.   The Nelder--Mead method is also featured in problem 4 on page 531, although the computations in that problem should be done by hand, to understand the inner workings of the method.  The file fminsearch.m can be found in
[Matlab root directory]/toolbox/matlab/fminsearch.m
linprog.m
Function to solve linear programming problems.  This function is available if you have Matlab's optimization toolbox.  If you do have that toolbox, linprog can be found in
[Matlab root directory]/toolbox/optim/linprog.m.
The function linprog is also documented well within Matlab's help system, while the in-line documentation for it can be displayed by typing "help linprog" from within Matlab's command window.  linprog is relevant to Section 9.4 of the text, and can be useful for problems 9, 11, and possibly problem 14 (if the instructor does not require it be done by hand), in Chapter 9 (Section 9.7).
Skelboe-Moore.m
A Matlab / INTLAB function implementing a simple branch and bound method for global optimization (Section 9.6.3 of the text),  available with the collected set of m-files at  http://www.siam.org/books/ot110  for the book Introduction to Interval Analysis, by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud, SIAM, Philadelphia, 2009.

#### Chapter 10

Matlab has functions to solve first-order boundary value problems, that is, for differential equations of the form  dy/dx = f(x,y);  see the Matlab help for the function bvp4c.  Matlab also has support for various types of higher-order boundary value problems for ordinary and partial differential equations.  Check the Matlab demo "Differential Equations in Matlab."  For teaching and illustration purposes, we provide the following.
two_point_BVP_with_finite_differences.m
Function that returns an approximate solution to the linear two-point boundary value problem d2u/dx2 = f(x),  u(a) = ua,  u(b) = ub, using a central difference approximation to the second derivative.  We also supply the Matlab script run_two_point_BVP_with_finite_differences.m to run two_point_BVP_with_finite_differences to study convergence of the central difference approximations., as well as the function  minus_32_x_cubed.m as an example function f.  Such finite difference methods are explained in Section 10.1.2, starting on page 540 of the text.
Fredholm_second_kind.m
Function that provides an approximate solution to the Fredholm integral equation of the second kind as described in Section 10.2.2.2 (Formula (10.73), page 562) of our text.  It uses the Trapezoidal rule for the quadrature rule, but has commented lines which, when uncommented, will implement Simpson's rule.  We provide a driver script run_Fredholm_second_kind.m that runs Fredholm_second_kind for various mesh sizes on an example problem,  to study the convergence rate.  We also provide the example functions example_Fredholm_kernel.m and example_Fredholm_g.m.