Section 4.1: Numerical Calculus

Andy Long, Spring, 2024

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]]
]
newtonCoefs=dividedDifference[{{0,0},{2,20.134},{3,30.452},{5,52.11}},True]

xpoints={1.25,1.6,1};
f[x_]:=Sin[Pi x]
ypoints=f[xpoints];
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],True]

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
]

interpolator[1.4,newtonCoefs,xpoints,True]
f[1.4]

Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

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

Show[
Plot[interpolator[x,newtonCoefs,xpoints,False],{x,Min[xpoints],Max[xpoints]}],
ListPlot[Transpose[{xpoints,ypoints}],PlotStyle->{PointSize->Large}],
PlotRange->All
]

Now let’s construct these functions:

NewtonsMethod[f_,a0_,eps_,nits_,verbose_,lots_] :=
Module[
{a=N[a0], (* Use a local variable (so as not to alter passed arguments) *)
fa = f[N[a0]],
i=0,
output (* Let's collect the output, since we've gone to the trouble of computing it. *)
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<=nits,
a = a-(fa/(f'[a]));
fa=f[a];
i=i+1;
output=Append[output,{i,a,fa}];  (* once again, collect output. *)
]; (* Now let's augment the naive function to include printing some nice formatted output: *)
If[verbose,
Module[{},
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"i","(Subscript[x, k])","f[(Subscript[x, k])]"}}],16]]; 
Print["  x  = ",NumberForm[a,16] ]; 
Print["f(x) = ",NumberForm[fa,16] ]; 
]
];
Return[
If[lots,Transpose[output],a]
]; (* Now I'll actually return the output, so that it can be plotted, etc. *)
]

fixedPoint code:

(* These are options that one can give at game time: 
-- epsilon: the value we consider "zero" (0 by default)
  -- tolerance: the width of the interval we consider "close enough" to the root
		-- verbose: print the details of the convergence
-- debug: print out a little more for the coder....
*)
Options[fixedPoint]={maxits->100,epsilon->0,tolerance->10^(-6),verbose->False,debug->False} ;

(* Here's the function itself: *)
fixedPoint[f_,a0_,OptionsPattern[]] :=

Module[{a=a0,c,fa=f[a0],fc,slope,width,output,i=1,its=OptionValue[maxits],eps=OptionValue[epsilon],tol=OptionValue[tolerance],verbose=OptionValue[verbose],debug=OptionValue[debug]},
width=10000;
(* Update and continue: let's start looping: *)
While[i<=its && width>tol,
(* c is the function value at a: *)
c=fa;
fc=f[c];
width=Abs[c-a];
slope=(fc-fa)/width;

(* Otherwise, store output and do it again: *)
If[i==1,
(* Let's collect the output, since we've gone to the trouble of computing it: *)
output={{a,c,fc,slope}},
output=Append[output,{a,c,fc,slope}]
]; 

(* Did we stomp on a root by good fortune? *)
If[width<=tol, 
If[verbose,
Print["Root at c."];printoutFixedPoint[output,a,c,fc,width,slope];Return[Transpose[output]],
Return[fc] (* if we're converging, this is better than c *)
]
];

(* update: *)
a=c;
fa=fc;
i=i+1;
]; 

(* We've popped out, so print some nice formatted output: *)
If[verbose,
printoutFixedPoint[output,a,c,fc,width,slope];
(* Return the output, so that it can be plotted, etc. *)
Return[Transpose[output]],
Return[fc]
];
]

Let’s generate the values our author uses in this section:

g[x_]:=x-2Pi/3Sin[x]
nooVals=Table[{x,NewtonsMethod[g,x,10^-6,1,False,False]},{x,93/70,105/70,2/70}];
MatrixForm[nooVals]

f[x_]:=2Pi/3Sin[x]
feeVals=Table[{x,fixedPoint[f,N[x],maxits->9]},{x,1/7,1,1/7}];
MatrixForm[feeVals]

First we generate quadratic approximations, using a subset of equally spaced points:

xpoints={nooVals[[1]][[1]],nooVals[[4]][[1]],nooVals[[7]][[1]]};
ypoints={nooVals[[1]][[2]],nooVals[[4]][[2]],nooVals[[7]][[2]]};
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],False];
nuP2[x_]=Simplify[newtonCoefs[[1]]+(x-xpoints[[1]])*(newtonCoefs[[2]]+(x-xpoints[[2]])*newtonCoefs[[3]])];
nuP2[x]

xpoints={feeVals[[1]][[1]],feeVals[[4]][[1]],feeVals[[7]][[1]]};
ypoints={feeVals[[1]][[2]],feeVals[[4]][[2]],feeVals[[7]][[2]]};
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],False];
feeP2[x_]=Simplify[
newtonCoefs[[1]]+(x-xpoints[[1]])*
(newtonCoefs[[2]]+(x-xpoints[[2]])*
newtonCoefs[[3]])];
feeP2[x]

Then we generate 6th degree polynomial approximations, using all of the equally spaced points:

{xpoints,ypoints}=Transpose[nooVals];
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],False];
nuP6[x_]=Simplify[
newtonCoefs[[1]]+(x-xpoints[[1]])*
(newtonCoefs[[2]]+(x-xpoints[[2]])*
(newtonCoefs[[3]]+(x-xpoints[[3]])*
(newtonCoefs[[4]]+(x-xpoints[[4]])*
(newtonCoefs[[5]]+(x-xpoints[[5]])*
(newtonCoefs[[6]]+(x-xpoints[[6]])*
newtonCoefs[[7]])))))];
nuP6[x]

{xpoints,ypoints}=Transpose[feeVals];
newtonCoefs=dividedDifference[Transpose[{xpoints,ypoints}],False];
feeP6[x_]=Simplify[
newtonCoefs[[1]]+(x-xpoints[[1]])*
(newtonCoefs[[2]]+(x-xpoints[[2]])*
(newtonCoefs[[3]]+(x-xpoints[[3]])*
(newtonCoefs[[4]]+(x-xpoints[[4]])*
(newtonCoefs[[5]]+(x-xpoints[[5]])*
(newtonCoefs[[6]]+(x-xpoints[[6]])*
newtonCoefs[[7]])))))];
feeP6[x]

Now let’s do some fake calculus...:)

noo[x_]=nuP2[x];
noo'[1.4]
fee[x_]=feeP2[x];
fee'[0.5]
Noo[x_]:=Integrate[noo[x],x]
Noo[x]
Fee[x_]:=Integrate[fee[x],x]
Fee[x]
Integrate[noo[x],{x,1.4,1.5}]
Integrate[fee[x],{x,0,1}]

noo[x_]=nuP6[x];
noo'[1.4]
fee[x_]=feeP6[x];
fee'[0.5]
Noo[x_]:=Integrate[noo[x],x]
Noo[x]
Fee[x_]:=Integrate[fee[x],x]
Fee[x]
Integrate[noo[x],{x,1.4,1.5}]
Integrate[fee[x],{x,0,1}]