Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 1.	Preliminaries
	Theorem 1.1    (Limits and Continuous Functions).
	Theorem 1.2    (Intermediate Value Theorem).
	Theorem 1.3    (Extreme Value Theorem for a Continuous Function).
	Theorem 1.4    (Differentiable function implies continuous function).
	Theorem 1.5    (Rolle's Theorem).
	Theorem 1.6    (Mean Value Theorem).
	Theorem 1.7    (Extreme Value Theorem for a Differentiable Function).
	Theorem 1.8    (Generalized Rolle's Theorem).
	Theorem 1.9    (First Fundamental Theorem)
	Theorem 1.10  (Second Fundamental Theorem).
	Theorem 1.11  (Mean Value Theorem for Integrals)
	Theorem 1.12  (Weighted Integral Mean Value Theorem).
	Theorem 1.13  (Taylor's Theorem).
	Theorem 1.14  (Horner's Method for Polynomial Evaluation).
	Theorem 1.15  (Geometric Series).
	Theorem 1.16  (Big "O" remainders for Taylor's Theorem).
	Theorem 1.17  (Remainder term for Taylor's Theorem).

Theorem 1.1  (Limits and Continuous Functions).
Assume that  f  is defined on the set S and x is an element of S.
The following statements are equivalent:
(i)   The function  f  is continuous at x0 .
(ii)  If        lim  xn  = x0         as   n->oo, 
       then   lim  f(xn) = f(x0)   as   n->oo. 
Section 1.1,  Review of Calculus, Page 4

Example for page 4.  Let  f(x) = sin(x).
          If       lim  x  = pi/4 ,          as   n->oo,  
verify that   lim  f(x ) = f(pi/4) ,   as   n->oo.

Clear[f,f0,X,X0,Y,y0,Y0,XY]
f[x_] = Sin[x]

X  = Table[N[Pi]/4 - 2.0^(-p),{p,1,15}];
Y  = Sin[X];
X0 = Transpose[{X,0 X}];
Y0 = Transpose[{0 Y,Y}];
XY = Transpose[{X,Y}];

Clear[grf,Lx0,Ly0,Lxy]
Lx0 = ListPlot[X0,Prolog->{PointSize[0.02]},
      DisplayFunction -> Identity];
Ly0 = ListPlot[Y0,Prolog->{PointSize[0.015]},
      DisplayFunction -> Identity];
Lxy = ListPlot[XY,Prolog->{PointSize[0.015]},
      DisplayFunction -> Identity];
grf = Plot[f[x],{x,0,0.8},PlotRange->{{0,0.8},{0,1.05}},
      DisplayFunction -> Identity];
Show[ grf,Lx0,Ly0,Lxy,
      PlotRange->{{0,0.8},{0,0.72}},
      Prolog->{PointSize[0.015]},
      Ticks->{Range[0,1,0.1],Range[0,1,0.1]},
      DisplayFunction -> $DisplayFunction ];

Print[{"Xk      ","Yk      "}];
For[i=1,i<=Length[X],i++,
    Print[XY[[i]]] ];
Print["The limiting values are:"];
Print[{N[Pi/4],N[f[Pi/4]]}];

y0 = Limit[ f[x], x->Pi/4 ];
f0 = f[Pi/4];
Print["lim f[xn] = ",y0,"    and    f[x0] = ",f0]
Print["Is f[x] continuous at x = x0 ?   ",y0 == f0]

Example 1.  Verify that    lim  sin(x) = sin(pi/6),   as   x -> pi/6.

Limit[Sin[x], x->Pi/6]

Sin[Pi/6]

Limit[Sin[x], x->Pi/6] == Sin[Pi/6]

Example 2.  Verify that    lim  sin(x)/x = 1,   as   x -> 0.

Limit[Sin[x]/x, x->0]

Example 3.  Verify that    lim  x exp(-x) = 0,   as   x -> oo

Limit[x Exp[-x], x->Infinity]

Theorem 1.2  (Intermediate Value Theorem).
     Assume that  f  is a continuous function on  [a,b],
and  L  is  any number between  f(a)  and  f(b).  Then
there exists a value  c  with  a < c < b  such that  f(c) = L.
Section 1.1,  Review of Calculus, Page 4

Example on page 4.  Consider the function  f(x) = cos(x-1),
which is continuous on the interval  [0,2.5]  and  L = 0.8
     Since  f(0) = 0.540302  <  L  <  1  =  f(1),
there exists a value  c1  with  0 < c1 < 1 so that  f(c1) = L.  
     Also, since  f(2.5) = 0.0707372  <  L  <  1  =  f(1),
there exists a value  c2  with  1 < c2 < 2.5  so that  f(c2) = L.

Clear[f]
f[x_] = Cos[x-1]

Clear[a1,a2,b1,b2,c,L]
a1 = 0;
b1 = 1;
a2 = 1;
b2 = 2.5;
L  = 0.8;
Print["f[x] = ",f[x]];
Print["   L = ",L];
Print["f[c] = ",L];
Solve[f[c] == L, c]

Clear[c1,c2,grf,ln0,ln1,ln2,y1,y2];
c1  = 1 - ArcCos[0.8];
y1  = Cos[c1-1];
c2  = 1 + ArcCos[0.8];
y2  = Cos[c2-1];
grf = Plot[f[x],{x,0,2.5},
      DisplayFunction -> Identity];
ln0 = Graphics[Line[{{0,L},{2.5,L}}]];
ln1 = Graphics[Line[{{c1,0},{c1,y1}}]];
ln2 = Graphics[Line[{{c2,0},{c2,y1}}]];
Show[ grf,ln0,ln1,ln2,
      DisplayFunction -> $DisplayFunction ];

Assume that  f  is a continuous on [a,b] and
L  is  any number between  f(a)  and  f(b).
There exists a value  c  with  a < c < b such that  f(c) = L.
Consider the two intervals and for this theorem.

Print["f(x) = ",f[x]];
Print["L = ",L];
Print["c = ",c1];
Print["f[c] = ",f[c1]];
Print["[a,b] = [",a1,",",b1,"]"];
Print[a1,"  <  ",c1"  <  ",b1];

Print["f(x) = ",f[x]];
Print["L = ",L];
Print["c = ",c2];
Print["f[c] = ",f[c2]];
Print["[a,b] = [",a2,",",b2,"]"];
Print[a2,"  <  ",c2"  <  ",b2];

Theorem 1.3  (Extreme Value Theorem 
                     for a Continuous Function).
Assume that  f  is a continuous function on the interval  [a,b].
Then there exists a lower bound M1 and an upper bound M2
and two numbers x1 and x2 in [a,b] such that
                M1  =  f(x1)  <=  f(x)  <=  f(x2)  =  M2     
whenever   a <= x <= b.
Section 1.1,  Review of Calculus, Page 5

Example for page 5.  Let  f(x) = 35 + 59.5x - 66.5x2  + 15x3.
Find the minimum and maximum of f(x) over the interval [0,3].

Solution 1.  Use Mathematica's built in function for finding
a minimum value.

Clear[f]
f[x_] = 35 + 59.5x - 66.5x^2 + 15x^3

Compute approximations for the minimum and maximum of f over [0,3]:

grf    = Plot[f[x],{x,0,3}];
minset = FindMinimum[ f[x],{x,2.5}];
maxset = FindMinimum[-f[x],{x,0.5}];
maxset[[1]] = - maxset[[1]];
Print["For the minimum we have:  ",minset];
Print[" "];
Print["For the maximum we have:  ",maxset];

Solution 2.  Since the given function f(x) is also differentiable,
                     locate the critical values where  f'(x) = 0.

f'[x];
solset = Solve[f'[x]==0,x];
Print["f'[x] = ",f'[x]];
Print["Solve f'[x] = 0"];
Print[solset];

xmin = solset[[2,1,2]];
ymin = f[xmin];
Print["For the minimum we have:"];
Print["xmin = ",xmin];
Print["ymin = ",ymin];

xmax = solset[[1,1,2]];
ymax = f[xmax];
Print["For the maximum we have:"];
Print["xmax = ",xmax];
Print["ymax = ",ymax];

epts = ListPlot[{{xmin,ymin},{xmax,ymax}},
       DisplayFunction -> Identity];
lmin = Graphics[Line[{{xmin+0.5,ymin},{xmin-0.5,ymin}}]];
lmax = Graphics[Line[{{xmax+0.5,ymax},{xmax-0.5,ymax}}]];
Show[  grf,epts,lmin,lmax,
       Prolog->{PointSize[0.02]},
       Ticks->{Range[0,3,0.5],Range[0,60,10]},
       DisplayFunction -> $DisplayFunction];
Print["The maximum value for the function f(x):"];
Print["f(x) = ",f[x]];
Print["x1 = ",xmax];
Print["max f(x) = f(",xmax,") = y1 = ",ymax];
Print[" "];
Print["The minimum value for the function f(x):"];
Print["f(x) = ",f[x]];
Print["x2 = ",xmin];
Print["min f(x) = f(",xmin,") = y2 = ",ymin];

Theorem 1.4  (Differentiable function implies 
                     continuous function).
    If  f(x)  is differentiable at  x = x0, 
then  f(x)  is continuous at  x = x0.
Section 1.1,  Review of Calculus, Page 5

Example for page 5.  Let  f(x) = sin(x+1).  Investigate the
derivative via limiting slopes of secant lines at  x0  = 0.

Clear[a,b,df,f]
f[x_]  = Sin[x+1];
df[x_] = Cos[x+1];
a = 0;
b = 1.2;
Print[" f[x] = ",f[x]];
Print["f'[x] = ",f'[x]];

Clear[ln0,lns,m,pt0,x0,x1,dx,y0,y1,yx]
ln0  = {};
pt0  = {};
For[j=1,j<=14,j++,
  x0 = N[a];
  y0 = N[f[x0]];
  x1 = N[a + 1/2^j];
  y1 = N[f[x1]];
  dx = x1-x0;
  dy = y1-y0;
  m  = dy/dx;
  ln0 = Append[ln0,Line[{{x0,y0},{b,y0+b*m}}]];
  pt0 = Append[pt0,{dx,dy,dy/dx}]; ];
  lns = Graphics[ln0];

grf = Plot[f[x],{x,a,b},
      DisplayFunction -> Identity];
Show[ grf,lns,
      PlotRange->{{0,1.2},{0,1.3}},
      DisplayFunction -> $DisplayFunction];
Print["Table of values for the slope."];
Print["  dx             dy             dy/dx"];
Print[TableForm[pt0]];
Print["The function f[x] = ",f[x]];
Print["is differentiable at x = ",a];

x00 = {};
y00 = {};
xy0 = {};
For[j = 1,j <= 14,j++,
  x1  = N[a + 1/2^j];
  y1  = N[f[x1]];
  x00 = Append[x00,{x1,0}];
  y00 = Append[y00,{0,y1}];
  xy0 = Append[xy0,{x1,y1}]; ];
dot0  = ListPlot[xy0,
        DisplayFunction -> Identity];
dotx  = ListPlot[x00,
        DisplayFunction -> Identity];
doty  = ListPlot[y00,
        DisplayFunction -> Identity];
Show[   grf,dot0,dotx,doty,
        PlotRange->{{0,0.5},{0,1.0}},
        Prolog->{PointSize[0.015]},
        DisplayFunction -> $DisplayFunction];
Print["  xn             yn"];
Print[TableForm[xy0]];
Print["The limiting values are:"];
Print[x0,"\t \t \t \t",y0];
Print["The function f[x] = ",f[x]];
Print["is continuous at x = ",a];

Theorem 1.5  (Rolle's Theorem).
    Assume that  f  is continuous on the interval  [a,b],
and that  f'(x)  exists for all  a < x < b.
     If  f(a) = f(b) = 0,  then there exists a value  c,
with  a < c < b,  such that  f'(c) = 0.
Section 1.1,  Review of Calculus, Page 5

Example for page 5.  The function  f(x) = sin(x) - 1/2
has two zeros in the interval [0,3].  
According to Rolle's theorem, there exists a 
value  c  with 0 < c < 3 such that  f'(c) = 0.
Enter the function  f  and find its zeros.

Clear[f,x]
f[x_] = Sin[x] - 0.5;
Print["Solve  ",f[x]," = 0"];
Solve[ f[x] == 0,x]

Or, two roots can obtained using the inverse function.

r1 = ArcSin[0.5]

r2 = r1 + 2 (N[Pi]/2 - r1)

Find the function  f'  and its zeros.

f'[x];
solset = Solve[ f'[x] == 0, x];
x1 = solset[[1,1,2]];
y1 = f[x1];
Print["         f[x] = ",f[x]];
Print["        f'[x] = ",f'[x]];
Print["Solve  ",f'[x]," = 0"];
Print["   x1  = ",x1];
Print[" f[x1] = ",y1];
Print["f'[x1] = ",f'[x1]];

grf = Plot[f[x],{x,0,3},
      DisplayFunction -> Identity];
cpt = ListPlot[{{x1,y1}},
      DisplayFunction -> Identity];
hor = Graphics[Line[{{x1+0.5,y1},{x1-0.5,y1}}]];
Show[ grf,cpt,hor,
      Prolog->{PointSize[0.02]},
      Ticks->{Range[0,3,0.5],Range[-1,1,0.5]},
      DisplayFunction -> $DisplayFunction ];
Print["The value c for Rolle's theorem for the function"];
Print["f(x) = ",f[x]];
Print["c = ",x1];
Print["f(c) = f(",x1,") = y1 = ",f[x1]];
Print["f'(c) = ",f'[x1]];

Theorem 1.6  (Mean Value Theorem).
    Assume that  f  is continuous on the interval  [a,b],
and that  f'(x)  exists for all  a < x < b.
     Then there exists a number  c,  with  a < c < b,  
such that    f'(c)  =  [f(b) - f(a)]/[b - a]  =  m.
Section 1.1,  Review of Calculus, Page 5

Example on page 6.  Consider f(x) = sin(x) over [0.1, 2.1].
According to the Mean Value theorem, there exists a 
value  c  with 0.1 < c < 2.1 such that
            f'(c)  = [f(2.1) - f(0.1)]/[2.1 - 0.1]  =  m.

Enter the function  f  and compute m:

Clear[a,b,f]
f[x_]  = Sin[x];
a  = 0.1;
b  = 2.1;
m = ( f[b] - a )/( b -  a );
Print["f[x] = ",f[x]];
Print["The interval is  [",a,", ",b,"]."];
Print["The slope is  m = ",m];

Use the function  f'  and where f'(c) = m:

solset = Solve[f'[x] == m, x];
c0 = solset[[1,1,2]];
y0 = f[c0];
ya = f[a];
yb = f[b];
dots  = Graphics[{Point[{a,f[a]}],
                 Point[{b,f[b]}],
                 Point[{c0,f[c0]}]}];
lines = Graphics[{Line[{{a,f[a]},{b,f[b]}}],
        Line[{{c0-1/2,f[c0]-m/2},{c0+1/2,f[c0]+m/2}}]}];
grf   = Plot[f[x],{x,0.1,2.1},
        DisplayFunction -> Identity];
Show[grf,dots,lines,
     PlotRange->{{0,2.2},{0,1.12}},
     Prolog->{PointSize[0.02]},
     Ticks->{Range[0,3,0.5],Range[0,1,0.5]},
     DisplayFunction -> $DisplayFunction ];
Print["The interval is [",a," , ",b,"]"];
Print["f`(x) = ",f'[x]];
Print["The value c where  f`(c) = m = ",m,"  is:"];
Print["c = ",c0];

Theorem 1.7  (Extreme Value Theorem 
                     for a Differentiable Function).
     Assume that  f  is continuous on the interval [a,b] and that
f'(x)  exists for all  a < x < b.  Then there exists a lower bound
M1 and an upper bound M2 and values x1, x2 in [a,b]  such that
     M1 = f(x1) <= f(x) <= f(x2) = M2     whenever   a <= x <= b.
The numbers  x1  and  x2  are critical points of  f  and
occur at either at the endpoints of  [a,b]  or where  f'(x) = 0.
Section 1.1,  Review of Calculus, Page 6

Example on page 5.  Let  f(x) = 35 + 59.5x - 66.5x2  + 15x3.
Find the minimum and maximum of f(x) over the interval [0,3].

f[x_] = 35 + 59.5x - 66.5x^2 + 15x^3

grf = Plot[f[x],{x,0,3}];

Solution.  Since the given function f(x) is also differentiable,
                locate the critical values where  f'(x) = 0.

f'[x]

The quadratic equation can be used to solve f'(x) = 0.

x1 = (133.0 - Sqrt[133.0^2 - 4*45*59.5])/(2*45);
x2 = (133.0 + Sqrt[133.0^2 - 4*45*59.5])/(2*45);
x1
f'[x1]
x2
f'[x2]

Or Mathematica's can be used to solve f'(x) = 0.

Print["Solve  ",f'[x]," = 0"];
solset = Solve[f'[x]==0,x]
xmin = solset[[2,1,2]]
f'[xmin]
ymin = f[xmin]
xmax = solset[[1,1,2]]
f'[xmax]
ymax = f[xmax]

grf  = Plot[f[x],{x,0,3},
       DisplayFunction -> Identity];
epts = ListPlot[{{xmin,ymin},{xmax,ymax}},
       DisplayFunction -> Identity];
lmin = Graphics[Line[{{xmin+0.5,ymin},{xmin-0.5,ymin}}]];
lmax = Graphics[Line[{{xmax+0.5,ymax},{xmax-0.5,ymax}}]];
Show[  grf,epts,lmin,lmax,
       Prolog->{PointSize[0.02]},
       Ticks->{Range[0,3,0.5],Range[0,60,10]},
       DisplayFunction -> $DisplayFunction ];
Print["The maximum value for the function f(x):"];
Print["f(x) = ",f[x]];
Print["x1 = ",xmax];
Print["max f(x) = f(",xmax,") = y1 = ",ymax];
Print[" "];
Print["The minimum value for the function f(x):"];
Print["f(x) = ",f[x]];
Print["x2 = ",xmin];
Print["min f(x) = f(",xmin,") = y2 = ",ymin];

Theorem 1.8  (Generalized Rolle's Theorem).
     Assume that  f  is continuous on the interval [a,b], and
that the derivatives  f'(x),  f''(x), ... ,f(n)(x)  exist on (a,b),
and that  x0 , x1, ... , xn  all lie in  [a,b].
     If  f(x0) = 0  for  j=0,1,...,n  then there exists a
value  c,  with  a < c < b  such that  f(n)(c) = 0.
Section 1.1,  Review of Calculus, Page 6

Example for Page 6.  The function  f(x) = x^3-7*x^2+14*x-8
has three zeros in the interval [0,5].  
According to the generalized Rolle's theorem, there exists a 
value  c  with 0 < c < 5 such that  f''(c) = 0.

Enter the function  f  and find its zeros.

Clear[a,b,c0,c1,c2,f,r1,r2,r3]
a=0;  
b = 5;
f[x_]  = x^3 - 7x^2 + 14x - 8;
solset = Solve[ f[x] == 0, x];
Print["Solve  ",f[x]," = 0"];
r1 = solset[[1,1,2]]
r2 = solset[[2,1,2]]
r3 = solset[[3,1,2]]

Use the function  f'  and find its zeros.
Mathematica can be used to find the roots.
These values are entered into c1 and c2.

f'[x];
Print["Solve  ",f'[x]," = 0"];
solset1 = Solve[ f'[x] == 0, x]
c1 = N[solset1[[1,1,2]]]
c2 = N[solset1[[2,1,2]]]

Use the function  f''  and find its root(s).

f''[x];
Print["Solve  ",f''[x]," = 0"];
solset2 = Solve[ f''[x] == 0, x]
c0 = N[solset2[[1,1,2]]]

grf0 = Plot[f[x],{x,a,b},
       DisplayFunction -> Identity];
dots = Graphics[{Point[{r1,f[r1]}],
                 Point[{r2,f[r2]}],
                 Point[{r3,f[r3]}]}];
Show[  grf0,dots,
       PlotRange->{{0,5},{-5,5}},
       Prolog->{PointSize[0.02]},
       Ticks->{Range[0,5,1],Range[-5,5,1]},
       DisplayFunction -> $DisplayFunction ];
Print["The interval is [",a," , ",b,"]"];
Print["f(x) = ",f[x]];
Print["There are 3 roots in this interval."];
Print[TableForm[{r1,r2,r3}]];

grf1 = Plot[f'[x],{x,a,b},
       DisplayFunction -> Identity];
dots = Graphics[{Point[{c1,f'[c1]}],
                 Point[{c2,f'[c2]}]}];
Show[  grf1,dots,
       PlotRange->{{0,5},{-3,9}},
       Prolog->{PointSize[0.02]},
       Ticks->{Range[0,5,1],Range[-4,15,2]},
       DisplayFunction -> $DisplayFunction ];
Print["The interval is [",a," , ",b,"]"];
Print["f'(x) = ",f'[x]];
Print["There are 2 roots in this interval."];
Print[TableForm[{c1,c2}]];

grf2 = Plot[f''[x],{x,a,b},
       DisplayFunction -> Identity];
dots = Graphics[{Point[{c0,f''[c0]}]}];
Show[  grf2,dots,
       PlotRange->{{0,5},{-9,9}},
       Prolog->{PointSize[0.02]},
       Ticks->{Range[0,5,1],Range[-14,15,2]},
       DisplayFunction -> $DisplayFunction ];
Print["The interval is [",a," , ",b,"]"];
Print["f''(x) = ",f''[x]];
Print["There is 1 root in this interval."];
Print[TableForm[{c0}]];

Theorem 1.9  (First Fundamental Theorem).
     If  f  is continuous over the interval [a,b], then there 
exists a function  F, called the antiderivative of  f, such that
b
º  f(x) dx   =   F(b) - F(a)     where  F'(x) = f(x).
a
Section 1.1,  Review of Calculus, Page 7

Example for Page 7.  Let f(x) = tan(x).  Show that
  p/4
   º     f(x) dx   =   F(pi/4) - F(-pi/3)     where  F'(x) = f(x).
-p/3

Clear[f,F]
f[x_] = Tan[x];
F[x_] = Integrate[f[x],x];
Print["f[x] = ",f[x]];
Print["The integral is:"];
Print["F[x] = ",F[x]];

a  = -Pi/3;
b  =  Pi/4;
Print["f[x] = ",f[x]];
Print["F[x] = ",F[x]];
Print["The interval is  [",a,",",b,"]"];
Print["The value of the integral is:"];
Print["F[b] - F[a] = ",F[b] - F[a],"  =  ",N[F[b] - F[a]]];

Or use Mathematica's integrate command with limits of integration:

a  = -Pi/3;
b  =  Pi/4;
Print["f[x] = ",f[x]];
v  = Integrate[f[x],{x,a,b}];
Print["The interval is  [",a,",",b,"]"];
Print["The value of the integral is:"];
Print[v,"  =  ",N[v]];

grf  =  Plot[f[x],{x,a,b},
        Ticks->{Range[-1,1,0.5],Range[-2,1,0.5]},
        DisplayFunction -> Identity];
lines = Graphics[{Line[{{a,0},{a,f[a]}}],
                  Line[{{b,0},{b,f[b]}}]}];
Show[   grf,lines,
        DisplayFunction -> $DisplayFunction ];
Print["The function is f[x] = ",f[x]];
Print["The value of the integral is over [",
       a,",",b,"] is:  ",N[v]];

Theorem 1.10  (Second Fundamental Theorem).
     If  f  is continuous over the interval  [a,b],  then
     d     x
     --    º  f(t) dt  =  f(x),     where  a < x < b.
     dx   a
Section 1.1,   Review of Calculus, Page 8

Example for Page 7.  Let   f(x)  =  2 x2  ex sin(x).   Show that
     d     x
     --    º  f(t) dt  =  f(x).
     dx   0

Clear[f,g,t,x]
f[t_] = 2 t^2 Exp[t] Sin[t];
Print["f[t] = ",f[t]];
g[x_] = Integrate[f[t],{t,0,x}];
Print["g[x] = ",g[x]];
Print["g'[x] = ",g'[x]];
Print["g'[x] = ",Expand[g'[x]]];
Print["f[x] = ",f[x]];

Example for Page 7.  Let  f(x) = 1/sqrt(2 pi) exp(-x2 /2).  Show that
      d     x
     --    º  f(t) dt  =  f(x).
     dx   0

Clear[f,g,t,x]
f[t_] = Exp[-t^2/2]/Sqrt[2 Pi];
Print["f[t] = ",f[t]];
g[x_] = Integrate[f[t],{t,0,x}];
Print["g[x] = ",g[x]];
Print["g'[x] = ",g'[x]];
Print["f[x] = ",f[x]];

Plot[f[x],{x,-4,4}];
Print["Graph of f[x] = ",f[x]];

Plot[g[x] + 1/2,{x,-4,4}];
Print["Graph of g[x] + 1/2 = ",g[x] + 1/2];

Theorem 1.11  (Mean Value Theorem for Integrals).
     If  f  is continuous over the interval  [a,b],  then 
there exists a number  c,  with  a < c < b,  such that
          1    b
       -----  º f(x) dx  =  f(c).
       b - a  a
Section 1.1,   Review of Calculus, Page 8

Example on Page 7.  Let  f(x) = sin(x) + sin(3x)/3,  find  c 
                   1    2.5
 so that     -----  º  (sin(x) + sin(3x)/3) dx  =  f(c).
                  2.5   0

Clear[a,b,c,f,g,h,L,v]
f[x_] = Sin[x] + 1/3 Sin[3 x]

g[x_] = Integrate[f[x],x]

a  = 0;
b  = 2.5;
v  = Integrate[f[x],{x,a,b}]
L  = 1/(2.5-0)*Integrate[f[x],{x,a,b}]
solset = FindRoot[f[x] == L, {x,1}]
c  = solset[[1,2]]
h  = f[c]

grf   = Plot[f[x],{x,a,b},
        Ticks->{Range[0,2.5,0.5],Range[0,1,0.2]},
        DisplayFunction -> Identity];
lines = Graphics[{Line[{{a,0},{a,f[a]}}],
                  Line[{{b,0},{b,f[b]}}]}];
line2 = Graphics[{Line[{{a,0},{a,h}}],
                  Line[{{a,h},{b,h}}],
                  Line[{{b,h},{b,0}}]}];
Show[   grf,lines,line2,
        DisplayFunction -> $DisplayFunction ];
Print["The function is f[x] = ",f[x]];
Print["The value of the integral is over [",
       a,",",b,"] is:  ",N[v]];
Print["The area under the curve is the same as"];
Print["the area of a rectangle of height h = ",h];

Theorem 1.12  (Weighted Integral Mean Value Theorem).

     Assume that both  f  and   g  are continuous over [a,b],
 and that  0 <= g(x)   for a <= x <= b.
     Then there exists a number c, with a < c < b, such that
        b                               b
        º f(x)g(x) dx   =   f(c) º g(x) dx.
        a                               a
Section 1.1,   Review of Calculus, Page 8

Example for age 8.  Let  f(x) = x^2  and  g(x) = 1 + sin(x)

and let [a, b] = [0, 3p/2].  find the value  c  so that
      3p/2                             3p/2
        º   f(x)g(x) dx   =   f(c)  º    g(x) dx.
       0                                  0

Clear[f,g,h]
f[x_] = x^2
g[x_] = 1 + Sin[x]
h[x_] = f[x] g[x]

L  = Integrate[h[x],{x,0,3Pi/2}] //Expand
R  = Integrate[g[x],{x,0,3Pi/2}]
fc = N[L/R]
c  = Sqrt[fc]

Print["   3¹/2                   3¹/2"];
Print["    º  f(x)g(x)dx  =  f(c)  º  g(x)dx"];
Print["    0                      0"];
Print["The value  c  and  f(c)  are given by:"];
Print["Computations for the left and right side are:"];
Print["c = ",c];
Print["f(c) =  ",f[c]];
Print["The left  side is:  ", N[L]];
Print["The right side is:  ", f[c] N[R]];

Theorem 1.13  (Taylor's Theorem).
     Assume that  f  and its derivatives f'(x),f''(x),...,f(n)(x)
are all continuous on [a,b].
     Suppose that x0  is a fixed value in the interval  [a,b].
Then       f(x)  =  Pn(x)  +  Rn(x),
where  Pn(x)  is the n-th degree Taylor polynomial expanded
about  x0  and  Rn(x)  is the Lagrange form of the remainder.
Section 1.1,   Review of Calculus, Page 10

Example on Page 10.  Let  f(x) = sin(x).
Find the Taylor polynomial of degree n = 9 for
f(x)  expanded about the point  x0  = 0.

Clear[f,n,Pn,x0]
n = 9;
x0 = 0;
f[x_] = Sin[x]
Pn[x_] = Series[f[x],{x,x0,n}] //Normal

Plot[{f[x],Pn[x]}, {x,0,2Pi}];
Print[" f[x] = ",f[x]];
Print["Pn[x] = ",Pn[x]];

Theorem 1.14  (Horner's Method for Polynomial Evaluation).
     Assume that  P(x) = anxn + an-1xn-1 +...+ a2x2 + a1x  + a0.
 Then P(x) can be computed recursively as follows:
 Set       bn  =  an , 
 and      bk  =  ak + x bk+1     for  k = n-1, n-2,...,1, 0.
 Then   P(x) = b0.
Section 1.1,   Review of Calculus, Page 12

Example on Page 13.  Use Horner's method to evaluate  P(3)
where  P(x) = x^5 - 6x^4 + 8x^3 + 8x^2 + 4x - 40

First, find Horner's nested form of the polynomial.

Clear[A,B,H,k,n,P,Q]
P[x_] = x^5 - 6 x^4 + 8 x^3 + 8 x^2 + 4 x - 40;
Q[x_] = - 40 + (4 + (8 + (8 + (-6 + x) x) x) x) x;
A = CoefficientList[P[x],x];
n = Length[A];
H[x_] := Module[{},
         B = Table[0,{n}];
         B[[n]] = A[[n]];
         For[k=n-1,1<=k,k--,
           B[[k]] =  A[[k]] + x B[[k+1]]];
         Return[B[[1]]] ];

P[3]
Q[3]
H[3]

Consider the computational time for vector arguments
for evaluating  P(x)  and the nested form  Q(x).

Clear[i,V]
V = Table[i,{i,1,100}];
Ptimeset = Timing[P[V]];
Qtimeset = Timing[Q[V]];
Ptimeset[[1]]
Qtimeset[[1]]

Remark. It might be interesting to perform the computation
              several times and see the different results.

V = Table[i,{i,1,200}];
Ptimeset = Timing[P[V]];
Qtimeset = Timing[Q[V]];
Ptimeset[[1]]
Qtimeset[[1]]

V = Table[i,{i,1,400}];
Ptimeset = Timing[P[V]];
Qtimeset = Timing[Q[V]];
Ptimeset[[1]]
Qtimeset[[1]]

V = Table[i,{i,1,800}];
Ptimeset = Timing[P[V]];
Qtimeset = Timing[Q[V]];
Ptimeset[[1]]
Qtimeset[[1]]

Theorem 1.15  (Geometric Series).
     If  |r| < 1, then
     1 + r + r2 + ... + rn  + ...   =   1/(1 - r).
Section 1.2,   Binary Numbers, Page 18

Example 1.2, Page 19.  Sum the infinite series

     1 + 1/4 + 1/16 + 1/64 + 1/256 + ... + 1/4^n + ...

First load the infinite series package.
( Make sure you do this only once! )

<< Algebra`SymbolicSum`

Clear[n,S,V]
V = {};
S = 0;
For[n = 0,n <= 12,n++,
    S = S + N[1/4^n];
    V = Append[V,{n, S}] ];
dots  = ListPlot[V,
        DisplayFunction -> Identity ];
line  = Graphics[{Line[{{0,4/3},{12,4/3}}]}];
Show[   dots,line,
        PlotRange->{{0,12},{0,1.35}},
        Prolog->{PointSize[0.015]},
        Ticks->{Range[0,15,1],Range[0,1.5,0.25]},
        DisplayFunction -> $DisplayFunction ];
Print["Partial sums for the series · 1/4^n."];

Mathematica  can find the sum of the infinite series.

Clear[n]
Sum[ 1/4^n, {n,0,Infinity}]

Example 1.3, Page 21.  Find the sum of the binary number
     S = 0.1011001100110...
                                          two
 Solution:
     S = 1/2 + 1/8*(1 + 1/16 + ... + 1/16n + ... )
                 + 1/16*(1 + 1/16 + ... + 1/16n + ... )
     S = 1/2 + 3/16*(1 + 1/16 + ... + 1/16n + ... )

Clear[n]
1/2 + 3/16 Sum[1/16^n, {n,0,Infinity}]

Theorem 1.16  (Big "O" remainders for Taylor's Theorem).
     Assume that  f(h)  =  p(h)  +  O(hn),  
                   and   g(h)  =  q(h)  +  O(hm),
                   and       p  =  min{m,n}.
Then
         f(h) + g(h)  =  p(h) + q(h)   +  O(hp),
         f(h) * g(h)  =  p(h) * q(h)   +  O(hp),
         f(h) / g(h)  =  p(h) / q(h)   +  O(hp).
Section 1.2,   Binary Numbers, Page 33

Example 1.11, Page 34.

First find Taylor expansions for exp(h) and cos(h):

Clear[m,n,p,s1,s2]
n = 5;
m = 7;
p = Min[n,m];
s1 = Series[Exp[h],{h,0,n-1}]
s2 = Series[Cos[h],{h,0,m-1}]

Then find the series for their sum.

s1 + s2
Series[ Exp[h] + Cos[h], {h,0,p-1}]

Then find the series for their product.

s1 s2
Series[ Exp[h] * Cos[h], {h,0,p-1}]

Then find the series for their quotient.

s1 / s2
Series[ Exp[h] / Cos[h], {h,0,p-1}]

Theorem 1.17  (Remainder term for Taylor's Theorem).
     Assume that  f(x)  and its derivatives  f'(x), f''(x) ,..., f(n)(x)
are all continuous on [a,b].  If both  x0   and  x = x0  + h
lie in the interval [a,b],  then
                f(x)  =  Pn(x)  +  O(hn),  where
 Pn(x)  is the n-th degree Taylor polynomial expanded about x0.
Section 1.2,   Error Analysis, Page 34

Example for Page 34.  Consider  f(x) = sin(x)  and the Taylor
polynomials of degree  n = 3, 5 and 7  expanded about x0 = 0:
n-th degree Taylor polynomial                    Order of the approximation

 P3(x)  =  x - x3/6                                         x5/120  =  O(x5)
 
 P5(x)  =  x - x3/6 + x5/120                          x7/5040  =  O(x7)
 
 P7(x)  =  x - x3/6 + x5/120 - x7/5040         x9/362880  =  O(x9)

(i)   Show that   |sin(x) - P3(x)|  <=  |x5/120|   over  [-1,1].
      REMARK. The factor 1.04 was added to separate the curves.

Clear[f,P3,E3];
f[x_]  = Sin[x];
P3[x_] = Series[f[x], {x,0,3}]  //Normal;
E3[x_] = x^5/5!;
maxerr = Abs[f[1.0]-P3[1.0]];
Plot[{Abs[f[x]-P3[x]], 1.04 Abs[E3[x]]},{x,-1,1},
      PlotRange->{{-1,1},{0,maxerr}}];
Print["f(x)  = ",f[x]];
Print["P3(x) = ",P3[x]];
Print["|sin(x) - P3(x)|  <=  |x^5/120|  over [-1,1]."];
Print["The maximum error over [-1,1] is:  ",maxerr];

(ii)  Show that   |sin(x) - P5(x)| <= |x^7/5040|   over  [-1,1].
       REMARK. The factor 1.06 was added to separate the curves.

Clear[f,P5,E5]
f[x_]  = Sin[x];
P5[x_] = Series[f[x], {x,0,5}]  //Normal;
E5[x_] = x^7/7!;
maxerr = Abs[f[1.0]-P5[1.0]];
Plot[{Abs[f[x]-P5[x]], 1.06 Abs[E5[x]]},{x,-1,1},
      PlotRange->{{-1,1},{0,maxerr}}];
Print["f(x)  = ",f[x]];
Print["P5(x) = ",P5[x]];
Print["|sin(x) - P5(x)|  <=  |x^7/5040|  over [-1,1]."];
Print["The maximum error over [-1,1] is:  ",maxerr];

(iii) Show that   |sin(x) - P7(x)|  <=  |x^9/362880|   over  [-1,1].
       REMARK. The factor 1.08 was added to separate the curves.

Clear[f,P7,E7]
f[x_]  = Sin[x];
P7[x_] = Series[f[x], {x,0,7}]  //Normal;
E7[x_] = x^9/9!;
maxerr = Abs[f[1.0]-P7[1.0]];
Plot[{Abs[f[x]-P7[x]], 1.08 Abs[E7[x]]},{x,-1,1},
      PlotRange->{{-1,1},{0,maxerr}}];
Print["f(x)  = ",f[x]];
Print["P7(x) = ",P7[x]];
Print["|sin(x) - P7(x)|  <=  |x^9/362880|  over [-1,1]."];
Print["The maximum error over [-1,1] is:  ",maxerr];