Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

CONTENTS

Chapter 5.   Curve Fitting
	Algorithm 5.1  Least Squares Line.
	Algorithm 5.2  Least Squares Polynomial.
	Algorithm 5.3  Non-linear Curve Fitting.
	Algorithm 5.4  Cubic Splines.
	Algorithm 5.5  Trigonometric Polynomials.

Clear["Global`*"];
Off[General::"spell1"];
Clear[a, a0, b, b0, c, c0, ChangeVariables, Coeffs, ComputeCoeff, CS, 
        Ctype, CurveFit, d, Differences, f, f0, h, i, j, k, LeastSqPoly, m, 
        mn, mn0, MakeCoeffs, MakePoints, MakeSpline, MakeTP, Modify, n, sum, 
        pow, prod, Regression, Regression2, RMS, SelectFunction, Stype, SumX, 
        SumXY, t, t0, TP, TriDiagonal, v0, x, xy, X, X0, X1, XY, XY0, Xmean, 
        y, Y, Y0, Y1, Ymean];
On[General::"spell1"];

Algorithm 5.1 Least Squares Line.  To construct the least squaresline  y = Ax + B  that fits the  n  data points  ((Subscript[x, 1]), (Subscript[y, 1])), ((Subscript[x, 2]), (Subscript[y, 2])), ⋯, 
      ((Subscript[x, n]), (Subscript[y, n])).
Remark.  The algorithm is computationally stable.  It gives reliable results in cases when the normal equations are ill-conditioned.
Section 5.1   Least-Squares Line   Page  264

Execute this cell to activate Least Squares Line.
Regression[XY0_]
XY0_ is the set of points.

Off[General::"spell1"]; Clear[Regression, XY0]; \(Regression[XY0_] := 
      Module[{a, b, n, SumX, SumXY, X, Y, XY, Xmean, Ymean}, 
\ \ \ 
        XY = XY0; 
\ \ \ n = Length[XY]; 
\ \ \ 
        X = \(Transpose[XY]\)[[1]]; 

        \ \ \ Y = 
          \(Transpose[XY]\)[[2]]; 

        \ \ \ Xmean = 0; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), 
          Xmean = Xmean + X[[k]]]; 

        \ \ \ Xmean = Xmean\/n; 
\ \ \ Ymean = 0; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), 
          Ymean = Ymean + Y[[k]]]; 

        \ \ \ Ymean = Ymean\/n; 
\ \ \ SumX = 0; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), \ 
          SumX = SumX + 
              \((X[[k]] - Xmean)\)\^2]; 
        
\ \ \ SumXY = 0; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), 
\ \ \ \ \ \ 
          SumXY = SumXY + 
              \((X[[k]] - Xmean)\)\ 
                \((Y[[k]] - Ymean)\)]; 

        \ \ \ a = SumXY\/SumX; 
\ \ \ b = Ymean - a\ Xmean; 
\ \ \ E2 = 0; 
        
\ \ \ For[k = 1, k ≤ n, \(k++\), 
\ \ \ \ \ \ 
          Z = a\ X[[k]] + b; 
          Err = Y[[k]] - Z; 
          E2 = N[E2 + Err\^2]]; 
\ \ \ E2 = \@\(E2\/n\); 
\ \ \ 
        Return[a\ x + b]; ]; 
\(On[General::"spell1"]; \)\)

Chapter 5  Example 5.2,  Page 261.  Find the "least squares" line for the data points:
(-1, 10), (0,  9), (1, 7), (2, 5), (3, 4), (4, 3), (5, 0), (6,  -1).

\(Clear[f, x, xy]; 

    xy = {{\(-1\), 10}, {0, 9}, {1, 7}, {2, 5}, {3, 4}, {4, 3}, {5, 0}, {6, 
          \(-1\)}}; 
f[x_] = N[Regression[xy]]; 

    Print["data points = ", xy]; 
Print["f[x] = ", f[x]]; \)

dots = 
      ListPlot[xy, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; \(gr = 
      Plot[f[x], {x, \(-1\), 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The least squares line: y = ", f[x]]; \)

The RMS error E2 for this example is:

\(\(Print[\*""The RMS error \!\(E\_2\) is: "", \ E2]; \)\)

Compare with the built-in Mathematica subroutine for finding the least squares line.

\(Print["f[x] = ", f[x]]; 
Print["f[x] = ", Fit[xy, {1, x}, x]]; 
    \)

Algorithm 5.2 Least Squares Polynomial.  To construct the "least squares" polynomial of degree  m  of the form:
(Subscript[P, m]) (x)  =  
        (Subscript[c, 1]) + (Subscript[c, 2]) x + (Subscript[c, 3]) (x^2) + (Subscript[c, 4]) (x^3) +  
          (Subscript[c, m]) (x^(m - 1)) +  (Subscript[c, m + 1]) (x^m)
that fits the  n  data points  ((Subscript[x, 1]), (Subscript[y, 1])), ((Subscript[x, 2]), (Subscript[y, 2])), ⋯, 
      ((Subscript[x, n]), (Subscript[y, n])).
Section  5.2   Curve Fitting   Page  278

Execute this cell to activate Least Squares Polynomial.
LeastSqPoly[XY0_,m_]
XY0_ is the set of points,
m_ is the degree.

Off[General::"spell1"]; Clear[LeastSqPoly, XY0, m, x]; \(LeastSqPoly[XY0_, m_] := 
      Module[{a, b, c, i, j, k, sum, pow, prod, X, Y, XY}, 
XY = XY0; 

        n = Length[XY]; 
a = Table[0, {m + 1}, {m + 1}]; 

        b = Table[0, {m + 1}]; 
c = Table[0, {m + 1}]; 
        pow = Table[0, {2\ m + 1}]; 

        X = \(Transpose[XY]\)[[1]]; 
        Y = \(Transpose[XY]\)[[2]]; 

        For[k = 1, k ≤ n, \(k++\), 
\ \ \ 
          YK = Y[[k]]; 
\ \ \ 
          XK = X[[k]]; 
\ \ \ 
          prod := 1; 
\ \ \ 
          For[r = 1, r ≤ m + 1, \(r++\), 
\ \ \ \ \ \ 
            b[[r]] = 
              b[[r]] + YK\ prod; 

            \ \ \ \ \ \ prod = prod\ XK]; \ ]; 

        For[j = 1, j ≤ 2\ m, \(j++\), 
\ \ \ 
          pow[[j + 1]] = 0]; 
\ \ \ 
        pow[[0 + 1]] := n; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), 
          XK = X[[k]]; 
\ \ \ \ \ \ 
          prod = X[[k]]; 
\ \ \ \ \ \ 
          For[j = 1, j ≤ 2\ m, \(j++\), 
\ \ \ \ \ \ \ \ \ 
            pow[[j + 1]] = 
              pow[[j + 1]] + prod; 

            \ \ \ \ \ \ \ \ \ prod = prod\ XK]; \ ]; 

        For[i = 1, i ≤ m + 1, \(i++\), 
\ \ \ 
          \(For[j = 1, j ≤ m + 1, \(j++\), 
\ \ \ \ \ \ 
            \(a[[i, j]] = 
              pow[[i + j - 1]]; \)\ ]; 
          \)\ ]; 
c = LinearSolve[a, b]; 

        sum = c[[m + 1]]; 

        For[k = m, 1 ≤ k, \(k--\), 
          sum = c[[k]] + sum\ x]; 

        sum = Expand[sum]; 
Return[sum]; ]; 
\(On[General::"spell1"]; 
    \)\)

Chapter 5  Example 5.6,  Page 277.  Find the "least squares" parabola for the data:
(-3, 3), (0, 1), (2, 1), (4, 3).

\(xy = {{\(-3\), 3}, {0, 1}, {2, 1}, {4, 3}}; 

    P[x_] = LeastSqPoly[xy, 2]; 
Print["data points = ", xy]; 

    Print["P[x] = ", P[x]]; 
Print["P[x] = ", N[P[x]]]; \)

\(dots = 
      ListPlot[xy, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[P[x], {x, \(-3\), 4}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, AxesLabel -> {"x", "y"}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The least squares polynomial is:"]; 

    Print["P[x] = ", P[x]]; 
Print["P[x] = ", N[P[x]]]; \)

Compare with the built-in Mathematica subroutine for finding the least squares polynomial.  
Remark. Mathematica returns only decimal approximations for the coefficients.

\(Print["P[x] = ", P[x]]; 
Print["P[x] = ", N[P[x]]]; 

    Print["P[x] = ", Fit[xy, {1, x, x\^2}, x]]; \)

Algorithm 5.3 Non-linear Curve Fitting.  To construct a non-linear fit  y = f(x)
given the  n  data points  ((Subscript[x, 1]), (Subscript[y, 1])), ((Subscript[x, 2]), (Subscript[y, 2])), ⋯, 
      ((Subscript[x, n]), (Subscript[y, n])).
 by using "data linearization."   Available curves are:
f[x] = a x + b,  f[x]=(a/x)+b,               f[x]=(d/(x + c)),         f[x]=(1/(a x + b)),
f[x]=(x/(a + b x)),  f[x] = a Log[x] + b,  f[x] = c (e^(a x)),
f[x] = c (x^a),       f[x]=(1/((((a x + b))^2))),       f[x] = c x (e^(-d x)),  f[x]=(L/(1 + c (e^(a x)))).
Section  5.2   Curve Fitting   Page  280

Execute this cell to activate Non-linear Curve Fitting.
CurveFit[Ctype_]
Ctype_ is the curve type.
Remark.  The data must be stored in  XY  and the limiting value  L  must be given.

Off[General::"spell1"]; Clear[SelectFunction, f, Ctype]; \(\(SelectFunction[Ctype_] := 
      Module[{}, 

        Which[
\ \ \ Ctype == 1, f[x_] = a\ x + b, 
\ \ \ Ctype == 2, 
          f[x_] = a\/x + b, 
\ \ \ Ctype == 3, f[x_] = d\/\(x + c\), 
\ \ \ 
          Ctype == 4, f[x_] = 1\/\(a\ x + b\), 
\ \ \ Ctype == 5, 
          f[x_] = x\/\(a + b\ x\), 
\ \ \ Ctype == 6, f[x_] = a\ Log[x] + b, 
          
\ \ \ Ctype == 7, f[x_] = c\ E\^\(a\ x\), 
\ \ \ 
          Ctype == 8, f[x_] = c\ x\^a, 
\ \ \ Ctype == 9, 
          f[x_] = 1\/\((a\ x + b)\)\^2, 
\ \ \ Ctype == 10, 
          f[x_] = c\ x\ E\^\(\(-d\)\ x\), 
\ \ \ Ctype == 11, 
          f[x_] = L\/\(1 + c\ E\^\(a\ x\)\)]]; \)\)

Clear[ChangeVariables, X1, Y1, Ctype]; \(\(ChangeVariables[X1_, Y1_, Ctype_] := 
      Module[{}, 

        Which[
\ \ \ Ctype == 1, X = X1; \(Y = Y1; \), 
\ \ \ Ctype == 2, 
          X = 1\/X1; \(Y = Y1; \), 
\ \ \ Ctype == 3, X = X1\ Y1; \(Y = Y1; 
          \), 
\ \ \ Ctype == 4, X = X1; \(Y = 1\/Y1; \), 
\ \ \ 
          Ctype == 5, X = 1\/X1; \(Y = 1\/Y1; \), 
\ \ \ Ctype == 6, 
          X = Log[X1]; \(Y = Y1; \), 
\ \ \ Ctype == 7, X = X1; 
          \(Y = Log[Y1]; \), 
\ \ \ Ctype == 8, X = Log[X1]; Y = Log[Y1], 

          \ \ \ Ctype == 9, X = X1; \(Y = 1\/\@Y1; \), 
\ \ \ Ctype == 10, 
          X = X1; \(Y = Log[Y1\/X1]; \), 
\ \ \ Ctype == 11, X = X1; 
          Y = Log[L\/Y1 - 1]; \ ]; 
\ \ \ X = N[X]; 
\ \ \ Y = N[Y]; \ ]; 
    \)\)

\(Clear[Regression2, X0, Y0]; 

    Regression2[X0_, Y0_] := 
      Block[{SumX, SumXY, X, Y, Xmean, Ymean}, 
\ \ \ X = N[X0]; 
\ \ \ 
        Y = N[Y0]; 
\ \ \ 
        Xmean = ReplaceAll[X, List -> Plus]\/Length[X]; 
\ \ \ 
        Ymean = ReplaceAll[Y, List -> Plus]\/Length[Y]; 
\ \ \ 
        SumX = ReplaceAll[\((X - Xmean)\)\^2, List -> Plus]; 
\ \ \ 
        SumXY = ReplaceAll[\((X - Xmean)\)\ \((Y - Ymean)\), 
            List -> Plus]; 
\ \ \ a = SumXY\/SumX; 
\ \ \ 
        b = Ymean - a\ Xmean; ]; \)

Clear[Coeffs, Ctype]; \(\(Coeffs[Ctype_] := 
      Block[{}, 

        \(Which[
\ \ \ Ctype == 3, c = \(-\(1\/a\)\); d = \(-\(b\/a\)\), 

          \ \ \ Ctype == 7, c = E\^b, 
\ \ \ Ctype == 8, 
          c = E\^b, 
\ \ \ Ctype == 10, 
          c = E\^b; d = \(-a\), 
\ \ \ Ctype == 11, 
          c = E\^b]; \)\ ]; \)\)

Clear[CurveFit, Ctype]; \(CurveFit[Ctype_] := 
      Block[{}, 
\ \ \ 
        X = \(Transpose[XY]\)[[1]]; 

        \ \ \ Y = 
          \(Transpose[XY]\)[[2]]; 

        \ \ \ ChangeVariables[X, Y, Ctype]; 
\ \ \ Regression2[X, Y]; 

        \ \ \ Coeffs[Ctype]; 
\ \ \ SelectFunction[Ctype]; 
\ \ \ 
        Return[Together[f[x]]]; ]; 
\(On[General::"spell1"]; \)\)

\(Clear[RMS]; 

    RMS[Ctype_] := 
      Block[{k}, 
\ \ \ E2 = 0; 
\ \ \ 
        X = \(Transpose[XY]\)[[1]]; 

        \ \ \ Y = 
          \(Transpose[XY]\)[[2]]; 

        \ \ \ n\  = \ Length[X]; 
\ \ \ 
        For[k = 1, k ≤ n, \(k++\), 
\ \ \ \ \ \ 
          Z = a\ X[[k]] + b; 
          Err = Y[[k]] - Z; 
          E2 = N[E2 + Err\^2]]; 
\ \ \ E2 = \@\(E2\/n\); 
\ \ \ 
        Return[E2]; ]; 
\(On[General::"spell1"]; \)\)

Chapter 5, Examples for curve fitting, Pages 272-274.
Use the following data and fit several curves.

\(XY = {{1, 2}, {2, 3}, {3, 10}, {6, 25}, {12, 32}}; 
\(L = 40; \)\)

You can find a particular curve fit:

\(\(\ i = 6; 
CurveFit[i]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[f[x], {x, 1, 12}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["y = f[x] = ", f[x]]; \)\)

Or, you can look at all of the functions, and the RMS error  (Subscript[E, 2]), and try to decide 
which fit is best.  Remark.  The line seems best, but it was not intentional.

\(\(For[i = 1, i ≤ 11, \(i++\), 
\ \ \ CurveFit[i]; 
\ \ \ 
      Print[" "]; 
\ \ \ Print["f"\_i, "[x] = ", f[x]]; 
\ \ \ 
      Print[\*""The RMS error \!\(E\_2\) is: "", \ RMS[i]]; \ ]; \)\)

\(\(\ i = 8; 
CurveFit[i]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[f[x], {x, 1, 12}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["y = f[x] = ", f[x]]; \)\)

Or, you can look at all of the graphs.

\(\(For[i = 1, i ≤ 11, \(i++\), 
\ \ \ CurveFit[i]; 
\ \ \ 
      dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
          DisplayFunction -> Identity]; 
\ \ \ 
      gr = Plot[f[x], {x, 1, 12}, PlotStyle -> {RGBColor[1, 0, 0]}, 
          DisplayFunction -> Identity]; 
\ \ \ 
      Show[gr, dots, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y    "}]; 
\ \ \ 
      Print["f"\_i, "[x] = ", f[x]]; \ ]; \)\)

Algorithm 5.4 Cubic Splines.  To construct and evaluate the cubic spline 
interpolant S(x) for n+1 data points  ((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])), ⋯, 
      ((Subscript[x, n]), (Subscript[y, n])).
Provision is made for the following choices of endpoint constraints:
(i)   "Clamped Cubic Spline"; specify  S'((Subscript[x, 0]))  and  S'((Subscript[x, n])).
(ii)  "Natural Cubic Spline"  (A "relaxed curve.").
(iii)  Extrapolate  S''(x)  to the endpoints.
(iv)  S''(x)  is constant near the endpoints.
(v)   Specify  S''(x)  at each endpoint, i.e.   S''((Subscript[x, 0]))  and  S''((Subscript[x, n])).
Section  5.3   Interpolation by Spline Functions   Page  297

Execute this cell to activate Cubic Splines.
MakeSpline[Stype_]
Stype_ is the type of spline.
Remark. The points must be stored in  XY.

Off[General::"spell1"]; Clear[Differences, XY0, n]; \(\(Differences[XY0_] := 
	
      Block[{k}, 
		XY = XY0; 
		n = Length[XY] - 1; 
		
        X = \(Transpose[XY]\)[[1]]; 

        		Y = \(Transpose[XY]\)[[2]]; 
        
		h = Table[0, {n}]; 
		d = Table[0, {n}]; 
		
        m = Table[0, {n + 1}]; 
		a = Table[0, {n - 1}]; 
		
        b = Table[0, {n - 1}]; 
		c = Table[0, {n - 1}]; 
		
        v = Table[0, {n - 1}]; 
		s = Table[0, {n}, {4}]; 
		
        h[[1]] = 
          X[[2]] - 
            X[[1]]; 
		
        d[[1]] = 
          \(Y[[2]] - 
              Y[[1]]\)\/h
              [[1]]; 
		
        For[k = 2, k ≤ n, \(k++\), 
			
          h[[k]] = 
            X[[k + 1]] - 
              X[[k]]; 
			
          d[[k]] = 
            \(Y[[k + 1]] - 
                Y[[k]]\)\/h
                [[k]]; 
			
          a[[k - 1]] = 
            h[[k]]; 
			
          b[[k - 1]] = 
            2\ \((h[[k - 1]] + 
                  h[[k]])\); 
			
          c[[k - 1]] = 
            h[[k]]; 
			
          v[[k - 1]] = 
            6\ \((d[[k]] - 
                  d[[k - 1]])\)]; ]; 
    \)\)

Clear[Modify,Stype];
Modify[Stype_] :=
Module[{},
Which[Stype==1,
  txt=StringForm["Enter the first derivative S'(``)",X[[1]]];
  Dx1 =Input[txt];
  txt=StringForm["Enter the first derivative S'(``)",X[[n+1]]];
  DxN1=Input[txt];
  b[[1]]   = b[[1]] - h[[1]]/2;
  v[[1]]   = v[[1]] - 3 (d[[1]] - Dx1);
  b[[n-1]] = b[[n-1]] - h[[n]]/2;
  v[[n-1]] = v[[n-1]] - 3 (DxN1 - d[[n]]),
  Stype==2,
  m[[1]]   = 0;
  m[[n+1]] = 0,
  Stype==3,
  b[[1]] = b[[1]]+h[[1]]+h[[1]]h[[1]]/h[[2]];
  c[[1]] = c[[1]]-h[[1]]h[[1]]/h[[2]];
  b[[n-1]] = b[[n-1]]+h[[n]]+h[[n]]h[[n]]/h[[n-1]];
  a[[n-2]] = a[[n-2]]-h[[n]]h[[n]]/h[[n-1]],
  Stype==4,
  b[[1]] = b[[1]] + h[[1]];
  b[[n-1]] = b[[n-1]] + h[[n]],
  Stype==5,
  txt=StringForm["Enter the first derivative S''(``)",X[[1]]];
  DDx1 =Input[txt];
  txt=StringForm["Enter the first derivative S''(``)",X[[n+1]]];
  DDxN1=Input[txt];
  v[[1]] = v[[1]] - h[[1]] DDx1;
  v[[n-1]] = v[[n-1]] - h[[n]] DDxN1 ]; ];

Clear[TriDiagonal, a0, b0, c0, v0]; \(\(TriDiagonal[a0_, b0_, c0_, v0_] := 
	
      Module[{k, t}, 
		a = a0; 
		b = b0; 
		c = c0; 
		v = v0; 
        
		For[k = 2, k ≤ n - 1, \(k++\), 
			
          t = a[[k - 1]]\/b
                [[k - 1]]; 
			
          b[[k]] = 
            b[[k]] - 
              t\ c[[k - 1]]; 
			
          v[[k]] = 
            v[[k]] - 
              t\ v[[k - 1]]; ]; 
		
        m[[n]] = 
          v[[n - 1]]\/b
              [[n - 1]]; 
		
        For[k = n - 2, 1 ≤ k, \(k--\), 
			
          \(m[[k + 1]] = 
            \(v[[k]] - 
                c[[k]]\ 
                  m[[k + 2]]\)\/b
                [[k]]; \)]; ]; \)\)

Clear[Coeffs, Stype]; \(\(Coeffs[Stype_] := 
	
      Block[{}, 
		
        \(Which[
			Stype == 1, 
          m[[1]] = 
            \(3\ \((d[[1]] - Dx1)\)\)\/h
                  [[1]] - 
              m[[2]]\/2; 
          m[[n + 1]] = 
            \(3\ \((DxN1 - d[[n]])\)\)\/h
                  [[n]] - 
              m[[n]]\/2, 
			
          Stype == 2, m[[1]] = 0; 
          m[[n + 1]] = 0, 
			
          Stype == 3, 
          m[[1]] = 
            m[[2]] - 
              \(h[[1]]\ 
                  \((m[[3]] - 
                      m[[2]])\)\)\/h
                  [[2]]; 
          m[[n + 1]] = 
            m[[n]] + 
              \(h[[n]]\ 
                  \((m[[n]] - 
                      m[[n - 1]])\)\)\/h
                  [[n - 1]], 
			
          Stype == 4, 
          m[[1]] = 
            m[[2]]; 
          m[[n + 1]] = 
            m[[n]], 
			Stype == 5, 
          m[[1]] = DDx1; 
          m[[n + 1]] = DDxN1]; \)]; 
    \)\)

Clear[ComputeCoeff, d, h, m, y]; \(\(ComputeCoeff[d_, h_, m_, y_] := 
	
      Block[{k}, 
		
        \(For[k = 1, k ≤ n, \(k++\), 
			
          s[[k, 1]] = 
            Y[[k]]; 
			
          s[[k, 2]] = 
            d[[k]] - 
              1\/6\ h[[k]]\ 
                \((2\ m[[k]] + 
                    m[[k + 1]])\); 

          			s[[k, 3]] = 
            m[[k]]\/2; 
			
          s[[k, 4]] = 
            \(m[[k + 1]] - 
                m[[k]]\)\/\(6\ 
                h[[k]]\)]; \)]; \)\)

Clear[CS, t0]; \(\(CS[t0_] := 
	
      Block[{j, sum}, 
		t = N[t0]; 
		
        For[j = 1, j ≤ n, \(j++\), 
			
          If[X[[j]] ≤ t && 
              t < X[[j + 1]], k = j]]; 

        		If[t < X[[1]], k = 1]; 
		
        If[X[[n + 1]] ≤ t, 
          k = n]; 
		w = t - X[[k]]; 
        
		sum = 
          \((\((s[[k, 4]]\ w + 
                        s[[k, 3]])\)\ w
                     + s[[k, 2]])\)\ w + 
            s[[k, 1]]; 
		
        Return[sum]]; \)\)

Clear[MakeSpline, Stype]; \(MakeSpline[Stype_] := 
	
      Block[{}, 
		Differences[XY]; 
		Modify[Stype]; 
		
        TriDiagonal[a, b, c, v]; 
		Coeffs[Stype]; 
		
        ComputeCoeff[d, h, m, Y]]; 
\(On[General::"spell1"]; \)\)

Chapter 5, Examples 5.7-11, Pages 292-295.
Five different splines will be constructed for the following data.
You must enter the values
    S'(0) = 0.2
    S'(3) = -1
    S''(0) = -0.3
    S''(3) = 3.3
when prompted.

\(XY = {{0., 0.}, {1., 0.5}, {2., 2.}, {3., 1.5}}; 
MakeSpline[1]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 3}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 0.}, {1., 0.5}, {2., 2.}, {3., 1.5}}; 
MakeSpline[2]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 3}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 0.}, {1., 0.5}, {2., 2.}, {3., 1.5}}; 
MakeSpline[3]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 3}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 0.}, {1., 0.5}, {2., 2.}, {3., 1.5}}; 
MakeSpline[4]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 3}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 0.}, {1., 0.5}, {2., 2.}, {3., 1.5}}; 
MakeSpline[5]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 3}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

Chapter 5, Exercise 15, Page 303.
Five different splines will be constructed for the following data.
You must enter the values
    S'(0) = -0.6
    S'(6) = -1.8
    S''(0) = 1
    S''(6) = -1
when prompted.

Make sure you enter the number or Mathematica might bomb!

\(XY = {{0., 1.}, {1., 0.}, {2., 0.}, {3., 1.}, {4., 2.}, {5., 2.}, {6., 
          1.}}; 
MakeSpline[1]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 1.}, {1., 0.}, {2., 0.}, {3., 1.}, {4., 2.}, {5., 2.}, {6., 
          1.}}; 
MakeSpline[2]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 1.}, {1., 0.}, {2., 0.}, {3., 1.}, {4., 2.}, {5., 2.}, {6., 
          1.}}; 
MakeSpline[3]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 1.}, {1., 0.}, {2., 0.}, {3., 1.}, {4., 2.}, {5., 2.}, {6., 
          1.}}; 
MakeSpline[4]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

\(XY = {{0., 1.}, {1., 0.}, {2., 0.}, {3., 1.}, {4., 2.}, {5., 2.}, {6., 
          1.}}; 
MakeSpline[5]; 

    dots = ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    gr = Plot[CS[x], {x, 0, 6}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    Show[gr, dots, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"x", "y    "}]; 

    Print["Spline y = ", S\_i, "[x]"]; 

    Print["data points = ", XY]; \)

Algorithm 5.5  Trigonometric Polynomials.  To construct and evaluate the trigonometric polynomial of order  m  of the form
P (x) = (((Subscript[a, 0]))/2) + (∑_(j = 1)^m)( (Subscript[a, j]) cos (j x)  +  (Subscript[b, j]) sin (j x) ),
based on the  n  equally spaced values  (Subscript[x, k]) = -π + ((2 π k)/n).  
The construction is possible provided that  2m+1 ≤  n . 
Section  5.4    Fourier Series and Trig. Polynomials  Page 311

Execute this cell to activate Trigonometric Polynomials.
MakePoints[f0_,a0_,b0_,n_,mn_]
MakeCoeffs[XY0_,mn0_]
MakeTP[m0_]
TP[j,x_]
f0_ is the function f(x),
a0_ is the left endpoint,
b0_ is the right endpoint,
n_ is the number of points,
mn_ is the degree,
XY0_ are the points,
mn0_ is the degree,
m0_ is the degree,
j_ is the degree of the Chebyshev polynomial,
x_ is the variable.

Off[General::"spell1"]; Clear[MakePoints, f0, a0, b0, n, mn]; \(\(MakePoints[f0_, a0_, b0_, n_, mn_] := 
	
      Block[{f}, 
		Set@@{g[x_], g0}; 
		Clear[X, Y]; 
		
        X = Table[0, {n + 1}]; 
		Y = Table[0, {n + 1}]; 
		
        L = \(b0 - a0\)\/2; 
		m = mn; 
		m0 = m; 
		
        max = Floor[\(n - 1\)\/2]; 
		If[m0 > max, m0 = max]; 
		
        h = \(b0 - a0\)\/n; 
		
        For[k = 1, k ≤ n, \(k++\), 
          X[[k]] = 
            N[a0 + \((k - 1)\)\ h]]; 
		
        X[[n + 1]] = N[b0]; 
		
        Print["The maximum degree of a trigonometric polynomial is ", 
          m0]; 
		
        For[k = 1, k ≤ n, \(k++\), 
          Y[[k]] = 
            N[f[X[[k]]]]]; 
		
        Y[[n + 1]] = 
          N[f[X[[n + 1]]]]; 
		
        XY = Transpose[{X, Y}]; 
		Return[XY]]; \)\)

Clear[MakeCoeffs, XY0, mn0]; \(\(MakeCoeffs[XY0_, mn0_] := 
	
      Block[{X, Y}, 
		XY = XY0; 
		n = Length[XY] - 1; 
		
        X = \(Transpose[XY]\)[[1]]; 

        		Y = \(Transpose[XY]\)[[2]]; 
        
		Y[[1]] = 
          1\/2\ \((
              Y[[1]] + 
                Y[[n + 1]])\); 
		
        Y[[n + 1]] = 
          Y[[1]]; 
		
        L = 1\/2\ 
            \((X[[n + 1]] - 
                X[[1]])\); 
		
        m = mn0; 
		max = Floor[\(n - 1\)\/2]; 
		m0 = Min[m, max]; 

        		m = m0; 
		If[Mod[n, 2] == 1 && m0 < m, m = m0 + 1]; 
		
        a = Table[0, {max + 2}]; 
		b = Table[0, {max + 2}]; 
		
        For[k = 1, k ≤ n, \(k++\), 
          a[[1]] = 
            N[a[[1]] + 
                Y[[k]]]; 
			
          For[j = 2, j ≤ m + 1, \(j++\), 
				
            t = \(\((j - 1)\)\ π\ 
                  X[[k]]\)\/L; 
				
            a[[j]] = 
              N[a[[j]] + 
                  Y[[k]]\ Cos[t]]; 

            				b[[j]] = 
              N[b[[j]] + 
                  Y[[k]]\ Sin[t]]; ]; ]; 
        
		For[j = 1, j ≤ m + 1, \(j++\), 
			
          a[[j]] = 
            \(2\ a[[j]]\)\/n; 
			
          b[[j]] = 
            \(2\ b[[j]]\)\/n; ]; 
		
        Print["The coefficients for the trigonometric polynomial are:"]; 
        
		Print[Chop[a]]; 
		Print[Chop[b]]]; \)\)

Clear[MakeTP]; \(MakeTP[m0_] := 
	
      Block[{j}, 
		Clear[TP]; 
		c = Table[0, {max + 1}]; 
		
        sum = a[[1]]\/2; 
		
        TP[0, x_] = Chop[N[sum]]; 
		
        For[j = 1, j ≤ m0, \(j++\), 
			
          z = \(j\ π\ x\)\/L; 
			
          sum = sum + 
              a[[j + 1]]\ Cos[z] + 
              b[[j + 1]]\ Sin[z]; 

          			TP[j, x_] = Chop[N[sum]]]; 
		
        If[Mod[n, 2] == 1 && m0 < m, 
			z = \(m\ π\ x\)\/L; 
			
          sum = sum + 1\/2\ A[m + 1]\ Cos[z]; 
			
          TP[m, x_] = Chop[N[sum]]; ]; 
		
        Print["The trigonometric polynomial of degree ", m0, 
          " has been formed."]; ]; 
\(On[General::"spell1"]; \)\)

Chapter 5, Exercise 15(b), Page 313.
Find the trigonometric polynomial of degree n=5  for the 12  equally spaced points in the interval  [-π, π]  for the function:
          f (x ) = {({{(π/2) + x   for    - π ≤ x < 0, }, {(π/2) - x   for       0 ≤ x 
                  ≤ π.}})

Clear[f, x]f[x_ /; -π ≤ x < 0] = (π/2) + xf[x_ /; 0 ≤ x ≤ π] = (π/2) - x

MakePoints[(x/2), -π, π, 12, 5]; MakeCoeffs[XY, 5]; MakeTP[5]; TP[5, x]

dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; \(gr = 
      Plot[TP[5, x], {x, \(-π\), π}, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    \(Print["A trigonometric polynomial interpolation."]; \)\)Show[gr, dots, DisplayFunction -> $DisplayFunction];

Now let us do it with Mathematica's built in Fit procedure.
Remark.  Since the data is "periodic" we do not use both the left and right endpoints when finding the trigonometric polynomial via "fitting."  If both endpoints were used then the "fitting" process would minimize the error at both the left and right endpoints and thus weight these points twice, compared to those in the middle of the interval.

TQ[x_] = 
      Chop[Fit[Take[XY, 12], {1, Cos[x], Sin[x], Cos[2 x], Sin[2 x], 
            Cos[3 x], Sin[3 x], Cos[4 x], Sin[4 x], Cos[5 x], 
            Sin[5 x]}, x]]

dots = 
      ListPlot[XY, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; \(gr = 
      Plot[TQ[x], {x, \(-π\), π}, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]; 

    \(Print["A trigonometric polynomial interpolation."]; \)\)Show[gr, dots, DisplayFunction -> $DisplayFunction];

Chapter 5, Exercise 17, Page 314.
Fit the following data with a trigonometric polynomial including terms up to  cos(((2 π x)/13))  and   sin(((2 π x)/13)).
The yearly temperature cycle for Fairbanks, Alaska is:

xts = {{0, -14}, {1, -9}, {2, 2}, {3, 15}, {4, 35}, {5, 52}, {
          6, 62}, {7, 63}, {8, 58}, {9, 50}, {10, 34}, {11, 12}, {12, 
          -5}, {13, -14}}

Now let us do it with Mathematica's built in Fit procedure.

Clear[f, t]; pi = N[π]; f[t_] = 
      Fit[Take[xts, 13], {1, Cos[((2 pi t)/13)], Sin[((2 pi t)/13)]}, t]

\(graph = 
      Plot[f[t], {t, 0, 13}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        PlotRange -> {\(-15\), 65}, DisplayFunction -> Identity]; 

    \(dots = ListPlot[xts, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; \)\)Show[graph, dots, 
      PlotRange -> {{-0.5, 13.5}, {-20, 70}}, 
      DisplayFunction -> $DisplayFunction];

vals = Transpose[{X, Y, f[X]}]; TableForm[N[vals, 3]]

Remark. If a higher degree approximation is required then we must include the next two terms  cos(((4π x)/13))  and   sin(((4π x)/13)).

Clear[f, t]; pi = N[π]; f[t_] = 
      Fit[Take[xts, 13], {1, Cos[((2 pi t)/13)], Sin[((2 pi t)/13)], 
          Cos[((4 pi t)/13)], Sin[((4 pi t)/13)]}, t]

\(graph = 
      Plot[f[t], {t, 0, 13}, PlotStyle -> {RGBColor[1, 0, 0]}, 
        PlotRange -> {\(-15\), 65}, DisplayFunction -> Identity]; 

    \(dots = ListPlot[xts, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; \)\)Show[graph, dots, 
      PlotRange -> {{-0.5, 13.5}, {-20, 70}}, 
      DisplayFunction -> $DisplayFunction];

vals = Transpose[{X,Y,f[X]}];
TableForm[N[vals,3]]