Excercise sheet 3

Interpolation

Tables

1. Make a table called DataA, that has as its elements ((Subscript[x, i]),  (Subscript[y, i])) the following table
	(Subscript[x, i])	(Subscript[y, i])
	1	4.0
	3	2.8
	4	2.5
	5	3.2
	6	3.9
	10	5.0
2. Make a table called DataB that has as its elements ((Subscript[x, i]),  ((Subscript[x, i])^2) + 10 (Subscript[x, i]) + 1), where (Subscript[x, i]) = 0, 1, 2,  ... , 5
3. Read a table from file "data.dat" and assign it to DataC.
4. Plot the tables DataA, DataB, and DataC
5. Add to the start of the table DataA the element (0, 6.0).
6. Add to the end of DataB the element (6,  0.0).
7. Insert (2,  3.4) between the third and fourth element of DataA.

Solutions

DataA = {{1.0, 4.0}, {3.0, 
          2.8}, {4, 2.5}, {5.0, 3.2}, {6.0, 3.9}, {10.0, 5.0}}

DataB = 
      Table[{x, x^2 + 10  x + 1}, {x, 0, 5, 1}]

DataC = 
      ReadList["c://documents and settings//jussi mattas//työp\
öytä//data.dat", {Number, Number}]

Note that ListPlot doesn't support the plotting of multiple sets of data points at once. Plotting one data set is done simply by ListPlot[data]. One way to plot them in one picture is to use the DisplayFunction option, which when given the value Identity suppresses the display of the resulting graphics object. The default value is $DisplayFunction, which has to be given as an option to the function Show to later display the graphics objects. Of course, plotting three data sets at once requires to distinguish them from oneanother, by color for example.

ListPlot[DataA, 
      PlotStyle -> {PointSize[0.01], Hue[0.0]}, 
      DisplayFunction -> Identity]

    ListPlot[DataB, PlotStyle -> {PointSize[0.01], Hue[0.3]}, 
      DisplayFunction -> Identity]

    ListPlot[DataC, PlotStyle -> {PointSize[0.01], Hue[0.6]}, 
      DisplayFunction -> Identity]

    Show[%, %%, %%%, DisplayFunction -> $DisplayFunction]

DataA = 
      Prepend[DataA, {0.0, 6.0}]

DataB = 
      Append[DataB, {6.0, 0.0}]

The function Insert[data,element,i] inserts element to be the ith element of the new list.

DataA = Insert[DataA, {2.0, 3.4}, 3]

Interpolating functions

First

1. Make a table that has as its elements ((Subscript[x, i]),  (Subscript[x, i]) exp(-(Subscript[x, i]))), where (Subscript[x, i]) = 0, 1/2, 1, 3/2.  .. , 5.  Assign the symbol A to the table (A=Table[...])
2. Print out the contents of table A using the TableForm[...] function.
3. Make an interpolating polynomial to table A. 
4. Calculate the value of the polynomial at points x = 0,  1,  2.5 and compare the results with the corresponding values in table A
5. Plot the interpolating polynomial and the function x (e^(-x)) into the same picture.
6. Plot the difference of the interpolating polynomial and the function x (e^(-x)).

Solutions

A = 
      Table[{x, x Exp[-x]}, {x, 0, 5, 0.5}]

TableForm[A]

poly = 
        InterpolatingPolynomial[A, x];

    Expand[%]

In the next excercise we'll do some tricks. The function Select takes as arguments a table and a criterion which it will apply to the table. In this case the criterion, #[[1]]==0.0& is in the form of a pure function, the argument of which, denoted by #, is a list. We test whether the first element of this list (which in this case is an element of A, equals to zero (or 1 or 2,5). The command & at the end just signifies that this really is a pure function. Select returns a list of elements of A which satisfy the criterion. We assume that there is only one element in A which satisfies our criterion and take the second element of the first element of the resulting list. This excercise can be done in a more simple way by just extracting A[[1,2]], A[[3,2]] and A[[6,2]].

Part[
          Select[A, #[[1]] == 0.0 &], 1, 2] - poly /. 
      x -> 0.0

    Part[Select[A, #[[1]] == 1.0 &], 1, 2] - poly /. 
      x -> 1.0

    Part[Select[A, #[[1]] == 2.5 &], 1, 2] - poly /. 
      x -> 2.5

ListPlot[A, 
        PlotStyle -> {PointSize[0.01], Hue[0.0]}, 
        DisplayFunction -> Identity];

    Plot[poly, {x, 0, 5}, 
        DisplayFunction -> Identity];

    Show[%, %%, DisplayFunction -> $DisplayFunction]

Plot[
      poly - x Exp[-x], {x, 0, 5}, 
      PlotRange -> {-6*10^(-6), 
          6*10^(-6)}]

Second

1. Read the file "data.dat" into table A.
2. Make an interpolating function to table A.
3. Plot the interpolating function.
4. Make an spline interpolating function to table A.
5. Plot the spline function.

Solutions

DataA = 
      ReadList["c://documents and settings//jussi mattas//työp\
öytä//data.dat", {Number,  Number}]

intf = 
      Interpolation[DataA]

Plot[
      intf[x], {x, 1, 6}]

The spline interpolation add-on package must be loaded before using it. There are three types of spline interpolating functions. We'll be primarily concerned of the cubic spline. Note the range of the interpolating spline. This is not the range of x or y, but of some parameter which gives the graph interpolating DataA as a parametric plot. It always starts from zero and the end point is the number of data points. Note that the usual graph {x,y[x]} can be presented as a parametric graph {x[t],y[t]} with x[t]=t.

<< 
      NumericalMath`SplineFit`

    splf = SplineFit[DataA, Cubic]

ParametricPlot[
      splf[x], {x, 0, 5}]

Third -- An evil example - this shows that you have to be careful with the interpolating polynomial!

1. Make a table ((Subscript[x, i]),  (1/(25 ((Subscript[x, i])^2) + 1))), where (Subscript[x, i]) = -1, -0.8,  ... ,  0.8,  1.
2. Make an interpolating polynomial to the table
3. Make an interpolating function to the table
4. Plot the interpolating function and polynomial in the same figure. What do you see? Try increseasing / decresing the step length when making the table. You can also try making a spline fit to the table. Try plotting that in the same figure also.

Solutions

A = 
        Table[{x, 1/(25  x^2 + 1)}, {x, -1, 1, 
            0.2}];

    poly = InterpolatingPolynomial[A, x];

    intf = Interpolation[A];

    Plot[{poly, intf[x]}, {x, -1, 1}, 
      PlotStyle -> {Hue[0], Hue[1/3]}]

It appears that the interpolating polynomial overshoots the function very badly, at least near the edges of the interpolation area. A finer spacing only makes the problem worse. A piecewise interpolation works much better for this case.

A = 
        Table[{x, 1/(25  x^2 + 1)}, {x, -1, 1, 
            0.1}];
poly = InterpolatingPolynomial[A, x];
intf = Interpolation[A];
Plot[{poly, intf[x]}, {x, -1, 1}]

Let's next compare the functions Interpolation and SplineFit:

<< NumericalMath`SplineFit`
splf = SplineFit[A, Cubic]
ParametricPlot[splf[x], {x, 0, 10}, PlotStyle -> {Hue[0]}, 
        DisplayFunction -> Identity];
Plot[intf[x], {x, -1, 1}, PlotStyle -> {Hue[2/3]}, 
        DisplayFunction -> Identity];
Show[%, %%, DisplayFunction -> $DisplayFunction]

Least squares fitting

Fit-command

1. Fit a second order polynomial to the data in the following table, and plot the data points and the fit.
({{-2, 2.84`}, {-1, 0.86`}, {0, 1.09`}, {1, 2.5100000000000002`}, {2, 7.17`}})

2. Fit the function a sin(x)  +  b cos(x)to the data in the following table, and plot the data points and the fit.
({{0, 1.04`}, {2, 0.44`}, {4, -1.51`}, {6, 0.6`}, {8, 0.78`}, {10, -1.45`}})

Solutions

Clear[x];

    A = {{-2, 2.84}, {-1, 0.86}, {0, 1.09}, {1, 2.51}, {2, 
            7.17}};

    fc = Fit[A, {1, x, x^2}, x];

    ListPlot[A, PlotStyle -> {Hue[0], PointSize[0.01]}, 
        DisplayFunction -> Identity];

    Plot[fc, {x, -2, 2}, PlotStyle -> Hue[1/3], 
        DisplayFunction -> Identity];

    Show[%, %%, DisplayFunction -> $DisplayFunction]

Just give the fitting functions as the second argument for Fit.

A = {{0, 1.04}, {2, 
            0.44}, {4, -1.51}, {6, 0.6}, {8, 
            0.78}, {10, -1.45}};

    fc = Fit[A, {Cos[x], Sin[x]}, x];

    ListPlot[A, PlotStyle -> {Hue[0], PointSize[0.01]}, 
        DisplayFunction -> Identity];

    Plot[fc, {x, 0, 10}, PlotStyle -> Hue[1/3], 
        DisplayFunction -> Identity];

    Show[%, %%, DisplayFunction -> $DisplayFunction]

Regress-command

1. Fit a straight line a x + b to following data, and get the standard error in a and b
({{-1, -0.5`}, {0, 1.2000000000000002`}, {1, 2.8000000000000003`}, {2, 4.9`}})

2. Fit the function a + b (e^x) + 
        c x (e^x) to following data, and get the standard error for a,  b and c
({{0, 1.6`}, {0.5`, 3.8000000000000003`}, {1.`, 6.300000000000001`}, {1.5`, 12.700000000000001`}, {2.`, 22.900000000000002`}})

Solutions

<< 
      Statistics`LinearRegression`

    Clear[x];

    A = {{-1, -0.5}, {0, 1.2}, {1, 2.8}, {2, 
            4.9}};

    fc = Fit[A, {1, x}, x];

    ParameterTable /. Regress[A, {1, x}, x];

    Print["Standard error for a: ", %[[1, 1, 2]]]

    Print["Standard error for b: ", %%[[1, 2, 2]]]

    ListPlot[A, PlotStyle -> {Hue[0], PointSize[0.02]}, 
        DisplayFunction -> Identity];

    Plot[fc, {x, -1, 2}, PlotStyle -> Hue[1/3], 
        DisplayFunction -> Identity];

    Show[%, %%, DisplayFunction -> $DisplayFunction]

At this point, an important notice: if you happen to make an error when handling a function which is not built-in but has to be loaded from a package, you might lose the usage of the function. For example, if Regress or ParameterTable ceases to work, you can reset them using the function Remove. Another thing that might seem weird is the use of three indices in obtaining the standard error. This is because we suppressed output when calculating the parameter table. Try different indices to see what's going on. The Print command simply prints its arguments. I have also used the function Fit because it's easier to extract the interpolating function from it than from Regress.

<< 
      Statistics`LinearRegression`

    A = {{0, 1.6}, {0.5, 3.8}, {1, 6.3}, {1.5, 12.7}, {2, 
            22.9}};

    fc = Fit[A, {1, Exp[x], x Exp[x]}, x];

    ParameterTable /. 
        Regress[A, {1, Exp[x], x Exp[x]}, x];

    Print["Standard error for a: ", %[[1, 1, 2]]]

    Print["Standard error for b: ", %%[[1, 2, 2]]]

    Print["Standard error for c: ", %%%[[1, 3, 
          2]]]

    ListPlot[A, PlotStyle -> {Hue[0], PointSize[0.02]}, 
        DisplayFunction -> Identity];

    Plot[fc, {x, 0, 2}, PlotStyle -> Hue[1/3], 
        DisplayFunction -> Identity];

    Show[%, %%, DisplayFunction -> $DisplayFunction]