Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Symbolic Polynomial Interpolator in C

Reply
Thread Tools

Symbolic Polynomial Interpolator in C

 
 
websnarf@gmail.com
Guest
Posts: n/a
 
      08-10-2005
http://www.velocityreviews.com/forums/(E-Mail Removed) wrote:
> I am looking for a quick C program that takes n+1 pairs of values
> (integers) (a_i, f(a_i)), i=0,...,n, generates the coefficients
> \alpha_i, i=0,...,n of the polynomial of degree n that fits these
> points and then outputs the polynomial with indeterminant in the form
> \sum_i=0^n\alpha_i X^i, where X is just some indeterminant.


In general, you need an n+1 degree polynomial to fit n points. You are
specifying the points as 0, ..., n which is already n+1 pointts. I'm
going to assume you mean to specify the points as 0, ..., n-1.

There is a well known interpolation mechanism called "Lagrange
Interpolation":

Let LF(X) = \prod_i=0^n-1\(X - a_i), and let F_i(X) = LF(X)/(X - a_i),
where F_I(a_i) is calculated by taking the limit as X->a_i; i.e.,
cancel out the factors.

Then the polynomial you are looking for is:

\sum_i=0,n-1\f(a_i)*F_i(X)/F_i(a_i)

The reason is that F_i(a_j) = 0 for any i not equal to j.

So then your problem reduces to how do you write a program to expand
the polynomials: f(a_i)/F_i(a_i) * F_i(X)

The way you would want start is to break this down into its
coefficients by writing a function that computed:

double coefficient_F(int i, int p, double * a, int n);

which is the pth coefficient of F_i as defined on the vector a of size
n. The idea is that you would do it recursively, by peeling off the
last factor one at a time:

double coefficient_F(int i, int p, double * a, int n) {
if (p == n-(i < n)) return 1.0;
if (p < n) return 0.0;
if (i = n-1) n--;
return coefficient_F(i, p-1, a, n-1)
- a[n-1] * coefficient_F(i, p, a, n-1);
}

Notice how the first termination condition basically looks for when you
have exactly p factors of (X - a_*) in this factor, the second just
deals with the case where you asked for coefficient that was negative,
or otherwise too low, and otherwise, it performs the one recursive
expansion corresponding to the last factor in the polynomial.

Then you need the actually need to be able to compute value F_i(a_i)
which is just a simple for loop (I'll leave that as an exercise for
you.)

Then you have what you need to perform a polynomial expansion of each
term of the sum, then you need to just sum the coefficients.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

 
Reply With Quote
 
 
 
 
Eric Sosman
Guest
Posts: n/a
 
      08-11-2005
(E-Mail Removed) wrote:
> (E-Mail Removed) wrote:
>
>>I am looking for a quick C program that takes n+1 pairs of values
>>(integers) (a_i, f(a_i)), i=0,...,n, generates the coefficients
>>\alpha_i, i=0,...,n of the polynomial of degree n that fits these
>>points and then outputs the polynomial with indeterminant in the form
>>\sum_i=0^n\alpha_i X^i, where X is just some indeterminant.

>
>
> In general, you need an n+1 degree polynomial to fit n points.


ITYM a degree n polynomial to fit n+1 points.

(Consider an easy case: To fit two points you need a
straight line, that is, a 1-degree polynomial. Using a
cubic would be redundantly superfluous overkill ...)

--
Eric Sosman
(E-Mail Removed)lid
 
Reply With Quote
 
 
 
 
Julian V. Noble
Guest
Posts: n/a
 
      08-12-2005
(E-Mail Removed) wrote:
>
> I am looking for a quick C program that takes n+1 pairs of values
> (integers) (a_i, f(a_i)), i=0,...,n, generates the coefficients
> \alpha_i, i=0,...,n of the polynomial of degree n that fits these
> points and then outputs the polynomial with indeterminant in the form
> \sum_i=0^n\alpha_i X^i, where X is just some indeterminant.
>
> Any hints?
>
> Thanks


Use Gram polynomials, which for integers in a given range (that is,
equally spaced x-values) are actually Chebyshev polynomials. See
Ralston, Intro to Numerical Analysis.

The chief virtue of Gram polynomials is you don't have to invert
matrices. The disadvantage is that the result is always in the
form

f\left( x \right) \approx
\sum\limits_{k = 1}^m {a_k C_k \left( x \right)}

--
Julian V. Noble
Professor Emeritus of Physics
(E-Mail Removed)
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/

"For there was never yet philosopher that could endure the
toothache patiently."

-- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1.
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Linear Interpolator aDawg VHDL 0 06-23-2006 04:48 PM
Looking for triangulator/interpolator Grant Edwards Python 13 05-28-2006 03:29 AM
polynomial S.Muralidharan VHDL 1 11-04-2004 12:44 PM
is there a inverse polynomial for the CRC generator polynomial Chen L. C Programming 2 07-06-2004 01:23 AM
polynomial division remainder Manfred Balik VHDL 5 05-18-2004 02:37 AM



Advertisments