Homework 3.3: Numerical Calculus

Andy Long, Spring, 2024

#1 and 2, p.  129

#1. Newton Polynomials: Divided differences is how we compute the coefficients of the Newton Polynomial:

Define divided differences:

dividedDifference[data_,verbose_]:=
Module[{n=Length[data],xs,ys,i,j,p},
p=ConstantArray[0,{n,n}];
{xs,ys}=Transpose[data];  (* The data is a list of coordinates *)
p[[All,1]]=ys; (* We initialize the first column of p with the y-values *)
For[j=2,j<=n+1,j++,
For[i=1,i<=n+1-j,i++, 
p[[i,j]]= (* Successive columns are filled with divided differences, 
from the previous column. *)
(p[[i+1,j-1]]-p[[i,j-1]])/(xs[[i+j-1]]-xs[[i]]);
];
];
If[verbose,Print[MatrixForm[p]]];
p[[1,All]]
]

#2. Let’s test out our code:

xpoints={0.0,1.0,2.0};
ypoints=E^xpoints;
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],True]

It appears to be working. At least it matches what’s given on p. 124.

This will be useful: Hornerize the Newton Polynomial

interpolator[x_,coefs_,xs_,verbose_]:=
Module[
{result,n=Length[xs]},
result=coefs[[n]];
For[i=n-1,i>0,i--,
If[verbose,Print[{xs[[i]],coefs[[i]]}]];
result=coefs[[i]]+(x-xs[[i]])result;
];
result
]

#3, p.  129

newtonCoefs=dividedDifference[{{0.1,0.12},{0.2,0.14},{0.3,0.13},{0.4,0.15}, {0.5,0.11}},True]

data={{0.1,0.12},{0.2,0.14},{0.3,0.13},{0.4,0.15}, {0.5,0.11}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#4, p.  129 

(solution p. 277 (288))

data={{1,2},{2,2},{3,0},{4,0}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#5, p.  129 

(solution p. 277 (288))

data={{0,10},{30,58},{1029,-32.0}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#6, p.  129

data={{1,0.987},{2.2,-0.123},{3,0.432}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#7, p.  129

data={{1.2,2.2},{1.4,2.1},{1.6,2.3}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#8, p.  129

data={{1,1},{2,3},{3,5},{4,7}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

Well, looking at the data, the expected degree is 1 -- so it depends on what one means by “expected”...:)

#9, p.  129

data={{2,3},{3,5},{4,4}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#10, p.  129

data={{0.0,-6.00000},{0.1,-5.89483},{0.3,-5.65014},{0.6,-5.17788},{1.0,-4.28172},{1.1,-3.9958}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

#11, p.  129

data={{8.3,17.56492},{8.6,18.50510}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
estimate=interpolator[8.4,newtonCoefs,xpoints,False]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
ListPlot[{{8.4,estimate}},PlotStyle->{Orange,PointSize->Large}],
PlotRange->All
]

data={{8.1,16.94410},{8.3,17.56492},{8.6,18.50510}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
estimate=interpolator[8.4,newtonCoefs,xpoints,False]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
ListPlot[{{8.4,estimate}},PlotStyle->{Orange,PointSize->Large}],
PlotRange->All
]

data={{8.1,16.94410},{8.3,17.56492},{8.6,18.50510},{8.7,18.82091}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
estimate=interpolator[8.4,newtonCoefs,xpoints,False]
Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
ListPlot[{{8.4,estimate}},PlotStyle->{Orange,PointSize->Large}],
PlotRange->All
]

#12, p.  129

Here’s our author’s solution (p. 278):

data={{1,0.987},{2.2,-0.123},{3,0.432}};
{xpoints,ypoints}=Transpose[data];
newtonCoefs=dividedDifference[data,True]
ystar=interpolator[2,newtonCoefs,xpoints,False]
Show[
Plot[{Abs[2(x−1)(x−2.2)(x−3)/3!],2/30},{x,1,3},PlotLabel->"Error bound at x=2.0"],
ListPlot[{{2,ystar}}],
PlotRange->All
]