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 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

Algorithm 5.1 (Least Squares Line). To construct the least squares
line  y = Ax + B  that fits the  n  data points  (x1,y1), (x2,y2), ..., (xn,yn).
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

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 = E2 + Err^2//N ];
  E2 = Sqrt[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:

xy = {{-1,10},{0,9},{1,7},{2,5},{3,4},{4,3},{5,0},{6,-1}}

Clear[f,x];
f[x_] = Regression[xy]//N

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
gr = Plot[f[x],{x,-1,6},DisplayFunction->Identity];
Print["The least squares line: y = ",f[x]];
Show[gr,dots,AxesLabel->{"x","y"},
     DisplayFunction->$DisplayFunction];

The RMS error E2 for this example is:

E2

f[x]

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

Fit[xy,{1,x},x]

Algorithm 5.2 (Least Squares Polynomial). To construct the 
least squares polynomial of degree  m  of the form:
		Pm(x) = c1 + c2x + c3x2 + c4x3 +... + cmxm-1 + cm+1xm 
that fits the  n  data points  (x1,y1), (x2,y2),..., (xn,yn). 
Section	5.2	Curve Fitting   Page  278

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,{2m+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:

xy = {{-3,3},{0,1},{2,1},{4,3}}

P[x_] = LeastSqPoly[xy,2]

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
gr = Plot[P[x],{x,-3,4},DisplayFunction->Identity];
Print["The least squares polynomial: y = ",P[x]];
Show[gr,dots,AxesLabel->{"x","y"},
     DisplayFunction->$DisplayFunction];

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

Fit[xy,{1,x,x^2},x]

P[x]//N

Algorithm 5.3 (Non-linear Curve Fitting) To construct a non-linear
fit  y = g(x)  given the  n  data points  (x1,y1), (x2,y2),..., (xn,yn)  by using "data linearization."   Available curves are:

y = A/x + B,  y = D/(x + C),  y = 1/(Ax + B),  y = x/(A + Bx),  y = A ln(x) + B,

y = C eAx,	y = C xA,	y = (Ax + B)-2,  y = C x e-Dx,	y = L/(1 + C eAx).
	
Section	5.2	Curve Fitting   Page  280

First, execute the following group of cells:

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 Exp[a x],
      Ctype==8, f[x_] = c x^a,
      Ctype==9, f[x_] = 1/(a x + b)^2,
      Ctype==10,f[x_] = c x Exp[- d x],
      Ctype==11,f[x_] = L/(1 + c Exp[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=Map[Log,X1]; Y=Y1;,
      Ctype==7, X=X1;          Y=Map[Log,Y1];,
      Ctype==8, X=Map[Log,X1]; Y=Map[Log,Y1],
      Ctype==9, X=X1     ;     Y=Map[Sqrt,1/Y1];,
      Ctype==10,X=X1;          Y=Map[Log,Y1/X1];,
      Ctype==11,X=X1;          Y=Map[Log,L/Y1-1]; ];
      X = X//N;
      Y = Y//N; ];

Clear[Regression2,X0,Y0];
Regression2[X0_,Y0_] := 
Block[{SumX,SumXY,X,Y,Xmean,Ymean},
   X = X0//N;
   Y = Y0//N;
   Xmean = (X /.List->Plus)/Length[X];
   Ymean = (Y /.List->Plus)/Length[Y];
   SumX  = (X-Xmean)^2 /.List->Plus;
   SumXY = (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 = Exp[b],
      Ctype==8,  c = Exp[b],
      Ctype==10, c = Exp[b]; d = -a,
      Ctype==11, c = Exp[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[f[x]//Together]; ];
 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[i,x],{x,1,12},DisplayFunction->Identity];
txt = StringForm["y = ``",F[i,x]];
Show[gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y    "}];

Or, you can look at them all:

For[i=1,i<=11,i++,
    F[i,x_] = CurveFit[i];
    Print["F[",i,",x] = ",F[i,x]];
    Print[" "] ];

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
For[i=1,i<=11,i++,
    gr = Plot[F[i,x],{x,1,12},DisplayFunction->Identity];
    txt = StringForm["y = ``",F[i,x]];
    Show[gr,dots,DisplayFunction->$DisplayFunction,
    PlotLabel->txt,AxesLabel->{"x","y    "}];];

Algorithm 5.4 (Cubic Splines) To construct and evaluate the cubic
spline interpolant S(x) for n+1 data points  (x0,y0), (x1,y1),..., (xn,yn).  
Provision is made for the following choices of endpoint constraints:
	(i)	 "Clamped Cubic Spline"; specify  S'(x0)  and  S'(xn).
	(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''(x0)  and  S''(xn).
Section	5.3	Interpolation by Spline Functions   Page  297

First, execute the following group of cells:

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]]-h[[k]](2 m[[k]]+m[[k+1]])/6;
    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 = t0//N;
  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.

Make sure you enter the number or Mathematica might bomb!

XY = {{0.,0.},
      {1.,0.5},
      {2.,2.},
      {3.,1.5}}

You can look at a particular spline:

i=2;
MakeSpline[i];
dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
gr = Plot[CS[x],{x,0,3},DisplayFunction->Identity];
txt = StringForm["Spline  y = S[``,x]",i];
Show[gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y    "}];

Or, you can look at all of the splines:

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
For[i=1,i<=5,i++,
  MakeSpline[i];
  gr = Plot[CS[x],{x,0,3},DisplayFunction->Identity];
  txt = StringForm["Spline  y = S[``,x]",i];
  Show[gr,dots,DisplayFunction->$DisplayFunction,
  PlotLabel->txt,AxesLabel->{"x","y    "}];];

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.}}

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
                DisplayFunction->Identity];
For[i=1,i<=5,i++,
  MakeSpline[i];
  gr = Plot[CS[x],{x,0,6},DisplayFunction->Identity];
  txt = StringForm["Spline  y = S[``,x]",i];
  Show[gr,dots,DisplayFunction->$DisplayFunction,
  PlotLabel->txt,AxesLabel->{"x","y    "}];];

Algorithm 5.5  (Trigonometric Polynomials).  To construct  and evaluate the trigonometric polynomial of order  m  of the form
		P(x)  =  a0/2  +  Sj=1j=mÊ(ajcos(jÊx)Ê+Êbjsin(jÊx)Ê), 
based on the  n  equally spaced values  xk = - p + 2pk/n .  
The construction is possible provided that  2m+1 <= n. 
Section 5.4 Fourier Series and Trig. Polynomials  Page 311

First, execute the following group of cells:

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]] = a0 + (k-1) h//N ];
  X[[n+1]] = b0//N;
  Print["The maximum degree of a trigonometric polynomial is ",m0];
  For[k=1,k<=n,k++,
    Y[[k]] = f[X[[k]]]//N ];
  Y[[n+1]] = f[X[[n+1]]]//N; 
  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]] = (Y[[1]]+Y[[n+1]])/2;
  Y[[n+1]] = Y[[1]];
  L =(X[[n+1]]-X[[1]])/2;
  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]] = a[[1]] + Y[[k]]//N;
    For[j=2,j<=m+1,j++,
      t = (j-1) Pi X[[k]]/L;
      a[[j]]=a[[j]]+Y[[k]]Cos[t]//N;
      b[[j]]=b[[j]]+Y[[k]]Sin[t]//N; ]; ];
  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[sum//N];
  For[j=1,j<=m0,j++,
    z = j Pi x/L;
    sum = sum + a[[j+1]]Cos[z] + b[[j+1]]Sin[z];
    TP[j,x_] = Chop[sum//N] ];
  If[Mod[n,2]==1 && m0<m, 
    z = m Pi x/L;
    sum = sum + A[m+1] Cos[z]/2;
    TP[m,x_] = Chop[sum//N]; ]; 
  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 [-p,p] for the function:

f[x_/;N[-Pi]<=x<0] = Pi/2 +x
f[x_/;0<=x<=N[Pi]] = Pi/2 -x

MakePoints[x/2,-Pi,Pi,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,-Pi,Pi},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_] = Fit[Take[XY,12],{1,Cos[x],Sin[x],
         Cos[2x],Sin[2x],Cos[3x],Sin[3x],
         Cos[4x],Sin[4x],Cos[5x],Sin[5x]},x]//Chop

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
gr = Plot[TQ[x],{x,-Pi,Pi},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(2px/13) and cos(2px/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=Pi//N;
f[t_] = Fit[Take[xts,13],
        {1,Cos[2pi t/13],Sin[2pi t/13]},t]

graph = Plot[f[t],{t,0,13},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(4px/13) and cos(4px/13).

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

graph = Plot[f[t],{t,0,13},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]]