# include # include # include # include # include # include # include "barycentric_interp_1d.h" # include "r8lib.h" /******************************************************************************/ double *lagcheby1_interp_1d ( int nd, double xd[], double yd[], int ni, double xi[] ) /******************************************************************************/ /* Purpose: LAGCHEBY1_INTERP_1D evaluates the Lagrange Chebyshev 1 interpolant. Discussion: The weight vector WD computed below is only valid if the data points XD are, as expected, the Chebyshev Type 1 points for [-1,+1], or a linearly mapped version for [A,B]. The XD values may be computed by: xd = r8vec_cheby1space ( nd, a, b ) for instance. Thanks to John Ferrier for pointing out that DENOM and NUMER needed to be initialized, 16 September 2013. Licensing: This code is distributed under the MIT license. Modified: 16 September 2013 Author: John Burkardt Reference: Jean-Paul Berrut, Lloyd Trefethen, Barycentric Lagrange Interpolation, SIAM Review, Volume 46, Number 3, September 2004, pages 501-517. Parameters: Input, int ND, the number of data points. ND must be at least 1. Input, double XD[ND], the data points. Input, double YD[ND], the data values. Input, int NI, the number of interpolation points. Input, double XI[NI], the interpolation points. Output, double LAGCHEBY1_INTERP_1D[NI], the interpolated values. */ { double *denom; int *exact; int i; int j; double *numer; double pi = 3.141592653589793; double t; double theta; double wd; double *yi; denom = ( double * ) malloc ( ni * sizeof ( double ) ); exact = ( int * ) malloc ( ni * sizeof ( int ) ); numer = ( double * ) malloc ( ni * sizeof ( double ) ); yi = ( double * ) malloc ( ni * sizeof ( double ) ); for ( i = 0; i < ni; i++ ) { exact[i] = -1; denom[i] = 0.0; numer[i] = 0.0; } for ( j = 0; j < nd; j++ ) { theta = ( double ) ( 2 * j - 1 ) * pi / ( double ) ( 2 * nd ); wd = r8_mop ( j + 1 ) * sin ( theta ); for ( i = 0; i < ni; i++ ) { if ( xi[i] == xd[j] ) { exact[i] = j; numer[i] = yd[j]; denom[i] = 1.0; } if ( exact[i] == -1 ) { t = wd / ( xi[i] - xd[j] ); numer[i] = numer[i] + t * yd[j]; denom[i] = denom[i] + t; } } } for ( i = 0; i < ni; i++ ) { yi[i] = numer[i] / denom[i]; } return yi; } /******************************************************************************/ double *lagcheby2_interp_1d ( int nd, double xd[], double yd[], int ni, double xi[] ) /******************************************************************************/ /* Purpose: LAGCHEBY2_INTERP_1D evaluates the Lagrange Chebyshev 2 interpolant. Discussion: The weight vector WD computed below is only valid if the data points XD are, as expected, the Chebyshev Type 2 points for [-1,+1], or a linearly mapped version for [A,B]. The XD values may be computed by: xd = r8vec_cheby2space ( nd, a, b ) for instance. Thanks to John Ferrier for pointing out that DENOM and NUMER needed to be initialized, 16 September 2013. Licensing: This code is distributed under the MIT license. Modified: 16 September 2013 Author: John Burkardt Reference: Jean-Paul Berrut, Lloyd Trefethen, Barycentric Lagrange Interpolation, SIAM Review, Volume 46, Number 3, September 2004, pages 501-517. Parameters: Input, int ND, the number of data points. ND must be at least 1. Input, double XD[ND], the data points. Input, double YD[ND], the data values. Input, int NI, the number of interpolation points. Input, double XI[NI], the interpolation points. Output, double LAGCHEBY2_INTERP_1D[NI], the interpolated values. */ { double *denom; int *exact; int i; int j; double *numer; double t; double wd; double *yi; denom = ( double * ) malloc ( ni * sizeof ( double ) ); exact = ( int * ) malloc ( ni * sizeof ( int ) ); numer = ( double * ) malloc ( ni * sizeof ( double ) ); yi = ( double * ) malloc ( ni * sizeof ( double ) ); for ( i = 0; i < ni; i++ ) { exact[i] = -1; denom[i] = 0.0; numer[i] = 0.0; } for ( j = 0; j < nd; j++ ) { wd = r8_mop ( j ); if ( j == 0 || j == nd - 1 ) { wd = 0.5 * wd; } for ( i = 0; i < ni; i++ ) { if ( xi[i] == xd[j] ) { exact[i] = j; numer[i] = yd[j]; denom[i] = 1.0; } if ( exact[i] == -1 ) { t = wd / ( xi[i] - xd[j] ); numer[i] = numer[i] + t * yd[j]; denom[i] = denom[i] + t; } } } for ( i = 0; i < ni; i++ ) { yi[i] = numer[i] / denom[i]; } return yi; } /******************************************************************************/ double *lageven_interp_1d ( int nd, double xd[], double yd[], int ni, double xi[] ) /******************************************************************************/ /* Purpose: LAGEVEN_VALUE_1D evaluates the Lagrange evenly-spaced interpolant. Discussion: The weight vector WD computed below is only valid if the data points XD are, as expected, evenly spaced in an interval [A,B] with spacing (B-A)/N. The XD values might be computed by: xd[i] = ( ( 2 * nd - 2 * i + 1 ) * a + ( 2 * i - 1 ) * b ) / ( 2 * nd ) for instance. Thanks to John Ferrier for pointing out that DENOM and NUMER needed to be initialized, 16 September 2013. Licensing: This code is distributed under the MIT license. Modified: 16 September 2013 Author: John Burkardt Reference: Jean-Paul Berrut, Lloyd Trefethen, Barycentric Lagrange Interpolation, SIAM Review, Volume 46, Number 3, September 2004, pages 501-517. Parameters: Input, int ND, the number of data points. ND must be at least 1. Input, double XD[ND], the data points. Input, double YD[ND], the data values. Input, int NI, the number of interpolation points. Input, double XI[NI], the interpolation points. Output, double LAGEVEN_INTERP_1D[NI], the interpolated values. */ { double *denom; int *exact; int i; int j; double *numer; double t; double wd; double *yi; denom = ( double * ) malloc ( ni * sizeof ( double ) ); exact = ( int * ) malloc ( ni * sizeof ( int ) ); numer = ( double * ) malloc ( ni * sizeof ( double ) ); yi = ( double * ) malloc ( ni * sizeof ( double ) ); for ( i = 0; i < ni; i++ ) { exact[i] = -1; denom[i] = 0.0; numer[i] = 0.0; } for ( j = 0; j < nd; j++ ) { wd = r8_mop ( j ) * r8_choose ( nd, j ); for ( i = 0; i < ni; i++ ) { if ( xi[i] == xd[j] ) { exact[i] = j; numer[i] = yd[j]; denom[i] = 1.0; } if ( exact[i] == -1 ) { t = wd / ( xi[i] - xd[j] ); numer[i] = numer[i] + t * yd[j]; denom[i] = denom[i] + t; } } } for ( i = 0; i < ni; i++ ) { yi[i] = numer[i] / denom[i]; } return yi; }