(************** Content-type: application/mathematica **************
                     CreatedBy='Mathematica 4.2'

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     32100,        876]*)
(*NotebookOutlinePosition[     32799,        900]*)
(*  CellTagsIndexPosition[     32755,        896]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["Polynomial Interpolation Examples", "Title"],

Cell[CellGroupData[{

Cell["Interpolation of population data by a single polynomial", "Section"],

Cell["World population data (in billions) from 1950 through 2000", "Text",
  FontSize->16],

Cell[BoxData[
    \(popdata = {{1950, 2.556}, {1960, 3.039}, {1970, 3.706}, {1980, 
          4.453}, {1990, 5.278}, {2000, 6.082}}\)], "Input"],

Cell[BoxData[
    \(rawdata = 
      ListPlot[popdata, PlotStyle \[Rule] PointSize[0.02]]\)], "Input"],

Cell[CellGroupData[{

Cell["Brute-force interpolation by setting up a linear system", "Subsection"],

Cell["\<\
Setting up and solving a linear system to feed a 5th degree polynomial \
through these data points.\
\>", "Text",
  FontSize->14],

Cell[BoxData[
    \(\(m\  = \ 
        Table[\((1940 + i*10)\)^j, {i, 1, 6}, {j, 0, 5}];\)\)], "Input"],

Cell[BoxData[
    \(m\  // \ MatrixForm\)], "Input"],

Cell[BoxData[
    \(p = {2.556, 3.039, 3.706, 4.453, 5.278, 6.082}\)], "Input"],

Cell[BoxData[
    \(coeffs\  = \ LinearSolve[m, p]\)], "Input"],

Cell[BoxData[
    \(interpolatedpop[t_]\  = \ 
      coeffs . {1, t, t^2, t^3, t^4, t^5}\)], "Input"],

Cell[BoxData[
    \(polyfit = Plot[interpolatedpop[t], {t, 1945, 2005}]\)], "Input"],

Cell[BoxData[
    \(Show[polyfit, rawdata]\)], "Input"],

Cell["\<\
Using this interpolating polynomial to predict some past or future \
populations\
\>", "Text",
  FontSize->14],

Cell[BoxData[
    \(interpolatedpop[1975]\)], "Input"],

Cell[BoxData[
    \(interpolatedpop[2010]\)], "Input"],

Cell[BoxData[
    \(interpolatedpop[2050]\)], "Input"],

Cell["\<\
Yikes -- extrapolated a bit too far past data here.  (Note that polynomial \
form may not be best for a given type of data.  Here, we may have better \
reason to believe in an exponential fit, which we will learn how to do \
later.)\
\>", "Text",
  FontSize->14]
}, Closed]],

Cell[CellGroupData[{

Cell["Lagrangian interpolation", "Subsection"],

Cell["\<\
Here's another computation of the same interpolating polynomial, this time \
using the Lagrangian interpolating polynomials for the computations.\
\>", "Text",
  FontSize->14],

Cell[BoxData[{
    \(\(tdata\  = \ {1950, 1960, 1970, 1980, 1990, 
          2000};\)\), "\[IndentingNewLine]", 
    \(\(pdata = {2.556, 3.039, 3.706, 4.453, 5.278, 
          6.082};\)\), "\[IndentingNewLine]", 
    \(lagrangepoly[j_, t_, n_]\  := \ 
      Product[\((t - tdata[\([i]\)])\)/\((tdata[\([j]\)] - 
                tdata[\([i]\)])\), {i, 1, j - 1}]*
        Product[\((t - tdata[\([i]\)])\)/\((tdata[\([j]\)] - 
                tdata[\([i]\)])\), {i, j + 1, n}]\), "\[IndentingNewLine]", 
    \(\)}], "Input"],

Cell["\<\
Here's a plot of the 2nd Lagrange polynomial -- it should equal 1 at t=1960 \
and vanish at 1950, 1970, 1980, 1990, and 2000, and it does.\
\>", "Text",
  FontSize->14],

Cell[BoxData[
    \(Plot[lagrangepoly[2, t, 6], {t, 1945, 2005}]\)], "Input"],

Cell["\<\
Now we use the Lagrange polynomials to fit the population data and plot it \
(finding the same curve as we found the first time).\
\>", "Text",
  FontSize->14],

Cell[BoxData[
    \(lagrangefit[t_]\  = \ 
      Simplify[Sum[
          pdata[\([k]\)]*lagrangepoly[k, t, 6], {k, 1, 6}]]\)], "Input"],

Cell[BoxData[
    \(lagrangepolyplot = Plot[lagrangefit[t], {t, 1945, 2005}]\)], "Input"],

Cell[BoxData[
    \(Show[rawdata, lagrangepolyplot, polyfit]\)], "Input"]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{

Cell["Errors in polynomial interpolation", "Section"],

Cell["\<\
Here's a function that polynomial interpolation has difficulty with:\
\>", "Text",
  FontSize->16],

Cell[BoxData[{
    \(Clear[x]; f[x_] = 1/\((1 + 16*x^2)\);\), "\[IndentingNewLine]", 
    \(\(xdata\  = \ {\(-4\), \(-3\), \(-2\), \(-1\), 0, 1, 2, 3, 
          4};\)\), "\[IndentingNewLine]", 
    \(\(ydata = {f[\(-4\)], f[\(-3\)], f[\(-2\)], f[\(-1\)], f[0], f[1], 
          f[2], f[3], f[4]};\)\), "\[IndentingNewLine]", 
    \(\(lagrangepoly[j_, x_, n_]\  := \ 
        Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, 1, j - 1}]*
          Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, j + 1, 
              n}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangefit[x_]\  = \ 
        Simplify[
          Sum[ydata[\([k]\)]*lagrangepoly[k, x, 9], {k, 1, 
              9}]];\)\), "\[IndentingNewLine]", 
    \(\)}], "Input"],

Cell["\<\
The true function is in red, the polynomial interpolant (based on true values \
at x=0,1,2,3,4,-1,-2,-3,-4) is in black:\
\>", "Text",
  FontSize->16],

Cell[BoxData[
    \(Plot[{lagrangefit[x], f[x]}, {x, \(-4\), 4}, 
      PlotStyle \[Rule] {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}, 
      PlotRange \[Rule] {\(-5\), 5}]\)], "Input"],

Cell["\<\
What's going on?  Why such large error?  We know that the error theorem says \
the error depends on the 10th derivative of f, so let's plot it to see if it \
might be getting large...\
\>", "Text",
  FontSize->16],

Cell[BoxData[{
    \(\(tenthderiv[x_]\  = \ D[f[x], {x, 10}];\)\), "\[IndentingNewLine]", 
    \(Plot[tenthderiv[x], {x, \(-4\), 4}]\)}], "Input"],

Cell["\<\
Yep, pretty enormous.  Even divided by 10! (approximately 3 X 10^6), this \
suggests a potentially very large error.\
\>", "Text",
  FontSize->14],

Cell["\<\
One might hope we could do better by taking more points, say from x=-8 \
through x=8.  However,...\
\>", "Text",
  FontSize->16],

Cell[BoxData[{
    \(Clear[x]; f[x_] = 1/\((1 + 16*x^2)\);\), "\[IndentingNewLine]", 
    \(\(xdata\  = \ Table[j, {j, \(-8\), 8}];\)\), "\[IndentingNewLine]", 
    \(\(ydata = Table[f[j], {j, \(-8\), 8}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangepoly[j_, x_, n_]\  := \ 
        Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, 1, j - 1}]*
          Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, j + 1, 
              n}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangefit[x_]\  = \ 
        Simplify[
          Sum[ydata[\([k]\)]*lagrangepoly[k, x, 17], {k, 1, 
              17}]];\)\), "\[IndentingNewLine]", 
    \(\)}], "Input"],

Cell[BoxData[
    \(Plot[{lagrangefit[x], f[x]}, {x, \(-8\), 8}, 
      PlotStyle \[Rule] {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}, 
      PlotRange \[Rule] {\(-5\), 5}]\)], "Input"],

Cell["Even worse than before.", "Text",
  FontSize->16],

Cell["\<\
Surely, at least we will do better by taking a finer grid of points, say from \
-4 to 4 again, but in steps of 0.5...\
\>", "Text",
  FontSize->16],

Cell[BoxData[{
    \(Clear[x]; f[x_] = 1/\((1 + 16*x^2)\);\), "\[IndentingNewLine]", 
    \(\(xdata\  = \ 
        Table[j, {j, \(-4\), 4, 0.5}];\)\), "\[IndentingNewLine]", 
    \(\(ydata = 
        Table[f[j], {j, \(-4\), 4, 0.5}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangepoly[j_, x_, n_]\  := \ 
        Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, 1, j - 1}]*
          Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, j + 1, 
              n}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangefit[x_]\  = \ 
        Simplify[
          Sum[ydata[\([k]\)]*lagrangepoly[k, x, 17], {k, 1, 
              17}]];\)\), "\[IndentingNewLine]", 
    \(\)}], "Input"],

Cell[BoxData[
    \(Plot[{lagrangefit[x], f[x]}, {x, \(-4\), 4}, 
      PlotStyle \[Rule] {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}, 
      PlotRange \[Rule] {\(-5\), 5}]\)], "Input"],

Cell["\<\
Sorry, still worse than the original, at least in some ways: better near the \
bump, but much worse on the flat part.  High-degree polynomials just like to \
oscillate.\
\>", "Text",
  FontSize->14],

Cell["\<\
All hail the Chebyshev nodes!  Let's try this last example, only shift from \
the evenly-spaced nodes to the Chebyshev ones:\
\>", "Text",
  FontSize->14],

Cell[BoxData[{
    \(Clear[x]; f[x_] = 1/\((1 + 16*x^2)\);\), "\[IndentingNewLine]", 
    \(\(xdata\  = \ 
        Table[\(-4\)*Cos[j*Pi/16], {j, 0, 16}];\)\), "\[IndentingNewLine]", 
    \(\(ydata = 
        Table[f[\(-4\)*Cos[j*Pi/16]], {j, 0, 
            16}];\)\), "\[IndentingNewLine]", 
    \(chebyplot = 
      ListPlot[Transpose[{xdata, ydata}], PlotStyle \[Rule] PointSize[0.02], 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(\(lagrangepoly[j_, x_, n_]\  := \ 
        Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, 1, j - 1}]*
          Product[\((x - xdata[\([i]\)])\)/\((xdata[\([j]\)] - 
                  xdata[\([i]\)])\), {i, j + 1, 
              n}];\)\), "\[IndentingNewLine]", 
    \(\(lagrangefit[x_]\  := \ 
        Simplify[
          Sum[ydata[\([k]\)]*lagrangepoly[k, x, 17], {k, 1, 
              17}]];\)\), "\[IndentingNewLine]", 
    \(\)}], "Input"],

Cell[BoxData[{
    \(pl = 
      Plot[{lagrangefit[x], f[x]}, {x, \(-4\), 4}, 
        PlotStyle \[Rule] {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[chebyplot, pl, DisplayFunction \[Rule] $DisplayFunction, 
      PlotRange \[Rule] {\(-5\), 5}]\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["\<\
Richardson extrapolation: using interpolation to improve numerical \
approximations to derivatives\
\>", "Section"],

Cell["\<\
Application of polynomial interpolation (really, extrapolation): Richardson \
extrapolation to get a better estimate for f '(a).\
\>", "Text",
  FontSize->18],

Cell["\<\
Consider the function f(x) = 100/x at x = 5.  We know that f '(x) = - \
100/x^2, which means that f '(5) = -100/25 = -4 exactly.  Recall our \
\"centered difference\" approximation to f '(5):\
\>", "Text",
  FontSize->16],

Cell[BoxData[{
    \(Clear[x]; f[x_]\  = \ 100/x;\), "\[IndentingNewLine]", 
    \(\(h\  = \ 0.1`25;\)\), "\[IndentingNewLine]", 
    \(ctrddifferenceapproxderiv\  = \ \(f[5 + h] - f[5 - h]\)\/\(2*h\)\)}], \
"Input"],

Cell[TextData[{
  "We claimed that this method guarantees an O(",
  Cell[BoxData[
      \(TraditionalForm\`h\^2\)]],
  ") error, which should be around ",
  Cell[BoxData[
      \(TraditionalForm\`\((0.1)\)\^2\)]],
  "= 0.01 in this case.  (In reality we do one digit better)."
}], "Text",
  FontSize->16],

Cell[TextData[{
  "Now, we do something clever.  Since we know the error is O(",
  Cell[BoxData[
      \(TraditionalForm\`h\^2\)]],
  "), if we compute two centered-difference approximations at different \
values of h, we should be able to fit them to a parabola of the form A + B ",
  Cell[BoxData[
      \(TraditionalForm\`h\^2\)]],
  ", with A being the approximate value of f '(5).  Working out the details \
of this fit, we find the Richardson extrapolation formula implemented below, \
and which is predicted by theory to give an error of O(",
  Cell[BoxData[
      \(TraditionalForm\`h\^4\)]],
  ") in approximating f '(5).  Let's see:"
}], "Text",
  FontSize->16],

Cell[BoxData[{
    \(Clear[x]; f[x_]\  = \ 100/x;\), "\[IndentingNewLine]", 
    \(\(h\  = \ 0.1`25;\)\), "\[IndentingNewLine]", 
    \(\(ctr1\  = \ \(f[5 + h] - f[5 - h]\)\/\(2*h\);\)\), \
"\[IndentingNewLine]", 
    \(\(ctr2\  = \ \(f[5 + h/2] - f[5 - h/2]\)\/\(2*h/2\);\)\), "\
\[IndentingNewLine]", 
    \(richardsonderivapprox\  = \ 4*ctr2/3 - ctr1/3\)}], "Input"],

Cell["\<\
Wow, that's impressively better, at only twice the computational cost!\
\>", "Text",
  FontSize->16]
}, Closed]],

Cell[CellGroupData[{

Cell["Linear splines", "Section"],

Cell["\<\
French-franc-to-US-dollar exchange rates from 1971 to 2001 (thanks \
forecasts.org)\
\>", "Text",
  FontSize->14],

Cell[BoxData[{
    \(\(tdata\  = \ Table[i, {i, 1971, 2001}];\)\), "\[IndentingNewLine]", 
    \(\(exdata\  = \ {5.519, 5.174, 5.084, 5.028, 4.369, 4.476, 4.973, 4.718, 
          4.243, 4.041, 4.645, 5.830, 6.773, 8.595, 9.704, 7.482, 6.201, 
          5.581, 6.254, 5.737, 5.125, 5.386, 5.475, 5.921, 5.291, 5.012, 
          5.415, 6.083, 5.659, 6.475, 6.996};\)\), "\[IndentingNewLine]", 
    \(dataplot\  = \ 
      ListPlot[Transpose[{tdata, exdata}], 
        PlotStyle \[Rule] PointSize[0.02]]\)}], "Input"],

Cell[BoxData[
    \(\(\(piecewiselinear[
        t_]\)\(\ \)\(:=\)\(\[IndentingNewLine]\)\(Block[{}, \
\[IndentingNewLine]j = 1; \[IndentingNewLine]While[tdata[\([j]\)] < t\ , 
          j = j + 1]; \[IndentingNewLine]exdata[\([j - 
                1]\)] + \((exdata[\([j]\)] - 
                exdata[\([j - 1]\)])\)*\((t - 
                  tdata[\([j - 1]\)])\)/\((tdata[\([j]\)] - 
                  tdata[\([j - 
                        1]\)])\)\[IndentingNewLine]]\)\(\ \)\)\)], "Input"],

Cell[BoxData[{
    \(piecewiselinearplot = 
      Plot[piecewiselinear[t], {t, 1971, 2001}, PlotPoints \[Rule] 200, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[dataplot, piecewiselinearplot, 
      DisplayFunction \[Rule] $DisplayFunction]\[IndentingNewLine]\), "\
\[IndentingNewLine]", 
    \(\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Quadratic splines", "Section"],

Cell["\<\
French-franc-to-US-dollar exchange rates from 1971 to 2001 (thanks \
forecasts.org)\
\>", "Text",
  FontSize->14],

Cell[BoxData[{
    \(\(tdata\  = \ Table[i, {i, 1971, 2001}];\)\), "\[IndentingNewLine]", 
    \(\(exdata\  = \ {5.519, 5.174, 5.084, 5.028, 4.369, 4.476, 4.973, 4.718, 
          4.243, 4.041, 4.645, 5.830, 6.773, 8.595, 9.704, 7.482, 6.201, 
          5.581, 6.254, 5.737, 5.125, 5.386, 5.475, 5.921, 5.291, 5.012, 
          5.415, 6.083, 5.659, 6.475, 6.996};\)\), "\[IndentingNewLine]", 
    \(dataplot\  = \ 
      ListPlot[Transpose[{tdata, exdata}], 
        PlotStyle \[Rule] PointSize[0.015]]\)}], "Input"],

Cell[BoxData[{
    \(\(zvec = Table[0, {31}];\)\), "\[IndentingNewLine]", 
    \(\(zvec[\([1]\)] = \((exdata[\([2]\)] - 
              exdata[\([1]\)])\)/\((tdata[\([2]\)] - 
              tdata[\([1]\)])\);\)\), "\[IndentingNewLine]", 
    \(Do[\[IndentingNewLine]\(zvec[\([i]\)] = \(-zvec[\([i - 1]\)]\) + 
            2*\((exdata[\([i]\)] - exdata[\([i - 1]\)])\)/\((tdata[\([i]\)] - 
                    tdata[\([i - 1]\)])\);\)\[IndentingNewLine], {i, 2, 
        31}\[IndentingNewLine]]\)}], "Input"],

Cell[BoxData[
    \(\(\(degtwospline[
        t_]\)\(\ \)\(:=\)\(\[IndentingNewLine]\)\(Block[{}, \
\[IndentingNewLine]j = 1; \[IndentingNewLine]While[tdata[\([j]\)] < t\ , 
          j = j + 1]; \[IndentingNewLine]exdata[\([j - 1]\)] + 
          zvec[\([j - 1]\)]*\((t - tdata[\([j - 1]\)])\) + \((zvec[\([j]\)] - 
                zvec[\([j - 1]\)])\)*\((t - 
                    tdata[\([j - 1]\)])\)^2/\((2*\((tdata[\([j]\)] - 
                        tdata[\([j - 
                              1]\)])\)^2)\)\[IndentingNewLine]]\)\(\ \)\)\)], \
"Input"],

Cell[BoxData[{
    \(degtwosplineplot = 
      Plot[degtwospline[t], {t, 1971, 2001}, PlotPoints \[Rule] 1000, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[dataplot, degtwosplineplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"],

Cell[BoxData[
    \(degtwosplineplot = 
      Plot[degtwospline[t], {t, 1971, 1981}, 
        PlotPoints \[Rule] 1000]\)], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Cubic splines", "Section"],

Cell["\<\
French-franc-to-US-dollar exchange rates from 1971 to 2001 (thanks \
forecasts.org)\
\>", "Text",
  FontSize->14],

Cell[BoxData[{
    \(\(tdata\  = \ Table[i, {i, 1971, 2001}];\)\), "\[IndentingNewLine]", 
    \(\(exdata\  = \ {5.519, 5.174, 5.084, 5.028, 4.369, 4.476, 4.973, 4.718, 
          4.243, 4.041, 4.645, 5.830, 6.773, 8.595, 9.704, 7.482, 6.201, 
          5.581, 6.254, 5.737, 5.125, 5.386, 5.475, 5.921, 5.291, 5.012, 
          5.415, 6.083, 5.659, 6.475, 6.996};\)\), "\[IndentingNewLine]", 
    \(dataplot\  = \ 
      ListPlot[Transpose[{tdata, exdata}], 
        PlotStyle \[Rule] PointSize[0.015]]\)}], "Input"],

Cell[BoxData[{
    \(\(hvec = 
        Table[tdata[\([i + 1]\)] - tdata[\([i]\)], {i, 1, 
            30}];\)\), "\[IndentingNewLine]", 
    \(\(uvec = 
        Table[\((hvec[\([i]\)] + hvec[\([i + 1]\)])\)*2, {i, 1, 
            29}];\)\), "\[IndentingNewLine]", 
    \(\(vvec = 
        Table[6*\((exdata[\([i + 2]\)] - exdata[\([i + 1]\)])\)/
                hvec[\([i + 1]\)] - 
            6*\((exdata[\([i + 1]\)] - exdata[\([i]\)])\)/hvec[\([i]\)], {i, 
            1, 29}];\)\), "\[IndentingNewLine]", 
    \(\(mat\  = \ Table[0, {31}, {31}];\)\), "\[IndentingNewLine]", 
    \(\(mat[\([1, 1]\)] = 1;\)\), "\[IndentingNewLine]", 
    \(\(mat[\([31, 31]\)] = 1;\)\), "\[IndentingNewLine]", 
    \(\(rhs = Table[0, {31}];\)\), "\[IndentingNewLine]", 
    \(\(Do[\[IndentingNewLine]mat[\([j, j - 1]\)] = 
          hvec[\([j - 1]\)]; \[IndentingNewLine]mat[\([j, j]\)] = 
          uvec[\([j - 1]\)]; \[IndentingNewLine]mat[\([j, j + 1]\)] = 
          hvec[\([j]\)]; \[IndentingNewLine]rhs[\([j]\)] = 
          vvec[\([j - 1]\)];\[IndentingNewLine], {j, 2, 
          30}\[IndentingNewLine]];\)\), "\[IndentingNewLine]", 
    \(\(zvec = LinearSolve[mat, rhs];\)\)}], "Input"],

Cell[BoxData[
    \(\(\(cubicspline[
        t_]\)\(\ \)\(:=\)\(\[IndentingNewLine]\)\(Block[{}, \
\[IndentingNewLine]j = 1; \[IndentingNewLine]While[tdata[\([j]\)] < t\ , 
          j = j + 1]; \[IndentingNewLine]h = 
          hvec[\([j - 
                1]\)]; \[IndentingNewLine]zvec[\([j]\)]*\((t - 
                    tdata[\([j - 1]\)])\)^3/\((6*h)\) + 
          zvec[\([j - 1]\)]*\((tdata[\([j]\)] - t)\)^3/\((6*
                  h)\) + \((exdata[\([j]\)]/h - h*zvec[\([j]\)]/6)\)*\((t - 
                tdata[\([j - 1]\)])\) + \((exdata[\([j - 1]\)]/h - 
                h*zvec[\([j - 1]\)]/6)\)*\((tdata[\([j]\)] - 
                t)\)\[IndentingNewLine]]\)\(\ \)\)\)], "Input"],

Cell[BoxData[{
    \(cubicsplineplot = 
      Plot[cubicspline[t], {t, 1971, 2001}, PlotPoints \[Rule] 1000, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[dataplot, cubicsplineplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"],

Cell[BoxData[
    \(cubicsplineplot = 
      Plot[cubicspline[t], {t, 1971, 1981}, 
        PlotPoints \[Rule] 1000]\)], "Input"],

Cell["Nicer than the quadratic spline, yes?", "Text",
  FontSize->18]
}, Closed]],

Cell[CellGroupData[{

Cell["Least squares", "Section"],

Cell[CellGroupData[{

Cell["Linear regression (least squares fit to a linear function)", \
"Subsection"],

Cell[BoxData[{
    \(\(tdata\  = \ Table[i, {i, 1971, 2001}];\)\), "\[IndentingNewLine]", 
    \(\(exdata\  = \ {5.519, 5.174, 5.084, 5.028, 4.369, 4.476, 4.973, 4.718, 
          4.243, 4.041, 4.645, 5.830, 6.773, 8.595, 9.704, 7.482, 6.201, 
          5.581, 6.254, 5.737, 5.125, 5.386, 5.475, 5.921, 5.291, 5.012, 
          5.415, 6.083, 5.659, 6.475, 6.996};\)\), "\[IndentingNewLine]", 
    \(dataplot\  = \ 
      ListPlot[Transpose[{tdata, exdata}], 
        PlotStyle \[Rule] PointSize[0.015]]\)}], "Input"],

Cell[BoxData[{
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%31 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%31 
                    tdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%31 
                    tdata[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%31 1\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"rhs", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%31 tdata[\([j]\)]*exdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%31 exdata[\([j]\)]\)}
              }], ")"}]}], 
      ";"}], "\[IndentingNewLine]", \(vec\  = \ 
        LinearSolve[a, rhs];\), "\[IndentingNewLine]", \(leastsquares[
        t_]\  = \ vec[\([1]\)]*t + vec[\([2]\)]\)}], "Input"],

Cell[BoxData[{
    \(lsplot = 
      Plot[leastsquares[t], {t, 1971, 2001}, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[dataplot, lsplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"],

Cell["\<\
Back to the population data, now fit by a least-squares line.\
\>", "Text",
  FontSize->18],

Cell[BoxData[{
    \(\(tdata = {1950, 1960, 1970, 1980, 1990, 
          2000};\)\), "\[IndentingNewLine]", 
    \(\(popdata = {2.556, 3.039, 3.706, 4.453, 5.278, 6.082};\)\)}], "Input"],

Cell[BoxData[
    \(rawdata = 
      ListPlot[Transpose[{tdata, popdata}], 
        PlotStyle \[Rule] PointSize[0.02]]\)], "Input"],

Cell[BoxData[{
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 1\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"rhs", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 tdata[\([j]\)]*popdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 popdata[\([j]\)]\)}
              }], ")"}]}], 
      ";"}], "\[IndentingNewLine]", \(vec\  = \ 
        LinearSolve[a, rhs];\), "\[IndentingNewLine]", \(leastsquares[
        t_]\  = \ vec[\([1]\)]*t + vec[\([2]\)]\)}], "Input"],

Cell[BoxData[{
    \(lsplot = 
      Plot[leastsquares[t], {t, 1950, 2000}, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[rawdata, lsplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Least-squares quadratic fit", "Subsection"],

Cell["Maybe a least-square parabola would be smarter:", "Text",
  FontSize->18],

Cell[BoxData[{
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^4\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 1\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"rhs", "=", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 tdata[\([j]\)]^2*popdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 tdata[\([j]\)]*popdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 popdata[\([j]\)]\)}
              }], ")"}]}], 
      ";"}], "\[IndentingNewLine]", \(vec\  = \ 
        LinearSolve[a, rhs];\), "\[IndentingNewLine]", \(leastsquares[
        t_]\  = \ 
      vec[\([1]\)]*t^2 + vec[\([2]\)]*t + vec[\([3]\)]\)}], "Input"],

Cell[BoxData[{
    \(lsplot = 
      Plot[leastsquares[t], {t, 1950, 2000}, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[rawdata, lsplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Numerical Considerations", "Subsection"],

Cell["Here's the original matrix and its condition number:", "Text",
  FontSize->18],

Cell[BoxData[{\(tdata = {1950. , 1960. , 1970. , 1980. , 1990. , 
          2000. };\), "\n", \(popdata = {2.556, 3.039, 3.706, 4.453, 5.278, 
          6.082};\), "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^4\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 1\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", \(<< 
      LinearAlgebra`MatrixManipulation`\), "\[IndentingNewLine]", \
\(MatrixConditionNumber[a]\)}], "Input"],

Cell["\<\
Now the result of rescaling time, so 1950 is t=0 and 2000 is t=1:\
\>", "Text",
  FontSize->18],

Cell[BoxData[{\(tdata = {0. , 0.2, 0.4, 0.6, 0.8, 
          1.0};\), "\n", \(popdata = {2.556, 3.039, 3.706, 4.453, 5.278, 
          6.082};\), "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^4\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^3\), \(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 
                      tdata[\([j]\)]^2\), \(\[Sum]\+\(j = 1\)\%6 
                    tdata[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 1\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", \(<< 
      LinearAlgebra`MatrixManipulation`\), "\[IndentingNewLine]", \
\(MatrixConditionNumber[a]\)}], "Input"],

Cell["\<\
Now the matrix for the Chebyshev basis and its condition number:\
\>", "Text",
  FontSize->18],

Cell[BoxData[{\(tdata = {0. , 0.2, 0.4, 0.6, 0.8, 
          1.0};\), "\n", \(popdata = {2.556, 3.039, 3.706, 4.453, 5.278, 
          6.082};\), "\[IndentingNewLine]", \(g0 = {1. , 1. , 1. , 1. , 1. , 
          1. };\), "\[IndentingNewLine]", \(g1 = 
        tdata;\), "\[IndentingNewLine]", \(g2 = 
        2*tdata^2 - 1;\), "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g0[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g0[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g0[\([j]\)]\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", \(<< 
      LinearAlgebra`MatrixManipulation`\), "\[IndentingNewLine]", \
\(MatrixConditionNumber[a]\)}], "Input"],

Cell["\<\
(A mild effect this time, but can be dramatic for larger data sets, \
higher-degree polynomials)\
\>", "Text",
  FontSize->18],

Cell["\<\
Now let's compute the solution using the Chebyshev basis, just to verify it's \
the same answer as before:\
\>", "Text",
  FontSize->18],

Cell[BoxData[{
    RowBox[{
      RowBox[{"a", " ", "=", " ", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*
                    g0[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*
                    g0[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g2[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g1[\([j]\)]\), \(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*
                    g0[\([j]\)]\)}
              }], ")"}]}], ";"}], "\[IndentingNewLine]", 
    RowBox[{
      RowBox[{"rhs", "=", 
        RowBox[{"(", GridBox[{
              {\(\[Sum]\+\(j = 1\)\%6 g2[\([j]\)]*popdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g1[\([j]\)]*popdata[\([j]\)]\)},
              {\(\[Sum]\+\(j = 1\)\%6 g0[\([j]\)]*popdata[\([j]\)]\)}
              }], ")"}]}], 
      ";"}], "\[IndentingNewLine]", \(vec\  = \ 
        LinearSolve[a, rhs];\), "\[IndentingNewLine]", \(leastsquares[
        t_]\  = \ 
      vec[\([1]\)]*\((2*\((\((t - 1950)\)/50)\)^2 - 1)\) + 
        vec[\([2]\)]*\((t - 1950)\)/50 + vec[\([3]\)]\)}], "Input"],

Cell[BoxData[{
    \(lsplot = 
      Plot[leastsquares[t], {t, 1950, 2000}, 
        DisplayFunction \[Rule] Identity]\), "\[IndentingNewLine]", 
    \(Show[rawdata, lsplot, 
      DisplayFunction \[Rule] $DisplayFunction]\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Nonlinear least squares", "Subsection"],

Cell[TextData[{
  "Fitting the population data to an equation of the form y = b ",
  Cell[BoxData[
      \(TraditionalForm\`\[ExponentialE]\^\(c(t - 1950)\)\)]]
}], "Text",
  FontSize->18],

Cell[BoxData[{
    \(\(tdata = {0, 10, 20, 30, 40, 50};\)\), "\n", 
    \(\(popdata = {2.556, 3.039, 3.706, 4.453, 5.278, 
          6.082};\)\), "\[IndentingNewLine]", 
    \(\(Clear[a, b, c];\)\), "\[IndentingNewLine]", 
    \(\(eqn1 = \ \[Sum]\+\(j = 1\)\%6\((\((popdata[\([j]\)] - 
                    b*\[ExponentialE]\^\(c*tdata[\([j]\)]\))\)*\[ExponentialE]\
\^\(c*tdata[\([j]\)]\))\) \[Equal] 0;\)\), "\[IndentingNewLine]", 
    \(\(eqn2 = \ \[Sum]\+\(j = 1\)\%6\((\((popdata[\([j]\)] - 
                    b*\[ExponentialE]\^\(c*tdata[\([j]\)]\))\)*b*
                tdata[\([j]\)]*\[ExponentialE]\^\(c*tdata[\([j]\)]\))\) \
\[Equal] 0;\)\), "\[IndentingNewLine]", 
    \(FindRoot[{eqn1, eqn2}, {b, 2.3}, {c, 0.01}]\)}], "Input"],

Cell[BoxData[{
    \(\(leastsquares[t_]\  = \ 
        2.60785*Exp[0.0172297*\((t - 1950)\)];\)\), "\[IndentingNewLine]", 
    \(\(lsplot = 
        Plot[leastsquares[t], {t, 1950, 2000}, 
          DisplayFunction \[Rule] Identity];\)\), "\n", 
    \(\(Show[rawdata, lsplot, 
        DisplayFunction \[Rule] $DisplayFunction];\)\), \
"\[IndentingNewLine]", 
    \(\)}], "Input"]
}, Closed]]
}, Closed]]
},
FrontEndVersion->"4.2 for Microsoft Windows",
ScreenRectangle->{{0, 1024}, {0, 685}},
WindowSize->{901, 638},
WindowMargins->{{0, Automatic}, {Automatic, 0}},
Magnification->1.5,
StyleDefinitions -> "Classroom.nb"
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 50, 0, 99, "Title"],

Cell[CellGroupData[{
Cell[1829, 55, 74, 0, 93, "Section"],
Cell[1906, 57, 90, 1, 48, "Text"],
Cell[1999, 60, 144, 2, 98, "Input"],
Cell[2146, 64, 102, 2, 72, "Input"],

Cell[CellGroupData[{
Cell[2273, 70, 77, 0, 78, "Subsection"],
Cell[2353, 72, 139, 4, 73, "Text"],
Cell[2495, 78, 103, 2, 72, "Input"],
Cell[2601, 82, 52, 1, 72, "Input"],
Cell[2656, 85, 79, 1, 72, "Input"],
Cell[2738, 88, 63, 1, 72, "Input"],
Cell[2804, 91, 101, 2, 72, "Input"],
Cell[2908, 95, 84, 1, 72, "Input"],
Cell[2995, 98, 55, 1, 72, "Input"],
Cell[3053, 101, 120, 4, 44, "Text"],
Cell[3176, 107, 54, 1, 72, "Input"],
Cell[3233, 110, 54, 1, 72, "Input"],
Cell[3290, 113, 54, 1, 72, "Input"],
Cell[3347, 116, 272, 6, 102, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[3656, 127, 46, 0, 52, "Subsection"],
Cell[3705, 129, 185, 4, 73, "Text"],
Cell[3893, 135, 522, 10, 202, "Input"],
Cell[4418, 147, 178, 4, 73, "Text"],
Cell[4599, 153, 77, 1, 72, "Input"],
Cell[4679, 156, 169, 4, 73, "Text"],
Cell[4851, 162, 135, 3, 98, "Input"],
Cell[4989, 167, 89, 1, 72, "Input"],
Cell[5081, 170, 73, 1, 72, "Input"]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{
Cell[5203, 177, 53, 0, 61, "Section"],
Cell[5259, 179, 108, 3, 48, "Text"],
Cell[5370, 184, 822, 16, 332, "Input"],
Cell[6195, 202, 160, 4, 81, "Text"],
Cell[6358, 208, 180, 3, 124, "Input"],
Cell[6541, 213, 223, 5, 114, "Text"],
Cell[6767, 220, 146, 2, 98, "Input"],
Cell[6916, 224, 156, 4, 73, "Text"],
Cell[7075, 230, 138, 4, 81, "Text"],
Cell[7216, 236, 732, 14, 306, "Input"],
Cell[7951, 252, 180, 3, 124, "Input"],
Cell[8134, 257, 55, 1, 48, "Text"],
Cell[8192, 260, 157, 4, 81, "Text"],
Cell[8352, 266, 760, 16, 306, "Input"],
Cell[9115, 284, 180, 3, 124, "Input"],
Cell[9298, 289, 208, 5, 73, "Text"],
Cell[9509, 296, 164, 4, 73, "Text"],
Cell[9676, 302, 958, 20, 358, "Input"],
Cell[10637, 324, 334, 6, 176, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[11008, 335, 125, 3, 102, "Section"],
Cell[11136, 340, 168, 4, 89, "Text"],
Cell[11307, 346, 231, 5, 114, "Text"],
Cell[11541, 353, 216, 4, 150, "Input"],
Cell[11760, 359, 304, 9, 81, "Text"],
Cell[12067, 370, 671, 15, 213, "Text"],
Cell[12741, 387, 369, 7, 236, "Input"],
Cell[13113, 396, 110, 3, 48, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[13260, 404, 33, 0, 61, "Section"],
Cell[13296, 406, 123, 4, 44, "Text"],
Cell[13422, 412, 515, 8, 228, "Input"],
Cell[13940, 422, 494, 9, 228, "Input"],
Cell[14437, 433, 345, 7, 228, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[14819, 445, 36, 0, 61, "Section"],
Cell[14858, 447, 123, 4, 44, "Text"],
Cell[14984, 453, 516, 8, 228, "Input"],
Cell[15503, 463, 506, 8, 254, "Input"],
Cell[16012, 473, 558, 10, 254, "Input"],
Cell[16573, 485, 282, 5, 150, "Input"],
Cell[16858, 492, 131, 3, 98, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[17026, 500, 32, 0, 61, "Section"],
Cell[17061, 502, 123, 4, 44, "Text"],
Cell[17187, 508, 516, 8, 228, "Input"],
Cell[17706, 518, 1182, 22, 462, "Input"],
Cell[18891, 542, 694, 12, 306, "Input"],
Cell[19588, 556, 279, 5, 150, "Input"],
Cell[19870, 563, 129, 3, 98, "Input"],
Cell[20002, 568, 69, 1, 52, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[20108, 574, 32, 0, 61, "Section"],

Cell[CellGroupData[{
Cell[20165, 578, 82, 1, 78, "Subsection"],
Cell[20250, 581, 516, 8, 228, "Input"],
Cell[20769, 591, 827, 18, 234, "Input"],
Cell[21599, 611, 237, 5, 124, "Input"],
Cell[21839, 618, 101, 3, 52, "Text"],
Cell[21943, 623, 186, 3, 98, "Input"],
Cell[22132, 628, 131, 3, 98, "Input"],
Cell[22266, 633, 823, 18, 234, "Input"],
Cell[23092, 653, 236, 5, 124, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[23365, 663, 49, 0, 52, "Subsection"],
Cell[23417, 665, 79, 1, 52, "Text"],
Cell[23499, 668, 1261, 26, 310, "Input"],
Cell[24763, 696, 236, 5, 124, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[25036, 706, 46, 0, 52, "Subsection"],
Cell[25085, 708, 84, 1, 52, "Text"],
Cell[25172, 711, 1033, 19, 261, "Input"],
Cell[26208, 732, 105, 3, 52, "Text"],
Cell[26316, 737, 1015, 19, 261, "Input"],
Cell[27334, 758, 104, 3, 52, "Text"],
Cell[27441, 763, 1288, 23, 365, "Input"],
Cell[28732, 788, 136, 4, 89, "Text"],
Cell[28871, 794, 146, 4, 89, "Text"],
Cell[29020, 800, 1415, 28, 397, "Input"],
Cell[30438, 830, 236, 5, 124, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[30711, 840, 45, 0, 52, "Subsection"],
Cell[30759, 842, 188, 5, 52, "Text"],
Cell[30950, 849, 740, 12, 322, "Input"],
Cell[31693, 863, 379, 9, 176, "Input"]
}, Closed]]
}, Closed]]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)