Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

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

Algorithm 4.1  (Evaluation of a Taylor Series).  To approximate

		P(x)  =  Sk=0k=¥ f(k)(x0)(x-x0)k/k!  =  Sk=0k=¥ D(k)(x-x0)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 .

First, execute the following group of cells:

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 = s0/.x->X+x0;
  c = CoefficientList[s1,X];
  Print["Coefficient list for the function"];
  Print[f[x],"  expanded about x = ",x0];
  Print[" "]; 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[sum]; ]; 
On[General::spell1];

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

TaylorCoeff[Sin[x],0,11]

SumTaylor[c,1,0,0.01]

SumTaylor[c,1,0,0.0001]

SumTaylor[c,1,0,0.000001]

SumTaylor[c,1,0,0.00000001]

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):

P[x_]= Normal[Series[Sin[x],{x,0,11}]]

P[1]

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]

1/39916800 //N

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

First, execute the following group of cells:

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.

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

Put the polynomial in nested form.

Clear[Q,x];
Q[x_] = Horner[P[x],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).

Q[1]

P[1]

Verify that the two forms of the polynomial are equivalent by forming the difference and simplifying.

Clear[R,x]
R[x_] = Q[x] - P[x]

Simplify[R[x]]

Now consider the derivative of the polynomial P(x).

P[x]

P'[x]

Put the derivative of the polynomial in nested form.

dP[x_] = DiffPoly[P[x],x]

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

dP[1]

P'[1]

Verify that the two forms of the derivative are equivalent.

P'[x] - dP[x]

P'[x] - dP[x]//Simplify

Now consider the integral of the polynomial P(x).

P[x]

Clear[V,x];
V[x_] = Integrate[P[x],x]

Put the integral of the polynomial in nested form.

iP[x_] = IntPoly[P[x],x]

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

iP[1]

V[1]

Verify that the two forms of the integral are equivalent.

V[x] - iP[x]

V[x] - iP[x]//Simplify

Algorithm 4.3  (Lagrange Approximation).  To evaluate the Lagrange polynomial

               P(x)  =  Sk=0k=nÊykLn,k(x),   of degree  n, 

based on the  n+1  points  (xk,yk)  for  k = 0,1, ...,n. 
Section	4.3	Lagrange Approximation   Page  224

(x - x0)...(x - xk-1)(x - xk+1)...(x - xn)
Ln,k(x) = ------------------------------------------- .
           (xk -x0)...(xk -xk-1)(xk -xk+1)...(xk -xn)

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,x];
f[x_] = Cos[x]
XY = {{0.0,f[0.0]},{1.2,f[1.2]}};
N[XY,10]

Clear[P1,x];
P1[x_] = N[Lagrange[XY,x],10]

Q1[x_] = N[Expand[P1[x]],10]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
gr1 = Plot[{f[x],Q1[x]},{x,0.0,1.2},
DisplayFunction->Identity];
Show[gr1,dots,PlotRange->{{0,1.3},{0.0,1.15}},
AxesLabel->{"x","y"},
PlotLabel->"f(x) and the Lagrange polynomial.",
DisplayFunction->$DisplayFunction];

Compare with Mathematica's built in Fit procedure.

Fit[XY,{1,x},x]

Q1[x]

Construct the other Lagrange polynomials of degree n=1.

Clear[f,x];
f[x_] = Cos[x]
XY = {{0.2,f[0.2]},{1.0,f[1.0]}};
N[XY,10]

Clear[P2,x];
P2[x_] = N[Lagrange[XY,x],10]

Q2[x_] = N[Expand[P2[x]],10]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
gr2 = Plot[{f[x],Q2[x]},{x,0.0,1.2},
DisplayFunction->Identity];
Show[gr2,dots,PlotRange->{{0,1.3},{0.0,1.15}},
AxesLabel->{"x","y"},
PlotLabel->"f(x) and the Lagrange polynomial.",
DisplayFunction->$DisplayFunction];

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

P1[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).

Clear[f,x];
f[x_] = Cos[x]
XY = {{0.0,f[0.0]},{0.6,f[0.6]},{1.2,f[1.2]}};
N[XY,10]

Clear[P2,x];
P2[x_] = N[Lagrange[XY,x],10]

Q2[x_] = N[Expand[P2[x]],10]

Compare with Mathematica's Fit procedure.

Fit[XY,{1,x,x^2},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},
DisplayFunction->Identity];
Show[gr3,dots,PlotRange->{{-0.6,2.1},{-0.8,1.15}},
AxesLabel->{"x","y"},
PlotLabel->"f(x) and the Lagrange polynomial.",
DisplayFunction->$DisplayFunction];

Clear[f,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]}};
N[XY,10]

Clear[P3,x];
P3[x_] = N[Lagrange[XY,x],10]

Q3[x_] = N[Expand[P3[x]],10]

Compare with Mathematica's Fit procedure.

Fit[XY,{1,x,x^2,x^3},x]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
gr4 = Plot[{f[x],Q3[x]},{x,-0.6,2.1},
DisplayFunction->Identity];
Show[gr4,dots,PlotRange->{{-0.6,2.1},{-0.8,1.15}},
AxesLabel->{"x","y"},
PlotLabel->"f(x) and the Lagrange polynomial.",
DisplayFunction->$DisplayFunction];

Algorithm 4.4 (Nested Multiplication with Multiple Centers). To evaluate the Newton form of a polynomial

P(x)  =  a0  +  a1(x - x0)  +  a2(x - x0)(x - x1)  +  a3(x - x0)(x - x1)(x - x2) 
			      + ... +  an(x - x0)(x - x1)...(x - xn-1), 

with multiple centers  x0,x1, ...,xn-1.
Section	4.4	Newton Polynomials   Page  233

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,t]

NestMult[A,X,2,t]

NestMult[A,X,3,t]

Chapter 4, Example 4.11, Page 228.

Using nested multiplication, evaluate

P3(x) = ((-0.1(x - 4) + 0.5)(x - 3) - 2)(x - 1)  +  5

Clear[A,X];
A = {5,-2,0.5,-0.1};
X = {1,3,4};;

NestMult[A,X,3,2.5]

Algorithm 4.5 (Newton Interpolation Polynomial). To construct and evaluate the Newton polynomial of  degree  ² n  that passes through 
(xk, yk) = (xk, f(xk))  for  k=0,1,...,n: 

P(x)  =  d0,0  +  d1,1(x - x0)  +  d2,2(x - x0)(x - x1) 
           + ... +  dn,n(x - x0)(x - x1)...(x - xn-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  x1, x2, ... , xn+1. Then it will simplify the diagonal coefficients.  Then it will take the limit to get derivatives of f(x) at x1.

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.

Off[General::spell1];
Clear[Newton];
Newton[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 = sum/.h->H;
  Return[sum]; ];
On[General::spell1];

Chapter 4, Exercise 7, Page 235.

Find the Newton polynomial of degree n=4 for the following data:

Clear[f,x];
f[x_] = 3.6/x
xy = Table[{x,f[x]},{x,1,5}]

Clear[P4,x];
P4[x_] = Newton[xy]

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,x];
f[x_] = Cos[x]
xy = Table[{i,N[f[i],8]},{i,0,1}]

Clear[P1,x];
P1[x_] = N[Newton[xy],10]

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
g1 = Plot[{f[x],P1[x]},{x,0,3.4},
DisplayFunction->Identity];
Show[g1,dots,PlotRange->{{0,3.4},{-1.1,1.1}},
PlotLabel->"f(x) and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

Clear[f,x];
f[x_] = Cos[x]
xy = Table[{i,N[f[i],8]},{i,0,2}]

Clear[P2,x];
P2[x_] = N[Newton[xy],10]

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
g2 = Plot[{f[x],P2[x]},{x,0,3.4},
DisplayFunction->Identity];
Show[g2,dots,PlotRange->{{0,3.4},{-1.1,1.1}},
PlotLabel->"f(x) and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

Clear[f,x];
f[x_] = Cos[x]
xy = Table[{i,N[f[i],8]},{i,0,3}]

Clear[P3,x];
P3[x_] = N[Newton[xy],10]

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
pointset = Map[Point,xy];
graphset = Prepend[pointset,PointSize[0.02]];
dots = Graphics[graphset];
Print["The function is f[x] = ",f[x]];
g3 = Plot[{f[x],P3[x]},{x,0,3.4},
DisplayFunction->Identity];
Show[g3,dots,PlotRange->{{0,3.4},{-1.1,1.1}},
PlotLabel->"f(x) and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

Algorithm 4.6  (Chebyshev Approximation). To construct and evaluate the Chebyshev interpolating polynomial of degree n over the interval [-1,1], where  P(x) =  Sj=0j=nÊcj Tj(x)  is based on the nodes
xk = cos( (2k+1)p/(2n+2) ).
Section	4.5	Chebyshev Polynomials   Page  246

First, execute the following group of cells:

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",Subscript[j],"(x) = ",T[j,x]];
    Print[" "] ]; ];

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 = Pi/(2n + 2);
  For[k=0,k<=n,k++,
    X[[k+1]] = Cos[(2k+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[sum//N]; ];
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 = Pi/(2n+2)//N;
  X = Table[Cos[(2k+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]] = c[[0+1]]/2;
  c = 2 c/(n+1);
  CP = Sum[c[[j+1]] ChebyshevT[j,x],{j,0,n}]//Expand;
  Return[CP] ];

First, generate some Chebyshev polynomials.

Cpolys[5]

Compare with Mathematica's built in function.

ChebyshevT[5,x]

Example 4.12,  Page 205.

Investigate a Chebyshev polynomial for exp(x).

Clear[f,x];
f[x_] = Exp[x]

Clear[P2,x];
P2[x_] = Chebyshev[f[x],2]

ChebyshevPoly[f[x],2]

dots = ListPlot[Transpose[{X,Y}],
       PlotStyle->{PointSize[0.02]},
       DisplayFunction->Identity];
Print["The function is f[x] = ",f[x]];
g2 = Plot[{f[x],P2[x]},{x,-1,1},
DisplayFunction->Identity];
Show[g2,dots,
PlotLabel->"f(x) and the Chebyshev polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

e[x_] = Exp[x] - P2[x]

gr = Plot[e[x],{x,-1,1}];

sol = N[FindMinimum[e[x],{x,0.5}],12]

The error at the right endpoint.

N[e[1.],12]

The error at the local minimum.

sol[[1]]

The error bound for the Chebyshev polynomial is:

Max[Abs[N[e[1.],12]],Abs[sol[[1]]]]

Compare with the error for the Taylor polynomial.

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

Plot[{Exp[x],T2[x]},{x,-1,1}];

e[x_] = Exp[x] - T2[x]

Plot[e[x],{x,-1,1},
PlotRange->{{-1,1},{-0.15,0.25}}];

N[e[-1.],12]

N[e[1.],12]

The error bound for the Taylor polynomial is:

Max[Abs[N[e[-1.],12]],Abs[N[e[1.],12]]]

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