403.sagews

#As per the discussion at http://mathforum.org/library/drmath/view/69525.html

var('a0 a1 a2 a3 sx sx2 sx3 sx4 sx5 sx6 sy sxy sx2y sx3y n')
eq1 = n*a0+sx*a1+sx2*a2+sx3*a3==sy
eq2 = sx*a0+sx2*a1+sx3*a2+sx4*a3==sxy
eq3 = sx2*a0+sx3*a1+sx4*a2+sx5*a3==sx2y
eq4 = sx3*a0+sx4*a1+sx5*a2+sx6*a3==sx3y

solve([eq1,eq2,eq3,eq4],a0,a1,a2,a3)

[[a0 == (-(sx2*(sx5^2 - sx4*sx6) + sx3^2*sx6 - 2*sx3*sx4*sx5 + sx4^3)*sy - sx3*(sx2*(-sx6*sxy - sx2y*sx5 + 2*sx3y*sx4) - sx4^2*sxy + sx*(sx3y*sx5 - sx2y*sx6)) - sx*((sx4*sx6 - sx5^2)*sxy + sx2y*sx4*sx5 - sx3y*sx4^2) - sx2*(sx4*sx5*sxy - sx2y*sx4^2) - sx3^2*(sx5*sxy + sx2y*sx4) - sx2^2*(sx2y*sx6 - sx3y*sx5) + sx3^3*sx3y)/(sx2*(n*(sx4*sx6 - sx5^2) - 2*sx*sx4*sx5) + sx^2*(sx5^2 - sx4*sx6) + sx3*(2*sx*sx2*sx6 + 2*n*sx4*sx5 + 2*sx2^2*sx5 + 2*sx*sx4^2) + sx3^2*(-n*sx6 - 2*sx*sx5 - 3*sx2*sx4) - sx2^3*sx6 - n*sx4^3 + sx2^2*sx4^2 + sx3^4), a1 == (-(sx*(sx4*sx6 - sx5^2) + sx3*(-sx2*sx6 - sx4^2) + sx2*sx4*sx5 + sx3^2*sx5)*sy - sx3*(sx2*(sx2y*sx4 - 2*sx5*sxy) + n*(sx2y*sx6 - sx3y*sx5) + sx*(sx2y*sx5 - sx3y*sx4)) - n*((sx5^2 - sx4*sx6)*sxy - sx2y*sx4*sx5 + sx3y*sx4^2) - sx2^2*(sx6*sxy - sx3y*sx4) - sx3^2*(sx4*sxy + sx2*sx3y) - sx*sx2*(sx3y*sx5 - sx2y*sx6) + sx2y*sx3^3)/(sx2*(n*(sx4*sx6 - sx5^2) - 2*sx*sx4*sx5) + sx^2*(sx5^2 - sx4*sx6) + sx3*(2*sx*sx2*sx6 + 2*n*sx4*sx5 + 2*sx2^2*sx5 + 2*sx*sx4^2) + sx3^2*(-n*sx6 - 2*sx*sx5 - 3*sx2*sx4) - sx2^3*sx6 - n*sx4^3 + sx2^2*sx4^2 + sx3^4), a2 == (-(sx3*(-sx*sx6 - sx2*sx5) + sx2^2*sx6 + sx*sx4*sx5 - sx2*sx4^2 + sx3^2*sx4)*sy - sx3*(n*(sx6*sxy - sx3y*sx4) + sx*(sx5*sxy - 2*sx2y*sx4) + sx2*sx4*sxy - sx2^2*sx3y) - sx2*(sx*(sx3y*sx4 - sx6*sxy) + n*(sx3y*sx5 - sx2y*sx6)) - n*(sx2y*sx4^2 - sx4*sx5*sxy) + sx3^3*sxy - sx^2*(sx2y*sx6 - sx3y*sx5) - sx3^2*(sx*sx3y + sx2*sx2y))/(sx2*(n*(sx4*sx6 - sx5^2) - 2*sx*sx4*sx5) + sx^2*(sx5^2 - sx4*sx6) + sx3*(2*sx*sx2*sx6 + 2*n*sx4*sx5 + 2*sx2^2*sx5 + 2*sx*sx4^2) + sx3^2*(-n*sx6 - 2*sx*sx5 - 3*sx2*sx4) - sx2^3*sx6 - n*sx4^3 + sx2^2*sx4^2 + sx3^4), a3 == ((sx3*(-sx*sx5 - 2*sx2*sx4) + sx2^2*sx5 + sx*sx4^2 + sx3^3)*sy + sx3*(n*(sx5*sxy + sx2y*sx4) + sx*sx4*sxy + 2*sx*sx2*sx3y + sx2^2*sx2y) + sx2*(sx*(-sx5*sxy - sx2y*sx4) + n*(sx3y*sx4 - sx2y*sx5)) + sx3^2*(-sx2*sxy - n*sx3y - sx*sx2y) - n*sx4^2*sxy + sx2^2*sx4*sxy + sx^2*(sx2y*sx5 - sx3y*sx4) - sx2^3*sx3y)/(sx2*(n*(sx4*sx6 - sx5^2) - 2*sx*sx4*sx5) + sx^2*(sx5^2 - sx4*sx6) + sx3*(2*sx*sx2*sx6 + 2*n*sx4*sx5 + 2*sx2^2*sx5 + 2*sx*sx4^2) + sx3^2*(-n*sx6 - 2*sx*sx5 - 3*sx2*sx4) - sx2^3*sx6 - n*sx4^3 + sx2^2*sx4^2 + sx3^4)]]

var ('a0 a1 a2 a3 x')
model(x) = a0+a1*x+a2*x^2+a3*x^3
data = [ [0,0],[1,1],[2,3],[3,15],[4,17],[6,19] ]
fit=find_fit(data, model, solution_dict=True)
print fit
model_fit=model.subs(fit)
plot(model_fit, (x, -1, 6))+plot(point(data))

{a1: -2.8359213250600614, a0: 0.025879917182081602, a3: -0.43115942029297494, a2: 3.5828157349952612}