Exercise 5.3 Solutions

2.

(a)

\(\(\(\(p[x_]\  = \ a\ x\^3\  + \ b\ x\^2\  + \ c\ x\  + \ d\ ;\)\ 

    \(Print["p[x] = ", p[x]\ ];\)\)\(\ \)\)\)

eq1  =     p[1]  ==  3;  
eq2  =  p'[1]  ==  -4;  
eq3  =     p[2]  ==  1;  
eq4  =  p'[2]  ==  2;  
eqns  =  {eq1, eq2, eq3, eq4};  
vars  =  {a, b, c, d};  
Print[TableForm[eqns] ];

(b)

solset  =  Solve[eqns, vars];

Print[solset];

f[x_]  =  
        ReplaceAll[p[x], 
          (Subscript[solset, [[1]]]) ];  
Needs["Graphics`Colors`"];  
xy = { {1, f[1]}, {2, f[2]} };  
dots  =  
        ListPlot[xy, PlotStyle -> {Red, PointSize[0.02]} , 
          DisplayFunction -> Identity];
gr  = 
        Plot[f[x], {x, 1, 2}, PlotStyle -> Magenta, 
          AspectRatio -> (2.2/3), 
          DisplayFunction -> Identity];
Show[gr, dots, DisplayFunction -> $DisplayFunction];
Print["y = f[x] = ", f[x] ];  
Print[" f[1] = ", f[1], ",      f[2] = ", 
        f[2] ];  
Print["f'[1] = ", f'[1] ,  ",     f'[2] = ", f'[
          2] ];

3a.

f1[x_] := (19/2) - (81/4) x + 
        15  (x^2) - (13/4) (x^3)

f2[x_] := (-77/2) + (207/4) x - 
        21  (x^2) + (11/4) (x^3)

f1[2]

f2[2]

f1'[2]

f2'[2]

f1''[2]

f2''[2]

3c.

f1[x_] := 18 - (75/2) x + 26  (x^2) - (11/2) (x^3)

f2[x_] := -70 + (189/2) x - 
        40  (x^2) + (11/2) (x^3)

f1[2]

f2[2]

f1'[2]

f2'[2]

f1''[2]

f2''[2]

5.

x0 = -3; x1 = -2; x2 = 1; x3 = 4; y0 = 2; y1 = 0; y2 = 3; 
    y3 = 1;

h0 = x1 - x0; h1 = x2 - x1; h2 = x3 - x2;

d0 = (y1 - y0)/h0; d1 = (y2 - y1)/h1; 
    d2 = (y3 - y2)/h2;

u1 = 6 (d1 - d0); u2 = 6 (d2 - d1);

m0 = 0; m3 = 0;

h0 m0  + 2 (h0 + h1) m1 + h1 m2  ==  u1

h1 m1  + 2 (h1  +  h2) m2  + h2 m3 == u2

solns = 
      Solve[{h0 m0  + 2 (h0 + h1) m1 + h1 m2  ==  u1, 
          h1 m1  + 2 (h1  +  h2) m2  + h2 m3 == u2}, {m1, 
          m2}]

s0[x_] := 
      y0 + (d0 - ((h0 (2  m0 + m1))/6)) (x - x0) + (m0/
                2) (((x - x0))^2) + (((m1 - m0)/(6 h0))) (((x - x0))^3) /. solns[[1]]

s1[x_] := 
      y1 + (d1 - ((h1 (2  m1 + m2))/6)) (x - x1) + (m1/
                2) (((x - x1))^2) + (((m2 - m1)/(6 h1))) (((x - x1))^3) /. solns[[1]]

s2[x_] := 
      y2 + (d2 - ((h2 (2  m2 + m3))/6)) (x - x2) + (m2/
                2) (((x - x2))^2) + (((m3 - m2)/(6 h2))) (((x - x2))^3) /. solns[[1]]

s0[x]

s1[x]

s2[x]

Expand[(41/87) (x + 3)^3 - ((215/87) (x + 3)) + 
        2]

s0[x0] == y0

s0[x1] == s1[x1]

8.

x0 = -3; x1 = -2; x2 = 1; x3 = 4; y0 = 2; y1 = 0; y2 = 3; 
    y3 = 1;

h0 = x1 - x0; h1 = x2 - x1; h2 = x3 - x2;

d0 = (y1 - y0)/h0; d1 = (y2 - y1)/h1; 
    d2 = (y3 - y2)/h2;

u1 = 6 (d1 - d0); u2 = 6 (d2 - d1);

m0 = -1; m3 = 2;

h0 m0  + 2 (h0 + h1) m1 + h1 m2  ==  u1

h1 m1  + 2 (h1  +  h2) m2  + h2 m3 == u2

solns = 
      Solve[{h0 m0  + 2 (h0 + h1) m1 + h1 m2  ==  u1, 
          h1 m1  + 2 (h1  +  h2) m2  + h2 m3 == u2}, {m1, 
          m2}]

s0[x_] := 
      y0 + (d0 - ((h0 (2  m0 + m1))/6)) (x - x0) + (m0/
                2) (((x - x0))^2) + (((m1 - m0)/(6 h0))) (((x - x0))^3) /. solns[[1]]

s1[x_] := 
      y1 + (d1 - ((h1 (2  m1 + m2))/6)) (x - x1) + (m1/
                2) (((x - x1))^2) + (((m2 - m1)/(6 h1))) (((x - x1))^3) /. solns[[1]]

s2[x_] := 
      y2 + (d2 - ((h2 (2  m2 + m3))/6)) (x - x2) + (m2/
                2) (((x - x2))^2) + (((m3 - m2)/(6 h2))) (((x - x2))^3) /. solns[[1]]

s0[x]

s1[x]

s2[x]

s0[x0] == y0

s0[x1] == s1[x1]

s1[x2] == s2[x2]

s0'[x1] == s1'[x1]

s1'[x2] == s2'[x2]

s0''[x1] == s1''[x1]

s1''[x2] == s2''[x2]

Functions Used in this notebook

NaturalSpline[XY0_] := Module[{XY = XY0}, 
Differences :=  
                  Module[{k},  
n  =  
                      Length[XY] - 
                        1;  
X  =  (Subscript[(Transpose[XY]), [[1]]]);  
Y  =  (Subscript[(Transpose[XY]), [[2]]]);  
h = d  =  
                        Table[0, {n}];  
m  =  
                      Table[0, {n + 
                            1}];  
a = b = c = v  =  Table[0, {n - 1}];  
s  =  
                      Table[0, {n}, {4}];  
(Subscript[h, [[1]]])  =  
                      (Subscript[X, [[2]]])  -  
                        (Subscript[X, [[1]]]);  
(Subscript[d, [[1]]])  =  (((Subscript[Y, [[2]]])  -  (Subscript[Y, [[1]]]))/(Subscript[h, [[1]]]));  
For[ k = 2,  
                      k ≤ 
                        n,  k++,  
(Subscript[h, [[k]]])  =  
                        (Subscript[X, [[k + 1]]])  -  
                          (Subscript[X, [[k]]]);  
(Subscript[d, [[k]]])  =  (((Subscript[Y, [[k + 1]]])  -  (Subscript[Y, [[k]]]))/(Subscript[h, [[k]]]));  
(Subscript[a, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[b, [[k - 1]]])  =  
                        2 ((Subscript[h, [[k - 1]]])  +  
                              (Subscript[h, [[k]]]));  
(Subscript[c, [[k - 1]]])  =  
                        (Subscript[h, [[k]]]);  
(Subscript[v, [[k - 1]]])  =  
                        6 ((Subscript[d, [[k]]])  -  
                              (Subscript[d, [[k - 1]]]))]; ]; 
Print[h];
Print[d];
TriDiagonal := 
                  Module[{k, 
                      t},  
(Subscript[m, [[1]]])  =  
                      0;  
(Subscript[m, [[n + 1]]])  =  0; 
For[ k = 2,  
                      k ≤ 
                        n - 1,  k++,  
t  =  
                        ((Subscript[a, [[k - 1]]])/(Subscript[b, [[k - 1]]]));  
(Subscript[b, [[k]]])   =   (Subscript[b, [[k]]])  -  
                          t  (Subscript[c, [[k - 1]]]);  
(Subscript[v, [[k]]])   =   (Subscript[v, [[k]]])  -  
                          t  (Subscript[v, [[k - 1]]]); ];  
(Subscript[m, [[n]]])   =   ((Subscript[v, [[n - 1]]])/(Subscript[b, [[n - 1]]]));  
For[ k = n - 2,  
                      1 ≤ 
                        k,  k--,  
(Subscript[m, [[k + 1]]])  =  (((Subscript[v, [[k]]])  -  (Subscript[c, [[k]]]) (Subscript[m, [[k + 2]]]))/(Subscript[b, [[k]]])); ]; ]; 
ComputeCoeff := 
                  Module[{k},  
For[ k = 1,  
                        k ≤ 
                          n,  k++,  
(Subscript[s, [[k, 1]]])  =  
                          (Subscript[Y, [[k]]]);  
(Subscript[s, [[k, 2]]])  =  
                          (Subscript[d, [[k]]])  -  
                            (1/6) (Subscript[h, [[k]]]) (2 (Subscript[m, [[k]]]) + 
                                  (Subscript[m, [[k + 1]]]));  
(Subscript[s, [[k, 3]]])  =  
                          ((Subscript[m, [[k]]])/2);  
(Subscript[s, [[k, 4]]])  =  (((Subscript[m, [[k + 1]]])  -  (Subscript[m, [[k]]]))/(6 (Subscript[h, [[k]]]))); ]; ]; 
CS[t_] := 
                  Module[{j},  
For[ j = 1,  
                      j ≤ 
                        n,  j++,  
If[ 
                          (Subscript[X, [[j]]]) ≤ t  &&  
                            t < 
                              (Subscript[X, [[j + 1]]]), k = j ]; ];  
If[
                      t < (Subscript[X, [[1]]]), k = 1 ];  
If[
                      (Subscript[X, [[n + 1]]]) ≤ t, 
                      k = n ];  
w  =  
                      t  -  
                        (Subscript[X, [[k]]]);  
Return[ (((Subscript[s, [[k, 4]]]) w + 
                                    (Subscript[s, [[k, 3]]])) w + 
                              (Subscript[s, [[k, 2]]])) w + 
                        (Subscript[s, [[k, 1]]]) ]; ]; 
 (*The main part of the program starts here*) ;  
Differences; 
TriDiagonal; 
ComputeCoeff; ];