Polynomial Curve Fitting with Least Squares Method

This worksheet is to teach you how to use the least squares method for finding the optimal polynomial

function for a given data set. In Example 1 we are given that the polynomial should be a cubic.

Then we plot the data and the polynomial. In Example 2 we are not given the degree of the polynomial to use. Thus, we plot the data and guess what degree we should use. Finally there are some problems for you to do.

Example 1.

>	with(LinearAlgebra):with(plots):

Below we are give 8 data points in the plane. The x coordinates are the integers from 0 to 7 and the y coordinates are various real numbers. We are to find the best cubic polynomial y=a+bx+cx^2+dx^3 that fits this data.

>	xdata:=Vector([0, 1, 2, 3, 4, 5, 6, 7]): ydata:=Vector([0.2,1.2,3.2,1.0,2.0,0.1,2.0,4.0]):

For each data point we get an equation. For example (2,3.2) yields a+2b+4c+8d = 3.2. We have eight equations and 4 unknowns.  

(Change the colons (:) to semicolons (;) and see what happens.

We shall use the VandermondeMatrix function to create the Vandermonde matrix A. The 8 and 4 determine the size of the Vandermonde matrix we want.

>	A:=VandermondeMatrix(xdata,8,4);

You can check that Ax=ydata has no solutions. Thus, this is a job for the least squares method!

We need to solve   A^T * A * x = A^T * ydata,   for x=[a,b,c,d]^T. (I used coeffvector for x below because I did not want to confuse it with xdata.)

>	coeffvector:=LinearSolve(Transpose(A).A,Transpose(A).ydata);

This tells us the values for the four coefficients of the cubic, a=0.00909091, etc. Next we plot th cubic and

the original data.

There are three commands below that will be executed together. The first plots the polynomial but stores it as CurvePlot rather than display it. The second plots the data, but due to a quirk in Malpe the data type must be changed from Vector to list. The convert cammond does this.

>	CurvePlot:=plot(coeffvector[1]+coeffvector[2]*x+coeffvector[3]*x^2+coeffvector[4]*x^3, x=0..7,color=green,thickness=3): PointPlot:=plot([[convert(xdata,list)[n],convert(ydata,list)[n]] $n=1..8],style=POINT): display([PointPlot,CurvePlot]);

Note: There is a direct way to do least squares curve fitting using the LeastSquares cammand.  

It is easier to use but less instructive. You can use if you like, but it is up to you to figure out how it works.

Example 2. This is similar to Example 1, but you are not given the degree of the polynomial to use.

Also, the x coordinates are not evenly spaced.

>	restart;with(LinearAlgebra):with(plots):

>	xdata:=Vector([0.2,1.2,2.1,3.4,5.1,7.5,10.1,11.2]): ydata:=Vector([0.3,0.6,1.7,2.5,2.3,1.9,1.4,1.1]):

We plot the data points.

>	PointPlot:=plot([[convert(xdata,list)[n],convert(ydata,list)[n]] $n=1..8],style=POINT): display(PointPlot);

It is certainly not a line. Looks more or less parabolic, so let's try f(x)=ax^+bx+c.

>	A:=VandermondeMatrix(xdata,8,3);

>	coeffvector:=LinearSolve(Transpose(A).A,Transpose(A).ydata);

>	CurvePlot:=plot(coeffvector[1]+coeffvector[2]*x+coeffvector[3]*x^2, x=0..12,color=green,thickness=3): PointPlot:=plot([[convert(xdata,list)[n],convert(ydata,list)[n]] $n=1..8],style=POINT): display([PointPlot,CurvePlot]);

Problem 1. Find a polynomial that best fits this data set:

 {(-4.0, -5.9), (-3.7, -3.5), (-3.2, -0.5), (-2.0, 3.8), (-0.9, 6.5), (0.2, 6.4), (1.7, 5.2), (2.5, 4.1), (4.6, 1.5), (5.4, 1.2), (5.7, 1.3)}.

Graph the data points and your polynomial on one graph. Delete the two examples and print your results.

Be sure to type in your name.