Homework, Chapter 5

Andy Long
MAT360
Spring, 2017

p. 190, #1, 2, 8
    p. 206, #2, 4
    p. 221, #7, 10

p. 190, #1, 2 and 8

Lagrange Polynomial Method -- p. 190, 1a

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=0
x2=1
x3=2
y1=0
y2=-1
y3=4

qL[x]

Plot[{qL[x]},{x,0,Pi}]

p1=Plot[{qL[x]},{x,Min[{x1,x2,x3}],Max[{x1,x2,x3}]}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3}},PlotMarkers->{Automatic,Large}];
Show[p1,p2]

Lagrange Polynomial Method -- p. 190, 1b

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

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

x1=1
x2=3
x3=6
x4=10
y1=1
y2=2
y3=3
y4=4

Plot[{qL[x]},{x,0,Pi}]

p1=Plot[{qL[x]},{x,Min[{x1,x2,x3,x4}],Max[{x1,x2,x3,x4}]}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3},{x4,y4}},PlotMarkers->{Automatic,Large}];
Show[p1,p2]

p. 190, #8:

If the points can be modelled by a lower degree polynomial (e.g. three co-linear points), then the Lagrange polynomial made up of polynomials of degree n+1 will not end up as n+1th degree. 

If the sum of the coefficients is 0, then the leading term will have coefficient 0, and so the polynomial will, in fact, be of lower degree.

p. 206, #2 and 4

#2:

The absolute error is bounded by |f’’(xi)/2!*(x-x0)(x-x1)|

f[x_]:=Log[x]
f''[x]

So we bound the largest absolute error by |(1/1.1)^2/2!*(1.13-1.1)(1.13-1.2)|

Abs[(1/1.1)^2/2!*(1.13-1.1)(1.13-1.2)]

Since we know that the actual absolute error is actually negative (by concavity) we bound the minimum absolute error by |(1/1.2)^2/2!*(1.13-1.1)(1.13-1.2)|

Abs[(1/1.2)^2/2!*(1.13-1.1)(1.13-1.2)]

the actual absolute error is

Plot[{(f[1.2]-f[1.1])/(1.2-1.1)*(x-1.1)+f[1.1],f[x]},{x,1.1,1.2}]
Abs[(f[1.2]-f[1.1])/(1.2-1.1)*(1.13-1.1)+f[1.1]-f[1.13]]

#4:

The relative error is given as |(e^x-p(x))/e^x|. According to our Theorem 5.3, p. 193, that’s 

Abs[(x-x0)(x-x1)(x-x2)(x-x3)*Exp[xi[x]]/(4!*Exp[x])]

Now we bound everything in sight: the polynomial in x has to have its maximum at the center of the interval between the two middle points. If the stepsize is h, then the polynomial piece contributes
(3h/2)(h/2)(h/2)(3h/2)

The exponentials are maximized if x is at one end of the middle interval, and xi(x) is at the opposite end, giving a difference of 2h in values, or .2. Cancelling we get e^(.2). 

All together I get a bound, for h=.1 of

h=.1;
Exp[.2]/Factorial[4]*9*h^4

Dropping out the Exp[0.2], we get

1/Factorial[4]*9*h^4/16

which is actually less than

1.0/60000

p. 221, #7 and 10

p. 221, #7: make sure that you plot and check your work. How well is your “interpolator” doing? The cool thing about Newton’s interpolator is that you can simply add an additional point to an existing interpolator very easily, using these divided difference tables. 

That’s a nice feature. You can easily build off what you’ve already got.

dividedDifference[list_]:=
Module[
{n=Length[list],b},
If[n==1,
b=list[[1]][[2]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
Print[list," ",b,"=",N[b]];
b
]
a=dividedDifference[{{3,2}}]
b=dividedDifference[{{3,2},{6,3}}]
c=dividedDifference[{{3,2},{6,3},{1,1}}]
d=dividedDifference[{{1,1},{3,2},{6,3},{10,4}}]
p[x_]:=a+(x-3)*(b+(x-6)*(c+(x-1)*d))
N[p[{1,3,6,10}]]
p1=Plot[{p[x]},{x,1,10}];
p2=ListPlot[{{1,1},{3,2},{6,3},{10,4}},PlotMarkers->{Automatic,Large}];
Show[p1,p2]

p. 221, #10

x1=2.0
x2=2.2
x3=2.4
x4=2.6
x5=2.8
x6=3.0
y1=.30103
y2=.34242
y3=.38021
y4=.41497
y5=.44716
y6=.47712

a=dividedDifference[{{x1,y1}}]
b=dividedDifference[{{x1,y1},{x2,y2}}]
c=dividedDifference[{{x1,y1},{x2,y2},{x3,y3}}]
d=dividedDifference[{{x1,y1},{x2,y2},{x3,y3},{x4,y4}}]
e=dividedDifference[{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5}}]
p[x_]:=a+(x-x1)*(b+(x-x2)*(c+(x-x3)*(d+(x-x4)*e)))

{a,b,c,d,e}
p[2.1]
Log[10,2.1]
(p[2.1]-Log[10,2.1])/Log[10,2.1]

p1=Plot[{p[x]},{x,x1,x6}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3},{x4,y4}},PlotMarkers->{Automatic,Large}];
Show[p1,p2]

dividedDifference[list_]:=
Module[
{n=Length[list],b},
If[n==1,
b=list[[1]][[2]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
Print[list," ",b];
b
]
dividedDifference[{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6}}];