# 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 > restart: > read(`/home/mat331/lsq_data`); defined line_pts(), quadratic_pts(), cubic_pts(), and circle_pts() > cubdata := cubic_pts(): > 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((cub(cubdata[i][1],a,b,c,d) - cubdata[i][2] > )^2, i=1..nops(cubdata)); G := (a, b, c, d) -> nops(cubdata) ----- \ 2 ) (cub(cubdata[i][1], a, b, c, d) - cubdata[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 -.3688032502 x + .3920904491 x + 2.656274173 x + 1.999084485 # Finally, lets make a picture, just to see how well it fits > with(plots): > cplot:= plot(cub(x,a,b,c,d), x=-3..3, 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: > > 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 = -.3688032508 x + .3920904505 x + 2.656274177 x + 1.999084480 # they match to 8 figures, which is all we could have hoped for. But of # course, they should.