Four Interpolating Methods

Given three points,  with (distinc) x values x1, x2, x3. Given function f, which determines y values y1, y2, y3. 
Find the interpolating quadratic.

Victoria’s Classic Method

Thanks for the start, Victoria!

Clear[x1,x2,x3,y1,y2,y3,a,b,c]

Solve[
{y1==a*(x1^2)+b*x1+c,
y2==a*(x2^2)+b*x2+c,
y3==a*(x3^2)+b*x3+c},
{a,b,c}]

a=-((-x2 y1+x3 y1+x1 y2-x3 y2-x1 y3+x2 y3)/((x2-x3) ((x1^2)-x1 x2-x1 x3+x2 x3)))
b=-(((x2^2) y1-(x3^2) y1-(x1^2) y2+(x3^2) y2+(x1^2) y3-(x2^2) y3)/((x1-x2) (x1-x3) (x2-x3)))
c=-((-(x2^2) x3 y1+x2 (x3^2) y1+(x1^2) x3 y2-x1 (x3^2) y2-(x1^2) x2 y3+x1 (x2^2) y3)/((x1-x2) (x1-x3) (x2-x3)))

Simplify[Solve[a*(x^2)+b*x+c==0,x]]

x1=1
x2=2
x3=3
f[x_]=Cos[x]
y1=N[f[x1]]
y2=N[f[x2]]
y3=N[f[x3]]

qC[x_]:=a*(x^2)+b*x+c

Plot[{f[x],qC[x]},{x,0,Pi}]
Solve[qC[x]==0,x]

Lagrange Polynomial Method -- p. 185

Clear[x1,x2,x3,y1,y2,y3,a,b,c]

qL[x_]:=y1*(x-x2)(x-x3)/((x1-x2)(x1-x3))+
y2*(x-x1)(x-x3)/((x2-x1)(x2-x3))+
y3*(x-x1)(x-x2)/((x3-x1)(x3-x2))

x1=1
x2=2
x3=3
f[x_]=Cos[x]
y1=N[f[x1]]
y2=N[f[x2]]
y3=N[f[x3]]

Plot[{f[x],qL[x]},{x,0,Pi}]
Solve[qL[x]==0,x]

Newton Polynomial Method - p. 209

Clear[x1,x2,x3,y1,y2,y3,a,b,c]

qN[x_]:=y1+(x-x1)((y2-y1)/(x2-x1)+(x-x2)*(y3-(x3-x1)(y2-y1)/(x2-x1)-y1)/((x3-x2)(x3-x1)))

x1=1
x2=2
x3=3
f[x_]=Cos[x]
y1=N[f[x1]]
y2=N[f[x2]]
y3=N[f[x3]]

Plot[{f[x],qN[x]},{x,0,Pi}]
Solve[qN[x]==0,x]

Horner’s Method (“Rule” - p. 212)

And this is the one that I chose for Muller’s method. I create a polynomial not in x, but rather in (x-x3). 
I choose x3, because I think of that as my best guess.

Clear[x1,x2,x3,y1,y2,y3,a,b,c]

This method uses divided differences, as I will derive in class. This is the recursive definition of p. 215.

dividedDifference[list_]:=
Module[
{n=Length[list]},
If[n==1,
b=list[[1]][[2]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
b
]

x1=1
x2=2
x3=3
f[x_]=Cos[x]
y1=N[f[x1]]
y2=N[f[x2]]
y3=N[f[x3]]

a=dividedDifference[{{x2,y2},{x3,y3},{x1,y1}}]
b=dividedDifference[{{x2,y2},{x3,y3}}]-(x2-x3)*a
c=y3
qH[x_]:=y3+(x-x3)*(b+(x-x3)*a)

Plot[{f[x],qH[x]},{x,0,Pi}]
Solve[qH[x]==0,x]

Another Quadratic Method (generalizing Newton’s method)

Muller’s generalizes Secant; we could generalize Newton’s method straightaway. Assume we know f, f’, f’’, and we have an initial guess. Then we can fit a quadratic to the point, its slope, and its curvature, then use the roots of that. The quadratic is an interpolant -- but of point, slope, and curvature data.

x1=1
f[x_]:=Cos[x] 
fp[x_]:=-Sin[x]
fpp[x_]:=-Cos[x] 
y=N[f[x1]]
yp=N[fp[x1]]
ypp=N[fpp[x1]]

qNC[x_]:=y+(x-x1)yp+(x-x1)^2ypp/2

Plot[{f[x],qNC[x]},{x,0,Pi}]
Solve[qNC[x]==0,x]

Newton Polynomial becomes easier using divided differences:

Clear[x1,x2,x3,y1,y2,y3,a,b,c]

qN[x_]:=y1+(x-x1)(dividedDifference[{{x2,y2},{x1,y1}}]+(x-x2)dividedDifference[{{x3,y3},{x2,y2},{x1,y1}}])

x1=1
x2=2
x3=3
f[x_]=Cos[x]
y1=N[f[x1]]
y2=N[f[x2]]
y3=N[f[x3]]

Plot[{f[x],qN[x]},{x,0,Pi}]
Solve[qN[x]==0,x]

qNN[x_]:=qN[x]+(x-x1)(x-x2)(x-x3)dividedDifference[{{x4,y4},{x3,y3},{x2,y2},{x1,y1}}]

x4=4;y4=N[f[x4]];

Plot[{f[x],qNN[x]},{x,0,5}]
Solve[qNN[x]==0,x]

x5=7;y5=N[f[x5]];
qNNN[x_]:=qNN[x]+(x-x1)(x-x2)(x-x3)(x-x4)dividedDifference[{{x5,y5},{x4,y4},{x3,y3},{x2,y2},{x1,y1}}]

Plot[{f[x],qNNN[x]},{x,0,8}]
Solve[qNNN[x]==0,x]

f[x_]:=Log[1+x]
D[f[x],{x,4}]