Sign in to follow this  
Geometrian

Polynomial Function Fitting

Recommended Posts

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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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...

Share this post


Link to post
Share on other sites
So I went ahead and wrote the following script (Python)
import numpy as np
points = [[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)-1

equations_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_constants

A = 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 equation
It treats the equation coefficients as the variables to be solved for. The result of the current parameters is
f(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

Share this post


Link to post
Share on other sites
Well, if we have a polynomial f of degree 4.

f(x) = ax^4 + bx^3 + cx^2 + dx + e

Its derivative is:

f'(x) = 4ax^3 + 3b^2 + 2cx + d

If you know points xi( i=1,...,5 ) for which you have:

f(xi) = 0
or
f'(xi) = 0

you 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]

Share this post


Link to post
Share on other sites
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][b]
[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.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this