• Choose a category
• All
• Classic-Fiction

By Root 1395 0

This is so that we are ready to compute what equates to the third row of the divided-difference table. (Nothing needs to be done for the first two rows because these entries are already stored in x and fx.) The final initialization step is to store the value of fx[0] in coeff[0] since this is the first coefficient of the interpolating polynomial.

The process of computing divided differences revolves around a single nested loop, which uses the formula for divided differences discussed earlier in the chapter. In terms of Figure 13.2 and Figure 13.3, the outer loop, k, counts the number of rows for which entries must be computed (excluding the rows for x and fx). The inner loop, i, controls which entry is being computed in the current row. As we complete the entries in each row, the value in table [0] becomes the next coefficient for the interpolating polynomial. Thus, we store this value in coeff[k]. Once we have determined all coefficients for the interpolating polynomial, we evaluate the polynomial at each point in z. The results are stored in pz.

The runtime complexity of interpol is O (mn 2), where m is the number of values in z (and values returned in pz), and n is the number of values in x (and fx). The factor n 2 comes from the following. For each iteration of the loop controlled by j, we multiply one factor more than the previous term into the current term. Thus, when j is 1, term requires one multiplication; when j is 2, term requires two multiplications, and so forth until when j is n - 1, term requires n - 1 multiplications. Effectively, this becomes a summation from 1 to n - 1, which results in a running time of T (n) = (n (n + 1)/2) - n, times some constant amount of time. (This is from the well-known formula for summing an arithmetic series from 1 to n.) In O-notation, this simplifies to O (n 2). The factor m in O (mn 2) comes from evaluating the interpolating polynomial once for each point in z. The first nested loop, in which divided differences are computed, is O (n 2). Thus, this term is not significant relative to mn 2, which has the additional factor m.

Example 13.1. Implementation of Polynomial Interpolation

/*****************************************************************************

* *

* ----------------------------- interpol.c ---------------------------------*

* *

*****************************************************************************/

#include

#include

#include "nummeths.h"

/*****************************************************************************

* *

* ------------------------------ interpol -------------------------------- *

* *

*****************************************************************************/

int interpol(const double *x, const double *fx, int n, double *z, double

*pz, int m) {

double term,

*table,

*coeff;

int i,

j,

k;

/*****************************************************************************

* *

* Allocate storage for the divided-difference table and coefficients. *

* *

*****************************************************************************/

if ((table = (double *)malloc(sizeof(double) * n)) == NULL)

return -1;

if ((coeff = (double *)malloc(sizeof(double) * n)) == NULL) {

free(table);

return -1;

}

/*****************************************************************************

* *

* Initialize the coefficients. *

* *

*****************************************************************************/

memcpy(table, fx, sizeof(double) * n);

/*****************************************************************************

* *

* Determine the coefficients of the interpolating polynomial. *

* *

*****************************************************************************/

coeff[0] = table[0];

for (k = 1; k < n; k++) {

for (i = 0; i < n - k; i++) {

j = i + k;

table[i] = (table[i + 1] - table[i]) / (x[j] - x[i]);

}

coeff[k] = table[0];

}

free(table);

/*****************************************************************************

* *

* Evaluate the interpolating