Thursday, February 13, 2014

La Grange and Neville's method implementation to fit an nth degree polynomial to a set of points


In numeric analysis class I had to create a solution for some curve fitting algorithms.  I choose to use C as well as implementing these functions in excel, because when I would like to automate data processing I would rather use a small program to perform the task as quickly as possible, rather than a manual process using Excel.

These functions can fit a curve to hundreds of points in hundredths of a second on a relatively slow net book.  The reason you would do this is that you may not know the function that governs a system, but you need to have a smooth continuous function that gives an accurate y value when you input an x value.  These curve fitting methods allow you to model an unknown function, inside the range of your data points, calculating the nth degree polynomial for your n points you know.  

The first is the Lagrange method, which builds a polynomial one term at a time, calling a second function to construct the parts of the polynomial for a specific x.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Lagrange.c

The second one was Neville's method and I didn't understand how this was working until I manually did a small problem using excel to calculate the result.  I finally worked it all out and came up with this very simple solution:

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Neville.c

Unfortunately, as easy to understand as the above program is, it has a serious problem where it recalculates the sub results many times.  This results in taking many minutes for relatively small point sets.

So I re-implemented the solution in the below code to build the solution from the "base of the triagle" rather than from the apex, so that each part of the solution is only calculated once.  This allows things that took minutes to run in sub second time intervals.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Neville-mod.c

The only tricky thing about this function is that I had to pass in the original set of points, unchanged, then use the depth that we are at in order to index the correct Point[i].x.  I tried just pushing a single value up, but the x values are selected diagonally so a single calculated value could need one of two returns for f(x).

The depth is really only an indication of the spread of the x points we need for the current calculation.

A benefit of this solution is that you can print out all the results of the sub functions and the data you used to generate the result, but uncommenting the printfs in this section.

This function also only calls itself n-1 times, and calculates just n-depth sub results in each call.
This value is stored in a net set of points and passed into the next lower depth of the function as it calls itself.

The function returns and unwinds from all the calls when it calls itself with just a single point.  All we do is return the value we last calculated and we are done. :D

No comments:

Post a Comment