# 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 GNU LGPL license. // // Modified: // // 08 July 2015 // // Author: // // John Burkardt // // Reference: // // Carl deBoor, // A Practical Guide to Splines, // Springer, 2001, // ISBN: 0387953663, // LC: QA1.A647.v27. // // Parameters: // // Input, int ND, the number of data points. // // Input, double XD[ND], the X values at which data was taken. // These values must be distinct. // // Input, 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 GNU LGPL license. // // Modified: // // 08 July 2015 // // Author: // // John Burkardt // // Reference: // // Carl deBoor, // A Practical Guide to Splines, // Springer, 2001, // ISBN: 0387953663, // LC: QA1.A647.v27. // // Parameters: // // Input, int ND, the order of the difference table. // // Input, double XD[ND], the X values of the difference table. // // Input, double CD[ND], the divided differences. // // Input, int NI, the number of interpolation points. // // Input, double XI[NI], the interpolation points. // // Output, double NEWTON_VALUE_1D[NI], the interpolation 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; } //****************************************************************************80 double *r8mat_copy_new ( int m, int n, double a1[] ) //****************************************************************************80 // // Purpose: // // R8MAT_COPY_NEW copies one R8MAT to a "new" R8MAT. // // Discussion: // // An R8MAT is a doubly dimensioned array of R8's, which // may be stored as a vector in column-major order. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 03 July 2008 // // Author: // // John Burkardt // // Parameters: // // Input, int M, N, the number of rows and columns. // // Input, double A1[M*N], the matrix to be copied. // // Output, double R8MAT_COPY_NEW[M*N], the copy of A1. // { double *a2; int i; int j; a2 = new double[m*n]; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a2[i+j*m] = a1[i+j*m]; } } return a2; } //****************************************************************************80 double *r8vec_copy_new ( int n, double a1[] ) //****************************************************************************80 // // Purpose: // // R8VEC_COPY_NEW copies an R8VEC to a new R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 03 July 2008 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of entries in the vectors. // // Input, double A1[N], the vector to be copied. // // Output, double R8VEC_COPY_NEW[N], the copy of A1. // { double *a2; int i; a2 = new double[n]; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } //****************************************************************************80 double *r8vec_linspace_new ( int n, double a_first, double a_last ) //****************************************************************************80 // // Purpose: // // R8VEC_LINSPACE_NEW creates a vector of linearly spaced values. // // Discussion: // // An R8VEC is a vector of R8's. // // 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. // // In other words, the interval is divided into N-1 even subintervals, // and the endpoints of intervals are used as the points. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 March 2011 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of entries in the vector. // // Input, double A_FIRST, A_LAST, the first and last entries. // // Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. // { double *a; int i; a = new double[n]; if ( n == 1 ) { a[0] = ( a_first + a_last ) / 2.0; } else { for ( i = 0; i < n; i++ ) { a[i] = ( ( double ) ( n - 1 - i ) * a_first + ( double ) ( i ) * a_last ) / ( double ) ( n - 1 ); } } return a; } //****************************************************************************80 double r8vec_max ( int n, double r8vec[] ) //****************************************************************************80 // // Purpose: // // R8VEC_MAX returns the value of the maximum element in an R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 22 August 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of entries in the array. // // Input, double R8VEC[N], a pointer to the first entry of the array. // // Output, double R8VEC_MAX, the value of the maximum element. This // is set to 0.0 if N <= 0. // { int i; double value; value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( value < r8vec[i] ) { value = r8vec[i]; } } return value; } //****************************************************************************80 double r8vec_min ( int n, double r8vec[] ) //****************************************************************************80 // // Purpose: // // R8VEC_MIN returns the value of the minimum element in an R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 02 July 2005 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of entries in the array. // // Input, double R8VEC[N], the array to be checked. // // Output, double R8VEC_MIN, the value of the minimum element. // { int i; double value; value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( r8vec[i] < value ) { value = r8vec[i]; } } return value; } //****************************************************************************80 double r8vec_norm_affine ( int n, double v0[], double v1[] ) //****************************************************************************80 // // Purpose: // // R8VEC_NORM_AFFINE returns the affine L2 norm of an R8VEC. // // Discussion: // // The affine vector L2 norm is defined as: // // R8VEC_NORM_AFFINE(V0,V1) // = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2 ) // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 27 October 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the dimension of the vectors. // // Input, double V0[N], the base vector. // // Input, double V1[N], the vector. // // Output, double R8VEC_NORM_AFFINE, the affine L2 norm. // { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + ( v1[i] - v0[i] ) * ( v1[i] - v0[i] ); } value = sqrt ( value ); return value; } //****************************************************************************80 void r8vec_print ( int n, double a[], string title ) //****************************************************************************80 // // Purpose: // // R8VEC_PRINT prints an R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 16 August 2004 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of components of the vector. // // Input, double A[N], the vector to be printed. // // Input, string TITLE, a title. // { int i; cout << "\n"; cout << title << "\n"; cout << "\n"; for ( i = 0; i < n; i++ ) { cout << " " << setw(8) << i << ": " << setw(14) << a[i] << "\n"; } return; } //****************************************************************************80 void r8vec2_print ( int n, double a1[], double a2[], string title ) //****************************************************************************80 // // Purpose: // // R8VEC2_PRINT prints an R8VEC2. // // Discussion: // // An R8VEC2 is a dataset consisting of N pairs of real values, stored // as two separate vectors A1 and A2. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2002 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of components of the vector. // // Input, double A1[N], double A2[N], the vectors to be printed. // // Input, string TITLE, a title. // { int i; cout << "\n"; cout << title << "\n"; cout << "\n"; for ( i = 0; i <= n - 1; i++ ) { cout << setw(6) << i << ": " << setw(14) << a1[i] << " " << setw(14) << a2[i] << "\n"; } return; } //****************************************************************************80 void timestamp ( ) //****************************************************************************80 // // Purpose: // // TIMESTAMP prints the current YMDHMS date as a time stamp. // // Example: // // 31 May 2001 09:45:54 AM // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 08 July 2009 // // Author: // // John Burkardt // // Parameters: // // None // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct std::tm *tm_ptr; std::time_t now; now = std::time ( NULL ); tm_ptr = std::localtime ( &now ); std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ); std::cout << time_buffer << "\n"; return; # undef TIME_SIZE }