# include # include # include # include # include # include using namespace std; # include "newton_interp_1d.hpp" //****************************************************************************80 double *newton_coef_1d ( int nd, double xd[], double yd[] ) //****************************************************************************80 // // Purpose: // // newton_coef_1d() computes coefficients of a Newton 1D interpolant. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 08 July 2015 // // Author: // // John Burkardt // // Reference: // // Carl deBoor, // A Practical Guide to Splines, // Springer, 2001, // ISBN: 0387953663, // LC: QA1.A647.v27. // // Input: // // int ND, the number of data points. // // double XD[ND], the X values at which data was taken. // These values must be distinct. // // double YD[ND], the corresponding Y values. // // Output: // // double NEWTON_COEF_1D[ND], the divided difference coefficients. // { double *cd; int i; int j; // // Copy the data values. // cd = new double[nd]; for ( i = 0; i < nd; i++ ) { cd[i] = yd[i]; } // // Make sure the abscissas are distinct. // for ( i = 0; i < nd; i++ ) { for ( j = i + 1; j < nd; j++ ) { if ( xd[i] - xd[j] == 0.0 ) { cerr << "\n"; cerr << "newton_coef_1d(): Fatal error!\n"; cerr << " Two entries of XD are equal!\n"; cerr << " XD[" << i << "] = " << xd[i] << "\n"; cerr << " XD[" << j << "] = " << xd[j] << "\n"; exit ( 1 ); } } } // // Compute the divided differences. // for ( i = 1; i <= nd - 1; i++ ) { for ( j = nd - 1; i <= j; j-- ) { cd[j] = ( cd[j] - cd[j-1] ) / ( xd[j] - xd[j-i] ); } } return cd; } //****************************************************************************80 double *newton_value_1d ( int nd, double xd[], double cd[], int ni, double xi[] ) //****************************************************************************80 // // Purpose: // // newton_value_1d() evaluates a Newton 1D interpolant. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 08 July 2015 // // Author: // // John Burkardt // // Reference: // // Carl deBoor, // A Practical Guide to Splines, // Springer, 2001, // ISBN: 0387953663, // LC: QA1.A647.v27. // // Input: // // int ND, the order of the difference table. // // double XD[ND], the X values of the difference table. // // double CD[ND], the divided differences. // // int NI, the number of interpolation points. // // double XI[NI], the interpolation points. // // Output: // // double NEWTON_VALUE_1D[NI], the interpolated values. // { int i; int j; double *yi; yi = new double[ni]; for ( j = 0; j < ni; j++ ) { yi[j] = cd[nd-1]; for ( i = 2; i <= nd; i++ ) { yi[j] = cd[nd-i] + ( xi[j] - xd[nd-i] ) * yi[j]; } } return yi; }