Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

CONTENTS

Chapter 4.   Interpolation and Polynomial Approximation
	Algorithm 4.1  Evaluation of a Taylor Series.
	Algorithm 4.2  Polynomial Calculus.
	Algorithm 4.3  Lagrange Approximation.
	Algorithm 4.4  Nested Multiplication with Multiple Centers.
	Algorithm 4.5  Newton Interpolation Polynomial.
	Algorithm 4.6  Chebyshev Approximation.

Clear["Global`*"];
Off[General::"spell1"];
	Clear[a, A, A0, c, c0, Chebyshev, Cpolys, DiffPoly, f, f0, Horner, 
        IntPoly, j, k, Lagrange, n, NestMult, NewtonP, p, prod, P, P0, s0, 
        s1, sum, SumTaylor, t, term, tol, TaylorCoeff, x, x0, X, X0, Y, XY, 
        XY0];
On[General::"spell1"];

Algorithm 4.1  Evaluation of a Taylor Series.   To approximate
P (x) = (∑_(k = 0)^∞)((Derivative[k][f]((Subscript[x, 0]))(((x - (Subscript[x, 0])))^k))/(k!)) = (∑_(k = 0)^∞)((D (k) (((x - (Subscript[x, 0])))^k))/(k!)),
by computing a partial sum with at most  n  terms. 
Section 4.1 Taylor Series & Calculation of Functions  Page 203

Only a finite number of terms can be accumulated. Hence, the summation process is terminated when consecutive partial sums differ by less than the preassigned value  Tol .

Execute this cell to activate Evaluation of a Taylor Series.
TaylorCoeff[f0_,x0_,n_]
SumTaylor[c0_,x_,x0_,tol_]
f0_ is the function f(x),
x0_ is the value for the center,
n_ is the degree of the polynomial,
c0_ is the coefficient list,
x_ is the variable for evaluating the polynomial,
tol_ is the tolerance.

Off[General::"spell1"]; Clear[TaylorCoeff, f0, x0, n]; \(\(TaylorCoeff[f0_, x0_, n_] := 
	
      Module[{f, s0, s1, X}, 
		Set@@{f[x_], f0}; 
\ \ \ \ \ \ 
        s0 = Normal[Series[f[x], {x, x0, n}]]; 
\ \ \ \ \ \ 
        s1 = ReplaceAll[s0, x -> X + x0]; 
\ \ \ \ \ \ 
        c = CoefficientList[s1, X]; 
\ \ \ \ \ \ 
        Print["Coefficient list for the function"]; 
\ \ \ \ \ \ 
        Print[f[x], "  expanded about x = ", x0]; 
\ \ \ \ \ \ 
        Print[c]; ]; \)\)

Clear[SumTaylor, c0, x, x0, tol]; \(SumTaylor[c0_, x_, x0_, tol_] := 
	
      Module[{k, n, prod, sum}, 
\ \ \ \ \ \ c = c0; 
		n = Length[c]; 

        		sum = 0; 
		prod = 1; 
		
        For[k = 1, k ≤ n, \(k++\), 
			
          \(Module[{absval, term}, 
				prod = prod\ \((x - x0)\); 

            \ \ \ \ \ \ \ \ \ 
            term = c[[k]]\ prod; 

            \ \ \ \ \ \ \ \ \ absval = N[Abs[term]]; 
\ \ \ \ \ \ \ \ \ 
            If[0 < absval && absval ≤ tol, Break[]]; 

            \ \ \ \ \ \ \ \ \ sum = sum + term]; \)]; 
\ \ \ \ \ \ \ \ \ 
        Return[N[sum, 20]]; ]; 
\(On[General::"spell1"]; \)\)

Chapter 4, Example for Taylor polynomials, Page 203-204.

TaylorCoeff[Sin[x], 0, 11]

SumTaylor[c, 1.0, 0, 0.01]

SumTaylor[c, 1.0, 0, 0.001]

SumTaylor[c, 1.0, 0, 0.0001]

SumTaylor[c, 1.0, 0, 0.00001]

The summation process was terminated s terminated when the magnitude of the term to be added was less than  0.001,  0.0001  and  0.0000001  respectively.
     If the goal it to evaluate the Taylor polynomial of a fixed degree, then it is easiest to use Mathematica's built in procedures.  For example, consider the 11-th degree Taylor polynomial for sin(x):

\(f[x_]\  = \ Sin[x]; 
n\  = \ 11; 

    P[x_] = Normal[Series[f[x], {x, 0, 11}]]; 
Print["  f[x] = ", f[x]]; 
    
Print[\*""\!\(P\_11\)[x] = "", P[x]]; 

    Print[\*""\!\(P\_11\)[1] = "", P[1]]; 

    Print[\*""\!\(P\_11\)[1] = "", N[P[1], 20]]; \)

Notice that all the terms in the Taylor polynomial were added in the computation of  P[1].  Compare this with the previous method:

SumTaylor[c, 1, 0, 0.000001]

SumTaylor[c, 1, 0, 0.00000001]

Notice that the magnitude of the last term is smaller than  0.000001  and larger than  0.00000001,  so that it was added in when we requested  SumTaylor[c,1,0,0.00000001], but was not added into the previous sum.

Algorithm 4.2 Polynomial Calculus. To evaluate the polynomial  P (x)   its derivative  P'(x),  and integral  ∫P (x) dx  by performing synthetic division.
Section  4.2   Introduction to Interpolation   Page  212

Execute this cell to activate Polynomial Calculus.
Horner[P0_,t_]
DiffPoly[P0_,t_]
IntPoly[P0_,t_]
P0_ is the polynomial,
t_ is the variable.

Off[General::"spell1"]; Clear[Horner, P0, t, x]; \(\(Horner[P0_, t_] := 
	
      Module[{c, P, k, n, X}, 
		Set@@{P[x_], P0}; 
		\ 
        X = \(Variables[P[x]]\)[[1]]; 

        		\ c = CoefficientList[P[x], X]; 
		\ n = Length[c] - 1; 

        		\ poly = c[[n + 1]]; 
		\ 
        For[k = n - 1, 0 ≤ k, \(k--\), 
			
          poly = c[[k + 1]] + poly\ t]; 
        
		Return[poly]; ]; \)\)

Clear[DiffPoly, P0, t, x]; \(\(DiffPoly[P0_, t_] := 
	
      Block[{c, P, k, n, X}, 
		Set@@{P[x_], P0}; 
		\ 
        X = \(Variables[P[x]]\)[[1]]; 

        		\ c = CoefficientList[P[x], X]; 
		\ n = Length[c] - 1; 

        		\ diff = n\ c[[n + 1]]; 

        		\ For[k = n - 1, 1 ≤ k, \(k--\), 
			\ 
          diff = k\ c[[k + 1]] + 
              diff\ t]; 
		Return[diff]; ]; \)\)

Clear[IntPoly, P0, t, x]; \(IntPoly[P0_, t_] := 
	
      Block[{c, P, k, n, X}, 
		Set@@{P[x_], P0}; 
		\ 
        X = \(Variables[P[x]]\)[[1]]; 

        		\ c = CoefficientList[P[x], X]; 
		\ n = Length[c] - 1; 

        		\ intpol = 
          c[[n + 1]]\/\(n + 1\); 
		\ 
        For[k = n, 1 ≤ k, \(k--\), 
			
          intpol = 
            c[[k]]\/k + intpol\ t]; 

        		intpol = intpol\ t; 
		\ Return[intpol]; ]; 

    \(On[General::"spell1"]; \)\)

This example is for illustration purposes only.  In practice one should modify the procedure to input the coefficient  c.  Put the following polynomial in nested form.

\(Clear[P, Q, x]; 
\(P[x_] = 3\ x\^3 - 2\ x\^2 + 5\ x - 7; 

    Q[x_] = Horner[P[x], x]; 
Print["P[x] = ", P[x]]; 

    Print["Q[x] = ", Q[x]]; \)\)

We can verify that Mathematica has stored the two forms of the polynomial differently by looking at them with FullForm.

FullForm[P[x]]

FullForm[Q[x]]

Evaluate the nested form Q(1) and compare with P(1).

\(Print["P[1] = ", P[1]]; 
Print["Q[1] = ", Q[1]]; \)

Verify that the two forms of the polynomial are equivalent by expanding Q(x)..

\(Print["Q[x] = ", Q[x]]; 
Print["Q[x] = ", Expand[Q[x]]]; 

    Print["P[x] = ", P[x]]; \)

Now consider the derivative of the polynomial P(x).
Put the derivative of the polynomial in nested form.

dP[x_] = DiffPoly[P[x], x];
Print[" P[x] = ", P[x]];
Print["
P'[x] = ",P'[x]];
Print["dP[x] = ", dP[x]];

Evaluate the derivative of P(x) at a specified value of x.

Print["dP[1] = ", dP[1]];
Print[" 
P'[1] = ",P'[1]];

Verify that the two forms of the derivative are equivalent.

Print["dP[x] = ", dP[x]];
Print["dP[x] = ", Expand[dP[x]]];
Print["
P'[x] = ",P'[x]];

Now consider the integral of the polynomial P(x).
Put the integral of the polynomial in nested form.

Clear[V, x]; \(V[x_] = ∫P[x] \[DifferentialD]x; 
iP[x_] = IntPoly[P[x], x]; 
    
Print["    P[x] = ", P[x]]; 

    Print[" ∫P[x]\[DifferentialD]x = ", V[x]]; 

    Print["    V[x] = ", V[x]]; 

    Print[" ∫P[x]\[DifferentialD]x = ", iP[x]]; 

    Print["   iP[x] = ", iP[x]]; \)

Evaluate the integral of P(x) at  x = 1.

\(Print["iP[1] = ", iP[1]]; 
Print[" V[1] = ", V[1]]; \)

Verify that the two forms of the integral are equivalent.

\(Print["iP[x] = ", iP[x]]; 
Print["iP[x] = ", Expand[iP[x]]]; 

    Print[" V[x] = ", V[x]]; \)

Algorithm 4.3  Lagrange Approximation.  To evaluate the Lagrange polynomial
	 P (x) = (∑_(k = 0)^n)(Subscript[y, k]) (Subscript[L, n, k]) (x)  of degree  n,  based on the  n+1  points  ((Subscript[x, k]), (Subscript[y, k]))  for  k = 0, 1, ⋯, n.
Section  4.3   Lagrange Approximation   Page  224

The Lagrange coefficient polynomials for degree  n  are:
(Subscript[L, n, k]) (x)=(((x - (Subscript[x, 0]))⋯(x - (Subscript[x, k - 1]))(x - (Subscript[x, k + 1]))⋯(x - (Subscript[x, n])))/( ((Subscript[x, k]) - (Subscript[x, 0]))⋯((Subscript[x, k]) - (Subscript[x, k - 1]))((Subscript[x, k]) - (Subscript[x, k + 1]))⋯((Subscript[x, k]) - (Subscript[x, n]))))   for  k = 0, 1, ⋯, n.

Execute this cell to activate Lagrange Approximation.
Lagrange[XY0_,t_]
XY0_ is the list of points,
t_ is the variable.

Off[General::"spell1"]; Clear[Lagrange, XY0, t]; \(Lagrange[XY0_, t_] := 
	
      Block[{j, k, n, prod, sum, term, X, Y, XY}, 
		XY = XY0; 
		\ 
        n = Length[XY] - 1; 
		\ 
        X = \(Transpose[XY]\)[[1]]; 

        		\ Y = 
          \(Transpose[XY]\)[[2]]; 

        		\ sum = 0; 
		\ Clear[P, L]; 
		\ 
        For[k = 0, k ≤ n, \(k++\), 
			prod = 1; 
			
          For[j = 0, j ≤ n, \(j++\), 
				
            term = Which[j == k, 1, j ≠ k, 
                \(t - X[[j + 1]]\)\/\(X
                      [[k + 1]] - 
                    X[[j + 1]]\)]; 

            				prod = prod\ term]; 
			L[k, x_] = prod; 
			
          sum = sum + Y[[k + 1]]\ prod]; 
        
		P[x_] = sum; 
		Return[P[x]]; ]; 

    \(On[General::"spell1"]; \)\)

Chapter 4, Example 4.7, Page 215.
Construct two linear Lagrange interpolating polynomials for the function
f(x) = cos(x)  over  [0.0, 1.2]  using different interpolation nodes.

\(Clear[f, P1, Q1, x]; 
f[x_] = Cos[x]; 

    XY = {{0.0, f[0.0]}, {1.2, f[1.2]}}; 
P1[x_] = Lagrange[XY, x]; 

    Q1[x_] = Expand[P1[x]]; 
Print["f[x] = ", f[x]]; 

    Print["Nodes = ", XY]; 
Print[\*""\!\(P\_1\)[x] = "", P1[x]]; 
    
Print[\*""\!\(Q\_1\)[x] = "", Q1[x]]; \)

\(\(\ dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr1 = Plot[{f[x], Q1[x]}, {x, 0.0, 1.2}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr1, dots, PlotRange -> {{0, 1.3}, {0.0, 1.15}}, 
      AxesLabel -> {"x", "y"}, 
      PlotLabel -> \*""f[x] and \!\(P\_1\)[x]"", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_1\)[x] = "", P1[x]]; \)\)

Compare with Mathematica's built in Fit procedure.

\(Print[\*""\!\(R\_1\)[x] = "", Fit[XY, {1, x}, x]]; 

    Print[\*""\!\(Q\_1\)[x] = "", Q1[x]]; \)

Construct the other Lagrange polynomials of degree n=1.

\(Clear[f, P2, Q2, x]; 
f[x_] = Cos[x]; 

    XY = {{0.2, f[0.2]}, {1.0, f[1.0]}}; 
P2[x_] = Lagrange[XY, x]; 

    Q2[x_] = Expand[P2[x]]; 
Print["f[x] = ", f[x]]; 

    Print["Nodes = ", XY]; 
Print[\*""\!\(P\_2\)[x] = "", P2[x]]; 
    
Print[\*""\!\(Q\_2\)[x] = "", Q2[x]]; \)

\(dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr2 = Plot[{f[x], Q2[x]}, {x, 0.0, 1.2}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr2, dots, PlotRange -> {{0, 1.3}, {0.0, 1.15}}, 
      AxesLabel -> {"x", "y"}, 
      PlotLabel -> \*""f[x] and \!\(P\_2\)[x]"", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

Notice that the two polynomials of degree n=1 were different.

\(Print[\*""\!\(P\_1\)[x] = "", P1[x]]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

Chapter 4, Example 4.8, Page 218.   Form several Lagrange polynomials 
of degree  n = 2  and  n = 3  for the function  f(x)  =  cos(x).

First, construct the Lagrange interpolation polynomial of degree n = 2.

\(Clear[f, P2, Q2, x]; 
f[x_] = Cos[x]; 

    XY = {{0.0, f[0.0]}, {0.6, f[0.6]}, {1.2, f[1.2]}}; 

    P2[x_] = Lagrange[XY, x]; 
Q2[x_] = Expand[P2[x]]; 

    Print["f[x] = ", f[x]]; 
Print["Nodes = ", XY]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; 

    Print[\*""\!\(Q\_2\)[x] = "", Q2[x]]; \)

Compare with Mathematica's Fit procedure.

\(Print[\*""\!\(R\_2\)[x] = "", Fit[XY, {1, x, x\^2}, x]]; 

    Print[\*""\!\(Q\_2\)[x] = "", Q2[x]]; \)

dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; Print["The function is f[x] = ", f[x]]; \(gr3 = 
      Plot[{f[x], Q2[x]}, {x, \(-0.6\), 2.1}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr3, dots, PlotRange -> {{\(-0.6\), 2.1}, {\(-0.6\), 1.05}}, 
      AxesLabel -> {"x", "y"}, 
      PlotLabel -> \*""f[x] and \!\(P\_2\)[x]"", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

Second, construct the Lagrange interpolation polynomial of degree n = 3.

\(Clear[f, P3, Q3, x]; 
f[x_] = Cos[x]; 

    XY = {{0.0, f[0.0]}, {0.4, f[0.4]}, {0.8, f[0.8]}, {1.2, f[1.2]}}; 

    P3[x_] = Lagrange[XY, x]; 
Q3[x_] = Expand[P3[x]]; 

    Print["f[x] = ", f[x]]; 
Print["Nodes = ", XY]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_3\)[x] = "", P3[x]]; 

    Print[\*""\!\(Q\_3\)[x] = "", Q3[x]]; \)

Compare with Mathematica's Fit procedure.

\(Print[\*""\!\(R\_3\)[x] = "", Fit[XY, {1, x, x\^2, x\^3}, x]]; 

    Print[\*""\!\(Q\_3\)[x] = "", Q3[x]]; \)

\(dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr4 = Plot[{f[x], Q3[x]}, {x, \(-0.6\), 2.1}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr4, dots, PlotRange -> {{\(-0.6\), 2.1}, {\(-0.6\), 1.05}}, 
      AxesLabel -> {"x", "y"}, 
      PlotLabel -> \*""f[x] and \!\(P\_3\)[x]"", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Lagrange polynomial is:"]; 

    Print[\*""\!\(P\_3\)[x] = "", P3[x]]; \)

Algorithm 4.4 Nested Multiplication with Multiple Centers. 
To evaluate the Newton form of a polynomial
\(P \((x)\)\  = \ 
        a\_0 + \(a\_1\) \((x - x\_0)\) + 
          \(a\_2\) \((x - x\_0)\)\ \((x - x\_1)\) + 
          \(a\_3\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + ⋯ + 
          \(a\_n\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)\ 
            ⋯ \((x - x\_\(n - 1\))\)\)
with multiple centers  (Subscript[x, 0]), (Subscript[x, 1]), (Subscript[x, 2]),  ⋯, (Subscript[x, n - 1]).
Section  4.4   Newton Polynomials   Page  233

Execute this cell to activate Nested Multiplication with Multiple Centers. 
NestMult[A0_,X0_,n_,t_]
A0_ is the list of coefficients,
X0_ is the list of centers,
n_ is the degree,
t_ is the variable.

Off[General::"spell1"]; Clear[NestMult, A0, X0, n, t]; \(NestMult[A0_, X0_, n_, t_] := 
	
      Block[{A, X, k}, 
		A = A0; 
		X = X0; 
        sum = A[[n + 1]]; 
		
        For[k = n, 1 ≤ k, \(k--\), 
			
          sum = sum\ \((t - X[[k]])\) + 
              A[[k]]]; 
		
        Return[sum]; ]; 
\(On[General::"spell1"]; \)\)

What is nested multiplication?  See Page 228.

Clear[A, X, a, x]; A = Table[a[i], {i, 5}]X = Table[x[i], {i, 4}]

NestMult[A, X, 1, x]

NestMult[A, X, 2, x]

NestMult[A, X, 3, x]

Chapter 4, Example 4.11, Page 228.   Using nested multiplication, evaluate
(Subscript[P, 3]) (x)  =  
        ((-0.1 (x - 4) + 0.5) (x  - 3) - 2) 
            (x - 1) + 5.

\(Clear[A, P, x, X]; 
\(A = {5, \(-2\), 0.5, \(-0.1\)}; \)\)\(X = {1, 3, 4}; 
P[x_]\  = \ NestMult[A, X, 3, x]; 

    Print["  P[x] = ", P[x]]; 
Print["P[2.5] = ", P[2.5]]; \)

Algorithm 4.5 Newton Interpolation Polynomial.   To construct and evaluate the 
Newton polynomial of  degree  ≤ n  that passes through  ((Subscript[x, k]), (Subscript[y, k]))  =  ((Subscript[x, k]), f ((Subscript[x, k]))),
for  k = 0, 1, ⋯, n: 
\(P \((x)\)\  = \ 
        d\_\(0, 0\) + \(a\_1\) \((x - x\_0)\) + 
          \(d\_\(1, 1\)\) \((x - x\_0)\)\ \((x - x\_1)\) + 
          \(d\_\(2, 2\)\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + ⋯ + 
          \(d\_\(n, n\)\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)\ 
            ⋯ \((x - x\_\(n - 1\))\)\)
where   dk,0  =  yk    and    dk,j  =  (dk,j-1 - dk-1,j-1)/(xk - xk-j). 
Section  4.4   Newton Polynomials   Page  234

The next loop will construct the divided difference table for f(x) based on the  n  points  
(Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, n + 1]).  Then it will simplify the diagonal coefficients.  Then it will take the 
limit to get derivatives  f'((Subscript[x, 1])).

In Mathematica, the coefficients are stored as "function" values  x[k]. The  following 
version is part of the Tangent Curves paper and it uses equally spaced abscissas.

Remark:  This is not the general Newton polynomial where the abscissas are arbitrary, i.e. not necessarily equally spaced.

Execute this cell to activate Newton Interpolation Polynomial.
NewtonP[XY0_]
XY0_ is the list of points.

Off[General::"spell1"]; \(Clear[NewtonP]; 

    NewtonP[XY0_] := 
	
      Block[{j, k, n, a, p, sum, X, Y}, 
		
        X = \(Transpose[XY0]\)[[1]]; 

        		Y = \(Transpose[XY0]\)[[2]]; 
        
		n = Length[XY0] - 1; 
		d = Table[0, {n + 1}, {n + 1}]; 

        		a = Table[0, {n + 1}]; 
		p = Table[0, {n + 1}]; 
		
        For[k = 0, k ≤ n, \(k++\), 
          \(d[[k + 1, 1]] = 
            Y[[k + 1]]; \)]; 
		
        For[j = 1, j ≤ n, \(j++\), 
			
          \(For[k = j, k ≤ n, \(k++\), 
				
            d[[k + 1, j + 1]] = 
              \(d[[k + 1, j]] - 
                  d[[k, j]]\)\/\(X
                    [[k + 1]] - 
                  X[[k - j + 1]]\)]; 
          \)]; 
		
        For[k = 0, k ≤ n, \(k++\), 
			
          \(a[[k + 1]] = 
            Together[d[[k + 1, k + 1]]]; 
          \)]; 
		p[[1]] = 1; 
		
        For[k = 1, k ≤ n, \(k++\), 
			
          p[[k + 1]] = 
            p[[k]]\ 
              \((x - X[[k]])\)]; 
		
        sum = 0; 
		
        For[k = 0, k ≤ n, \(k++\), 
			
          \(sum = sum + 
              a[[k + 1]]\ 
                p[[k + 1]]; \)]; 
		
        sum = ReplaceAll[sum, h -> H]; 
		Return[sum]; ]; 

    \(On[General::"spell1"]; \)\)

Chapter 4, Exercise 7, Page 235.
Find the Newton polynomial of degree n=4, and evaluate it at  x = 2.5,   
for the data points  (1,  3.6), (2, 1.8), (3, 1.2), (4, 0.9), (5, 0.72)

\(Clear[P4, x, xy]; 

    xy = {{1, 3.6}, {2, 1.8}, {3, 1.2}, {4, 0.9}, {5, 0.72}}; 

    Print["Points = ", xy]; 
Print["The Newton Polynomial is:"]; 

    P4[x_] = NewtonP[xy]; 
Print[\*""\!\(P\_4\)[x] = "", P4[x]]; 

    Print[\*""\!\(P\_4\)[2.5] = "", P4[2.5]]; \)

The divided difference table for this example is:

TableForm[d]

Chapter 4, Example 4.13, Page 231.
Construct Netwon polynomials of degree  n = 1,2,3  for f(x) = cos(x).

\(Clear[f, P1, i, x, xy]; 
f[x_] = Cos[x]; 

    xy = Table[{i, N[f[i]]}, {i, 0, 1}]; 
Print["f[x] = ", f[x]]; 

    Print["Points = ", xy]; 
Print["The Newton Polynomial is:"]; 

    P1[x_] = NewtonP[xy]; 
Print[\*""\!\(P\_1\)[x] = "", P1[x]]; \)

\(dots = 
      ListPlot[xy, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    g1 = Plot[{f[x], P1[x]}, {x, 0, 3.4}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[g1, dots, PlotRange -> {{0, 3.4}, {\(-1.1\), 1.1}}, 
      PlotLabel -> \*""f[x] and \!\(P\_1\)[x]"", 
      AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Newton polynomial is:"]; 

    Print[\*""\!\(P\_1\)[x] = "", P1[x]]; \)

\(Clear[f, P2, i, x, xy]; 
f[x_] = Cos[x]; 

    xy = Table[{i, N[f[i]]}, {i, 0, 2}]; 
Print["f[x] = ", f[x]]; 

    Print["Points = ", xy]; 
Print["The Newton Polynomial is:"]; 

    P2[x_] = NewtonP[xy]; 
Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

\(dots = 
      ListPlot[xy, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    g2 = Plot[{f[x], P2[x]}, {x, 0, 3.4}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[g2, dots, PlotRange -> {{0, 3.4}, {\(-1.1\), 1.1}}, 
      PlotLabel -> \*""f[x] and \!\(P\_2\)[x]"", 
      AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Newton polynomial is:"]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

\(Clear[f, P3, i, x, xy]; 
f[x_] = Cos[x]; 

    xy = Table[{i, N[f[i]]}, {i, 0, 3}]; 
Print["f[x] = ", f[x]]; 

    Print["Points = ", xy]; 
Print["The Newton Polynomial is:"]; 

    P3[x_] = NewtonP[xy]; 
Print[\*""\!\(P\_3\)[x] = "", P3[x]]; \)

dots = 
      ListPlot[xy, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; pointset = Point/@xy; graphset = Prepend[pointset, PointSize[0.02]]; \(dots = Graphics[graphset]; 

    g3 = Plot[{f[x], P3[x]}, {x, 0, 3.4}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[g3, dots, PlotRange -> {{0, 3.4}, {\(-1.1\), 1.1}}, 
      PlotLabel -> \*""f[x] and \!\(P\_3\)[x]"", 
      AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Newton polynomial is:"]; 

    Print[\*""\!\(P\_3\)[x] = "", P3[x]]; \)

Algorithm 4.6  Chebyshev Approximation. To construct and evaluate the Chebyshev interpolating polynomial of degree  n  over the interval  [-1, 1],  where  P (x) = (∑_(j = 0)^n)(Subscript[c, j]) (Subscript[T, j]) (x)  is based on the nodes   (Subscript[x, k])  =  cos ((((2  k + 1) π)/(2  n + 2))).
Section  4.5   Chebyshev Polynomials   Page  246

First, execute the following group of cells:

Execute this cell to activate Chebyshev Approximation.
Cpolys[n_]
Chebyshev[f0_,n_]
ChebyshevPoly[f0_,n_]
f0_ is the function f(x),
n_ is the degree.

Off[General::"spell1"]; Clear[Cpolys, n]; \(\(Cpolys[n_] := 
	
      Block[{j}, 
		T[0, x_] = 1; 
		T[1, x_] = x; 
		
        For[j = 1, j ≤ n - 1, \(j++\), 
			
          Tj1x = 2\ x\ T[j, x] - T[j - 1, x]; 
			
          T[j + 1, x_] = Simplify[Tj1x]]; 
		
        For[j = 0, j ≤ n, \(j++\), 
			
          \(Print["T"\_j, "(x) = ", T[j, x]]; \)]; ]; \)\)

\(Clear[Chebyshev, f0, n]; 

    Chebyshev[f0_, n_] := 
	
      Block[{f, j, k, sum}, 
		Set@@{f[x_], f0}; 
		Clear[X, Y, c, T]; 
        
		X = Table[0, {n + 1}]; 
		Y = Table[0, {n + 1}]; 
		
        h = π\/\(2\ n + 2\); 
		
        For[k = 0, k ≤ n, \(k++\), 
			
          X[[k + 1]] = 
            Cos[\((2\ k + 1)\)\ h]; 
			
          Y[[k + 1]] = 
            f[X[[k + 1]]]; 
			
          c[k] = 0]; 
		
        For[k = 0, k ≤ n, \(k++\), 
			Z = \((2\ k + 1)\)\ h; 
          
			For[j = 0, j ≤ n, \(j++\), 
				
            c[j] = c[j] + 
                Y[[k + 1]]\ Cos[j\ Z]]]; 
        
		c[0] = c[0]\/\(n + 1\); 
		
        For[j = 1, j ≤ n, \(j++\), 
			
          c[j] = \(2\ c[j]\)\/\(n + 1\)]; 
		sum = 0; 
		T[0, x_] = 1; 
        
		T[1, x_] = x; 
		
        For[j = 1, j ≤ n - 1, \(j++\), 
			
          Tj1x = 2\ x\ T[j, x] - T[j - 1, x]; 
			
          T[j + 1, x_] = Simplify[Tj1x]]; 
		
        For[j = 0, j ≤ n, \(j++\), 
			
          sum = sum + c[j]\ T[j, x]]; 
		sum = Expand[sum]; 
		
        Return[N[sum]]; ]; 
\(On[General::"spell1"]; \)\)

The subroutine given above uses our recursive formulas for the Chebyshev polynomials.  The following version uses Mathematica's built-in functions.  Both give the same results.

Clear[ChebyshevPoly, f0, n]; \(\(ChebyshevPoly[f0_, n_] := 
      Block[{f, h, j, k}, 
\ \ \ Set@@{f[x_], f0}; 
\ \ \ Clear[X, Y, c]; 

        \ \ \ h = N[π\/\(2\ n + 2\)]; 
\ \ \ 
        X = Table[Cos[\((2\ k + 1)\)\ h], {k, 0, n}]; 
\ \ \ Y = f[X]; 

        \ \ \ c = Table[0, {k, 0, n}]; 
\ \ \ 
        For[k = 0, k ≤ n, \(k++\), 
\ \ \ \ \ \ 
          Z = \((2\ k + 1)\)\ h; 
\ \ \ \ \ \ 
          cosZ = Table[Cos[\((j - 1)\)\ Z], {j, 1, n + 1}]; 
\ \ \ 
          c = c + Y[[k + 1]]\ cosZ]; 

        \ \ \ c[[0 + 1]] = 
          1\/2\ c[[0 + 1]]; 
\ \ \ 
        c = \(2\ c\)\/\(n + 1\); 
\ \ \ 
        CP = Expand[
            ∑\+\(j = 0\)\%n
                 c[[j + 1]]\ 
                ChebyshevT[j, x]]; 
\ \ \ Return[CP]]; \)\)

First, generate some Chebyshev polynomials.

Cpolys[5]

Compare with Mathematica's built in function.

\(\(Print[\*""\!\(T\_5\)(x) = "", ChebyshevT[5, x]]; \)\)

Example 4.12,  Page 205.  Investigate a Chebyshev polynomial for  (e^x).

\(Clear[f, P2, T2, x]; 
f[x_] = Exp[x]; 
P2[x_] = Chebyshev[f[x], 2]; 

    Print[" f[x] = ", f[x]]; 
Print["The Chebyshev polynomial is:"]; 
    
Print[\*""\!\(P\_2\)[x] = "", P2[x]]; 

    Print[\*""\!\(P\_2\)[x] = "", ChebyshevPoly[f[x], 2]]; \)

\(dots = 
      ListPlot[Transpose[{X, Y}], PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    g2 = Plot[{f[x], P2[x]}, {x, \(-1\), 1}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[g2, dots, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Chebyshev polynomial is:"]; 

    Print[\*""\!\(P\_2\)[x] = "", P2[x]]; \)

\(e[x_] = Exp[x] - P2[x]; 

    gr = Plot[e[x], {x, \(-1\), 1}, PlotStyle -> {RGBColor[0, 0, 1]}]; 

    Print["Error for the Chebyshev approximation."]; 

    Print[\*""f[x]-\!\(P\_2\)[x] = "", e[x]]; \)

Find the error bound for the Chebyshev approximation.
The error at the right endpoint.

\(xr\  = \ 1.0; 
yr\  = \ e[xr]; 

    Print["f[", xr, \*""] - \!\(P\_2\)["", xr, "] = ", yr]; \)

The error at the local minimum.

\(sol = FindMinimum[e[x], {x, 0.5}]; 

    xm\  = \ sol[[2, 1, 2]]; 

    ym\  = \ sol[[1]]; 

    Print["f[", xm, \*""] - \!\(P\_2\)["", xm, "] = ", ym]; \)

The error bound for the Chebyshev polynomial is:

\(eset\  = \ {Abs[ym], Abs[yr]}; 
bdd\  = \ Max[Abs[ym], Abs[yr]]; 

    Print["Error bound for the Chebyshev approximation."]; 

    Print["Error Bound = max ", eset]; 

    Print["Error Bound = ", bdd]; \)

Compare with the Taylor polynomial approximation, and compare the error bounds.

\(T2[x_] = Normal[Series[Exp[x], {x, 0, 2}]]; 

    Print[" f[x] = ", f[x]]; 
Print["The Taylor polynomial is:"]; 

    Print[\*""\!\(T\_2\)[x] = "", T2[x]]; \)

\(dots = 
      ListPlot[Transpose[{X, Y}], PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    g2 = Plot[{f[x], T2[x]}, {x, \(-1\), 1}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity]; 

    Show[g2, dots, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is f[x] = ", f[x]]; 

    Print["The Taylor polynomial is:"]; 

    Print[\*""\!\(T\_2\)[x] = "", T2[x]]; \)

\(e[x_] = Exp[x] - T2[x]; 

    gr = Plot[e[x], {x, \(-1\), 1}, PlotStyle -> {RGBColor[0, 0, 1]}, 
        PlotRange -> {{\(-1\), 1}, {\(-0.15\), 0.25}}]; 

    Print["Error for the Taylor approximation."]; 

    Print[\*""f[x]-\!\(T\_2\)[x] = "", e[x]]; \)

Find the error bound for the Chebyshev approximation.
The error at the left endpoint.

\(xl\  = \ \(-1.0\); 
yl\  = \ e[xl]; 

    Print["f[", xl, \*""] - \!\(T\_2\)["", xl, "] = ", yl]; \)

The error at the right endpoint.

\(xr\  = \ 1.0; 
yr\  = \ e[xr]; 

    Print["f[", xr, \*""] - \!\(T\_2\)["", xr, "] = ", yr]; \)

The error bound for the Taylor polynomial is:

\(eset\  = \ {Abs[yl], Abs[yr]}; 
bdd\  = \ Max[Abs[yl], Abs[yr]]; 

    Print["Error bound for the Taylor approximation."]; 

    Print["Error Bound = max ", eset]; 

    Print["Error Bound = ", bdd]; \)

The error the Chebyshev polynomial is less than the error for the Taylor polynomial.