(*^ ::[paletteColors = 128; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L3, e8, 24, "New York"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "New York"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 14, "New York"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 14, "New York"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 12, "New York"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 10, "New York"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "New York"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; ; fontset = name, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = Left Header, inactive, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; ; fontset = Left Footer, inactive, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Geneva"; ; fontset = clipboard, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = completions, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special1, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special2, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, center, M7, 12, "New York"; ; fontset = special3, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, right, M7, 12, "New York"; ; fontset = special4, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special5, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ;] :[font = text; inactive; locked; preserveAspect; center; ] 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/ ;[s] 6:0,1;67,0;92,1;157,3;159,1;232,2;486,-1; 4:1,17,12,New York,0,12,0,0,0;3,23,17,New York,0,18,0,0,0;1,19,14,New York,0,14,0,0,0;1,27,18,New York,32,14,0,0,0; :[font = text; inactive; locked; preserveAspect; center; ] CONTENTS ;[s] 2:0,1;8,0;9,-1; 2:1,17,12,New York,0,12,0,0,0;1,23,17,New York,1,18,0,0,0; :[font = text; inactive; locked; preserveAspect; ] 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 ;[s] 2:0,1;53,0;323,-1; 2:1,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 25:0,1;47,0;76,2;77,5;80,6;82,7;83,0;85,4;88,0;90,3;91,0;96,3;97,0;98,4;99,0;107,2;108,5;111,6;113,7;114,0;123,3;124,0;125,4;126,0;187,1;250,-1; 8:9,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;2,26,18,Symbol,1,18,0,0,0;3,25,16,New York,64,12,0,0,0;3,25,16,New York,32,12,0,0,0;2,20,13,New York,64,9,0,0,0;2,20,13,New York,32,9,0,0,0;2,27,17,Symbol,32,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] 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 . :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] 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]; ]; :[font = input; dontPreserveAspect; endGroup; ] 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]; :[font = text; inactive; preserveAspect; ] Chapter 4, Example for Taylor polynomials, Page 203-204. ;[s] 3:0,0;1,1;57,0;59,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; preserveAspect; ] TaylorCoeff[Sin[x],0,11] :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.01] :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.0001] :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.000001] :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.00000001] :[font = text; inactive; preserveAspect; ] 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): ;[s] 3:0,0;255,1;266,0;363,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] P[x_]= Normal[Series[Sin[x],{x,0,11}]] :[font = input; preserveAspect; ] P[1] :[font = text; inactive; preserveAspect; ] Notice that all the terms in the Taylor polynomial were added in the computation of P[1]. Compare this with the previous method: :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.000001] :[font = input; preserveAspect; ] SumTaylor[c,1,0,0.00000001] :[font = input; preserveAspect; ] 1/39916800 //N :[font = text; inactive; preserveAspect; ] 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. :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 5:0,1;36,0;107,2;108,0;152,1;206,-1; 3:2,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;1,17,12,Symbol,1,14,0,0,0; :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] 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]; ]; :[font = input; dontPreserveAspect; ] 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]; ]; :[font = input; dontPreserveAspect; endGroup; ] 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]; :[font = text; inactive; preserveAspect; ] This example is for illustration purposes only. In practice one should modify the procedure to input the coefficient c. :[font = input; dontPreserveAspect; ] Clear[P,x]; P[x_] = 3x^3 - 2x^2 + 5x - 7 :[font = text; inactive; preserveAspect; ] Put the polynomial in nested form. :[font = input; dontPreserveAspect; ] Clear[Q,x]; Q[x_] = Horner[P[x],x] :[font = text; inactive; preserveAspect; ] We can verify that Mathematica has stored the two forms of the polynomial differently by looking at them with FullForm. ;[s] 3:0,0;20,1;31,0;122,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] FullForm[P[x]] :[font = input; preserveAspect; ] FullForm[Q[x]] :[font = text; inactive; preserveAspect; ] Evaluate the nested form Q(1) and compare with P(1). :[font = input; dontPreserveAspect; ] Q[1] :[font = input; preserveAspect; ] P[1] :[font = text; inactive; preserveAspect; ] Verify that the two forms of the polynomial are equivalent by forming the difference and simplifying. :[font = input; preserveAspect; ] Clear[R,x] R[x_] = Q[x] - P[x] :[font = input; dontPreserveAspect; ] Simplify[R[x]] :[font = text; inactive; preserveAspect; ] Now consider the derivative of the polynomial P(x). :[font = input; dontPreserveAspect; ] P[x] :[font = input; dontPreserveAspect; ] P'[x] :[font = text; inactive; preserveAspect; ] Put the derivative of the polynomial in nested form. :[font = input; dontPreserveAspect; ] dP[x_] = DiffPoly[P[x],x] :[font = text; inactive; preserveAspect; ] Evaluate the derivative of P(x) at a specified value of x. :[font = input; dontPreserveAspect; ] dP[1] :[font = input; preserveAspect; ] P'[1] :[font = text; inactive; preserveAspect; ] Verify that the two forms of the derivative are equivalent. :[font = input; dontPreserveAspect; ] P'[x] - dP[x] :[font = input; dontPreserveAspect; ] P'[x] - dP[x]//Simplify :[font = text; inactive; preserveAspect; ] Now consider the integral of the polynomial P(x). :[font = input; dontPreserveAspect; ] P[x] :[font = input; dontPreserveAspect; ] Clear[V,x]; V[x_] = Integrate[P[x],x] :[font = text; inactive; preserveAspect; ] Put the integral of the polynomial in nested form. :[font = input; dontPreserveAspect; ] iP[x_] = IntPoly[P[x],x] :[font = text; inactive; preserveAspect; ] Evaluate the integral of P(x) at x = 1. :[font = input; dontPreserveAspect; ] iP[1] :[font = input; preserveAspect; ] V[1] :[font = text; inactive; preserveAspect; ] Verify that the two forms of the integral are equivalent. :[font = input; dontPreserveAspect; ] V[x] - iP[x] :[font = input; preserveAspect; ] V[x] - iP[x]//Simplify :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 15:0,1;40,0;103,2;104,3;107,4;110,0;112,3;113,0;114,3;117,0;169,3;170,0;172,3;173,0;198,1;245,-1; 5:6,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;1,26,18,Symbol,1,18,0,0,0;5,20,13,New York,64,9,0,0,0;1,20,13,New York,32,9,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] (x - x0)...(x - xk-1)(x - xk+1)...(x - xn) Ln,k(x) = ------------------------------------------- . (xk -x0)...(xk -xk-1)(xk -xk+1)...(xk -xn) ;[s] 31:0,0;1,1;18,3;19,1;29,3;32,1;39,3;42,1;52,3;53,1;56,3;59,1;124,3;125,1;128,3;129,1;135,3;136,2;137,1;139,3;142,1;145,3;146,2;147,1;149,3;152,1;158,3;159,2;160,1;162,3;163,1;176,-1; 4:1,17,12,New York,0,12,0,0,0;14,17,12,Monaco,0,12,0,0,0;3,25,16,Monaco,64,12,0,0,0;13,19,12,Monaco,64,9,0,0,0; :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] 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. ;[s] 3:0,0;1,1;34,0;178,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] Clear[f,x]; f[x_] = Cos[x] XY = {{0.0,f[0.0]},{1.2,f[1.2]}}; N[XY,10] :[font = input; dontPreserveAspect; ] Clear[P1,x]; P1[x_] = N[Lagrange[XY,x],10] :[font = input; dontPreserveAspect; ] Q1[x_] = N[Expand[P1[x]],10] :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in Fit procedure. ;[s] 3:0,0;14,1;25,0;53,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] Fit[XY,{1,x},x] :[font = input; preserveAspect; ] Q1[x] :[font = text; inactive; preserveAspect; ] Construct the other Lagrange polynomials of degree n=1. :[font = input; dontPreserveAspect; startGroup; ] Clear[f,x]; f[x_] = Cos[x] XY = {{0.2,f[0.2]},{1.0,f[1.0]}}; N[XY,10] :[font = input; dontPreserveAspect; endGroup; ] Clear[P2,x]; P2[x_] = N[Lagrange[XY,x],10] :[font = input; dontPreserveAspect; ] Q2[x_] = N[Expand[P2[x]],10] :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] Notice that the two polynomials of degree n=1 were different. :[font = input; preserveAspect; ] P1[x] :[font = input; preserveAspect; ] P2[x] :[font = text; inactive; preserveAspect; ] 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). ;[s] 3:0,0;1,1;35,0;130,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] 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] :[font = input; dontPreserveAspect; ] Clear[P2,x]; P2[x_] = N[Lagrange[XY,x],10] :[font = input; dontPreserveAspect; ] Q2[x_] = N[Expand[P2[x]],10] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's Fit procedure. ;[s] 3:0,0;14,1;25,0;44,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] Fit[XY,{1,x,x^2},x] :[font = input; dontPreserveAspect; ] 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]; :[font = input; dontPreserveAspect; ] 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] :[font = input; dontPreserveAspect; ] Clear[P3,x]; P3[x_] = N[Lagrange[XY,x],10] :[font = input; dontPreserveAspect; ] Q3[x_] = N[Expand[P3[x]],10] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's Fit procedure. ;[s] 3:0,0;14,1;25,0;44,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] Fit[XY,{1,x,x^2,x^3},x] :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 37:0,1;60,0;116,2;117,0;123,2;124,0;130,2;131,0;138,2;139,0;145,2;146,0;153,2;154,0;161,2;162,0;168,2;169,0;176,2;177,0;184,2;185,0;207,2;208,0;214,2;215,0;222,2;223,0;233,2;236,0;265,2;266,0;268,2;269,0;276,2;279,0;281,1;324,-1; 3:18,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;17,20,13,New York,64,9,0,0,0; :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] What is nested multiplication? See Page 228. :[font = input; dontPreserveAspect; ] Clear[A,X,a,x]; A = Table[a[i],{i,5}] X = Table[x[i],{i,4}] :[font = input; dontPreserveAspect; ] NestMult[A,X,1,t] :[font = input; dontPreserveAspect; ] NestMult[A,X,2,t] :[font = input; dontPreserveAspect; ] NestMult[A,X,3,t] :[font = text; inactive; preserveAspect; ] 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 ;[s] 5:0,0;1,1;35,0;77,2;78,0;131,-1; 3:3,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0;1,20,13,New York,64,10,0,0,0; :[font = input; dontPreserveAspect; ] Clear[A,X]; A = {5,-2,0.5,-0.1}; X = {1,3,4};; :[font = input; dontPreserveAspect; ] NestMult[A,X,3,2.5] :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 47:0,1;48,0;109,3;110,0;137,2;138,0;141,2;142,0;148,2;149,0;154,2;155,0;189,2;192,0;198,2;201,0;207,2;208,0;215,2;218,0;224,2;225,0;232,2;233,0;257,2;260,0;266,2;267,0;274,2;275,0;285,2;288,0;302,2;305,0;311,2;312,0;324,2;327,0;334,2;339,0;343,2;350,0;354,2;355,0;359,2;362,0;366,1;409,-1; 4:23,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;21,20,13,New York,64,9,0,0,0;1,14,9,Times,0,12,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] 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. ;[s] 11:0,0;93,2;94,0;97,2;98,0;107,2;110,0;221,2;222,0;228,1;239,0;393,-1; 3:6,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0;4,20,13,New York,64,9,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] Remark: This is not the general Newton polynomial where the abscissas are arbitrary, i.e. not necessarily equally spaced. :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] Chapter 4, Exercise 7, Page 235. Find the Newton polynomial of degree n=4 for the following data: ;[s] 3:0,0;1,1;33,0;101,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; preserveAspect; ] Clear[f,x]; f[x_] = 3.6/x xy = Table[{x,f[x]},{x,1,5}] :[font = input; dontPreserveAspect; ] Clear[P4,x]; P4[x_] = Newton[xy] :[font = input; preserveAspect; ] P4[2.5] :[font = text; inactive; preserveAspect; ] The divided difference table for this example is: :[font = input; dontPreserveAspect; ] TableForm[d] :[font = text; inactive; preserveAspect; ] Chapter 4, Example 4.13, Page 231. Construct Netwon polynomials of degree n = 1,2,3 for f(x) = cos(x). ;[s] 3:0,0;1,1;35,0;107,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] Clear[f,x]; f[x_] = Cos[x] xy = Table[{i,N[f[i],8]},{i,0,1}] :[font = input; dontPreserveAspect; ] Clear[P1,x]; P1[x_] = N[Newton[xy],10] :[font = input; dontPreserveAspect; ] 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]; :[font = input; dontPreserveAspect; ] Clear[f,x]; f[x_] = Cos[x] xy = Table[{i,N[f[i],8]},{i,0,2}] :[font = input; dontPreserveAspect; ] Clear[P2,x]; P2[x_] = N[Newton[xy],10] :[font = input; dontPreserveAspect; ] 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]; :[font = input; dontPreserveAspect; ] Clear[f,x]; f[x_] = Cos[x] xy = Table[{i,N[f[i],8]},{i,0,3}] :[font = input; dontPreserveAspect; ] Clear[P3,x]; P3[x_] = N[Newton[xy],10] :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; locked; preserveAspect; cellOutline; ] 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 ;[s] 18:0,1;41,0;160,2;161,3;164,4;167,0;169,3;171,0;172,3;173,0;201,3;202,0;210,5;216,6;217,0;218,5;224,0;228,1;274,-1; 7:7,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;1,26,18,Symbol,1,18,0,0,0;4,20,13,New York,64,9,0,0,0;1,20,13,New York,32,9,0,0,0;2,14,10,New York,0,9,0,0,0;1,19,13,Symbol,0,12,0,0,0; :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] 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[" "] ]; ]; :[font = input; dontPreserveAspect; ] 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]; :[font = text; inactive; preserveAspect; ] 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. ;[s] 3:0,0;115,1;126,0;179,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; endGroup; ] 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] ]; :[font = text; inactive; preserveAspect; ] First, generate some Chebyshev polynomials. :[font = input; dontPreserveAspect; ] Cpolys[5] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in function. ;[s] 3:0,0;14,1;25,0;48,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; preserveAspect; ] ChebyshevT[5,x] :[font = text; inactive; dontPreserveAspect; ] Example 4.12, Page 205. Investigate a Chebyshev polynomial for exp(x). ;[s] 3:0,0;1,1;25,0;75,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; preserveAspect; ] Clear[f,x]; f[x_] = Exp[x] :[font = input; dontPreserveAspect; ] Clear[P2,x]; P2[x_] = Chebyshev[f[x],2] :[font = input; dontPreserveAspect; ] ChebyshevPoly[f[x],2] :[font = input; dontPreserveAspect; ] 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]; :[font = input; dontPreserveAspect; ] e[x_] = Exp[x] - P2[x] :[font = input; dontPreserveAspect; ] gr = Plot[e[x],{x,-1,1}]; :[font = input; dontPreserveAspect; ] sol = N[FindMinimum[e[x],{x,0.5}],12] :[font = text; inactive; preserveAspect; ] The error at the right endpoint. :[font = input; dontPreserveAspect; ] N[e[1.],12] :[font = text; inactive; preserveAspect; ] The error at the local minimum. :[font = input; preserveAspect; ] sol[[1]] :[font = text; inactive; preserveAspect; ] The error bound for the Chebyshev polynomial is: :[font = input; preserveAspect; ] Max[Abs[N[e[1.],12]],Abs[sol[[1]]]] :[font = text; inactive; preserveAspect; ] Compare with the error for the Taylor polynomial. :[font = input; dontPreserveAspect; ] T2[x_] = Normal[Series[Exp[x],{x,0,2}]] :[font = input; dontPreserveAspect; ] Plot[{Exp[x],T2[x]},{x,-1,1}]; :[font = input; dontPreserveAspect; ] e[x_] = Exp[x] - T2[x] :[font = input; dontPreserveAspect; ] Plot[e[x],{x,-1,1}, PlotRange->{{-1,1},{-0.15,0.25}}]; :[font = input; dontPreserveAspect; ] N[e[-1.],12] :[font = input; dontPreserveAspect; ] N[e[1.],12] :[font = text; inactive; preserveAspect; ] The error bound for the Taylor polynomial is: :[font = input; preserveAspect; ] Max[Abs[N[e[-1.],12]],Abs[N[e[1.],12]]] :[font = text; inactive; preserveAspect; ] The error the Chebyshev polynomial is less than the error for the Taylor polynomial. ^*)