   Polynomial Function Fitting
  Hi,Does anyone remember those assignments we did in Calculus, where you start with a list of givens, and find a cubic, quadratic, etc. equation that fits them?  Like, find a, b, c, and d of the cubic equation f(x)=ax3+bx2+cx+d that satisfies f(0)=0, f(1)=1, f'(0)=1, f'(1)=0?I run into these problems quite a bit while programming, and rather than work them out by hand, I was wondering if someone has written a program somewhere to expedite finding the coefficients--or do I have to go and write it myself?Thanks,-G   
  You can certainly use software like Mathematica to solve this type of problem. You can also do the work of turning the problem into a system of linear equations, express it in matrix format and invert the matrix using some online calculator.;  I'd just use MATLAB (since I have a license) or Octave (if I didn't).  It has a function called "polyfit" which will give you a least-squares polynomial approximation to some given data.  Naturally, when you ask it to give you a polynomial whose degree is the same as the number of data points you have, this reduces to exact interpolation.If I wanted to write some code to do it myself, I'd probably just use Lagrange form and never solve for the coeffients at all.  If I needed the coefficients, I'd probably just do the naive thing and construct the Vandemonde matrix and invert it, though there are some algorithms that can do it faster than this.  And as usual, you can probably learn everything you need to about the subject just by reading the Wikipedia article.Actually, I'll be honest; I'd probably not use a single interpolating polynomial at all, and go with bsplines instead.  But now I'm getting off-topic...   ;  So I went ahead and wrote the following script (Python)import numpy as nppoints = [[0,0],[0.1,1],[1,0]] #Of the form [x,f (x)]derivatives = [[0.1,0,1],[1,0,1]] #Of the form [x,f^(n) (x),n]order = len(points)+len(derivatives)-1equations_coefficients = []equations_constants = []for p in points:    coeffs_str = ""    coefficients = []    for term in xrange(order+1):   coefficients.append(p[0]**(order-term))   coeffs_str += str(coefficients[-1])+"*x"*(order-term)+"+"    equations_coefficients.append(coefficients)    equations_constants.append([p[1]])##    print str(p[1])+" = "+coeffs_str[:-1]def power_rule(const,exp,num):    if num == 1:   if exp == 0:  return [0,0]   else:  return [const*exp,exp-1]    else:   derivative = power_rule(const,exp,1)   return power_rule(derivative[0],derivative[1],num-1)for d in derivatives:    d_num = d[2]    coefficients = []    for term in xrange(order+1):   term = power_rule(1,order-term,d_num)   coefficients.append(term[0]*(d[0]**term[1]))    equations_coefficients.append(coefficients)    equations_constants.append([d[1]])for p in points:    print "f("+str(p[0])+") = "+str(p[1])for d in derivatives:    print "f"+"'"*d[2]+"("+str(d[0])+") = "+str(d[1])print "Order",order##print "Coefficients:"##print equations_coefficients##print "Constants:"##print equations_constantsA = np.matrix( equations_coefficients )x = np.matrix( equations_constants )cs = []for coeff in np.linalg.solve(A, x).tolist():    cs.append(str(coeff[0]))equation = "f(x) = "for term in xrange(order+1):    xstring = ""    for x in xrange(order-term):   xstring += "*x"    equation += cs[term]+xstring    if term != order:   equation += " + "print equationIt treats the equation coefficients as the variables to be solved for.  The result of the current parameters isf(x) = -96.0219478738*x*x*x*x + 213.991769547*x*x*x + -139.917695473*x*x + 21.9478737997*x + 0.0. . . which, when graphed here looks wrong.  Although the equation satisfies all conditions, it doesn't look the way I want it to (which in this case is akin to a strongly right skewed chi-squared distribution).  It seems that f'(x) is getting some extra zeros.  What I'd really like is for f(x) and f'(x) to have only the zeros specified . . .See where I'm going with this?  Anyone know how to get the right function?Thanks,-G   ;  Well, if we have a polynomial f of degree 4.f(x) = ax^4 + bx^3 + cx^2 + dx + eIts derivative is:f'(x) = 4ax^3 + 3b^2 + 2cx + dIf you know points xi( i=1,...,5 ) for which you have:f(xi) = 0orf'(xi) = 0you can solve your system easily. If you have 5 equations of that type, you just replace x, f(x) and f'(x) by their values and you obtain 5 linear equations with 5 unknowns.P.S:* of course, you're searching for a polynomial of degree 3...So it's even faster to find the polynomial.* what you want to achieve is what is called "polynomial interpolation". Since it is an interpolation, it has its pitfalls.[Edited by - johnstanp on March 15, 2010 6:32:39 AM];  What you are encountering is an effect known as "over-fitting."  In essence, a degree 4 polynomial fitting your data points is too much.  Although you have 4 points, you actually want  to fit a smaller (or simpler) function to your data-set.  The easy way to do this is to fit a 2-d or 3-d polynomial instead of a 4-d one, although since you said the shape you want isn't a polynomial at all, maybe you should use that basis instead.To fit an m-d polynomial to n data points, where n > m, you do the exact same linear algebra you do to fit the square case (a square vandemonde matrix), except you end up with an OVERDETERMINED system of linear equations.  For example,  fitting the best 1-d (line) polynomial y=mx+b to a series of 4 points gives you[y1] [x1 1][m][y2] [x2 1][y3]=[x3 1][y4] [x4 1][y5] [x5 1]You can solve this system y=Ax using LEAST SQUARES, which, without going into proof, is basically solving the new system A'y=(A'A)x.;  Hi,Thanks; that sounds like a plausible solution.  I'm familiar with statistical methods, but I'm still not sure what you mean by solving matrices using least squares.  Perhaps I'm missing something.Thanks,-G   ; Least squares.   