Let's now try to fit a cubic polynomial to some data. Again, maple's built-in package can do this for us, but we like to work harder than we need to, so we understand more. First, let's get some data, using the provided function
cubdata := read(`/home/mat331/lsq_data`): plot(cubdata,style=point, axes=boxed);
As in the linear case, we will do the work "by hand".
To simplify typing, let's define the type of function we want to fit. We know (somehow) that our goal is to fit a cubic polynomial, so
cub := (x,a,b,c,d) -> a*x^3 + b*x^2 + c*x + d; 3 2 cub := (x, a, b, c, d) -> a x + b x + c x + d
Now we define an "error functional" which gives us the distance from our cubic to our data. Note that we use maple's "nops" command so that we don't need to explicitly say how many points there are.
G := (a,b,c,d) -> sum((cubdata(pts[i][1],a,b,c,d) - cubdata[i][2] )^2, i=1..nops(cubdata)); nops(cubdata) ----- \ 2 G := (a, b, c) -> ) (cub(cubdata[i][1], a, b, c,d) - cubdat[i][2]) / ----- i = 1
Now we just ask maple to find the (unique) critical point of G, which is our minimum. Since we used assign to set the coefficients, we won't see a result.
assign(solve( {diff( G(a,b,c,d), a)=0, diff( G(a,b,c,d), b)=0, diff( G(a,b,c,d), c)=0}, diff( G(a,b,c,d), d)=0}, {a,b,c,d} ));
but we can readily see what the answer is:
cub(x,a,b,c,d); 3 2 x - .4806519894 x - 10.94459406 x + 10.28284644
Finally, lets make a picture, just to see how well it fits
with(plots): cplot := plot(cub(x,a,b,c,d), x=-4..4, axes=boxed): pplot:= plot(cubdata, style=point): display({cplot, pplot});
Just for comparison, we can let maple do this for us, using the built-in least-squares package:
with(stats): fit[leastsquare[[x,y],y=A*x^3+B*x^2+C*x+D]] ([[cubdata[i][1]$i=1..nops(cubdata)], [cubdata[i][2] $i=1..nops(cubdata)]]); 3 2 y = x - .4806519892 x - 10.94459410 x + 10.28284642