Velocity Reviews > Symbolic Polynomial Interpolator in C

# 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/

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

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."