C 代码 使用除差计算多项式插值 到给定的数据集.rar

# include <stdlib.h> # include <stdio.h> # include <math.h> # include <time.h> # include <string.h> # include "divdif.h" /******************************************************************************/ double *cheby_t_zero ( int n ) /******************************************************************************/ /* Purpose: CHEBY_T_ZERO returns zeroes of the Chebyshev polynomial T(N)(X). Discussion: The I-th zero of T(N)(X) is cos((2*I-1)*PI/(2*N)), I = 1 to N Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Output, double CHEBY_T_ZERO[N], the zeroes of T(N)(X). */ { double angle; int i; static double pi = 3.141592653589793; double *z; z = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { angle = ( double ) ( 2 * i + 1 ) * pi / ( double ) ( 2 * n ); z[i] = cos ( angle ); } return z; } /******************************************************************************/ double *cheby_u_zero ( int n ) /******************************************************************************/ /* Purpose: CHEBY_U_ZERO returns zeroes of the Chebyshev polynomial U(N)(X). Discussion: The I-th zero of U(N)(X) is cos((I-1)*PI/(N-1)), I = 1 to N Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2011 Author: John Burkardt Parameters: Input, int N, the order of the polynomial. Output, double CHEBY_U_ZERO[N], the zeroes of U(N)(X). */ { double angle; int i; static double pi = 3.141592653589793; double *z; z = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { angle = ( double ) ( i + 1 ) * pi / ( double ) ( n + 1 ); z[i] = cos ( angle ); } return z; } /******************************************************************************/ void data_to_dif ( int ntab, double xtab[], double ytab[], double diftab[] ) /******************************************************************************/ /* Purpose: DATA_TO_DIF sets up a divided difference table from raw data. Discussion: Space can be saved by using a single array for both the DIFTAB and YTAB dummy parameters. In that case, the difference table will overwrite the Y data without interfering with the computation. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of pairs of points (XTAB[I],YTAB[I]) which are to be used as data. Input, double XTAB[NTAB], the X values at which data was taken. These values must be distinct. Input, double YTAB[NTAB], the corresponding Y values. Output, double DIFTAB[NTAB], the divided difference coefficients corresponding to the input (XTAB,YTAB). */ { int i; int j; /* Copy the data values into DIFTAB. */ for ( i = 0; i < ntab; i++ ) { diftab[i] = ytab[i]; } /* Make sure the abscissas are distinct. */ for ( i = 0; i < ntab; i++ ) { for ( j = i + 1; j < ntab; j++ ) { if ( xtab[i] - xtab[j] == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_DIF - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal!\n" ); fprintf ( stderr, " XTAB[%d] = %f\n", i, xtab[i] ); fprintf ( stderr, " XTAB[%d] = %f\n", j, xtab[j] ); exit ( 1 ); } } } /* Compute the divided differences. */ for ( i = 1; i <= ntab - 1; i++ ) { for ( j = ntab - 1; i <= j; j-- ) { diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] ); } } return; } /******************************************************************************/ void data_to_dif_display ( int ntab, double xtab[], double ytab[], double diftab[] ) /******************************************************************************/ /* Purpose: DATA_TO_DIF_DISPLAY sets up a divided difference table and prints out intermediate data. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Parameters: Input, int NTAB, the number of pairs of points (XTAB[I],YTAB[I]) which are to be used as data. Input, double XTAB[NTAB], the X values at which data was taken. These values must be distinct. Input, double YTAB[NTAB], the corresponding Y values. Output, double DIFTAB[NTAB], the divided difference coefficients corresponding to the input (XTAB,YTAB). */ { int i; int j; if ( !r8vec_distinct ( ntab, xtab ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_DIF_DISPLAY - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal!\n" ); exit ( 1 ); } printf ( "\n" ); printf ( " The divided difference table:\n" ); printf ( "\n" ); printf ( " " ); for ( i = 0; i < ntab; i++ ) { printf ( "%10f ", xtab[i] ); } printf ( "\n" ); printf ( "\n" ); printf ( "%6d ", 0 ); for ( i = 0; i < ntab; i++ ) { printf ( "%10f ", ytab[i] ); } printf ( "\n" ); /* Copy the data values into DIFTAB. */ for ( i = 0; i < ntab; i++ ) { diftab[i] = ytab[i]; } /* Compute the divided differences. */ for ( i = 1; i <= ntab-1; i++ ) { printf ( "%6d ", i ); for ( j = ntab-1; i <= j; j-- ) { diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] ); } for ( j = i; j < ntab; j++ ) { printf ( "%10f ", diftab[j] ); } printf ( "\n" ); } return; } /******************************************************************************/ void data_to_r8poly ( int ntab, double xtab[], double ytab[], double c[] ) /******************************************************************************/ /* Purpose: DATA_TO_R8POLY computes the coefficients of a polynomial interpolating data. Discussion: Space can be saved by using a single array for both the C and YTAB parameters. In that case, the coefficients will overwrite the Y data without interfering with the computation. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 February 2011 Author: John Burkardt Reference: Carl deBoor, A Practical Guide to Splines, Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. Parameters: Input, int NTAB, the number of data points. Input, double XTAB[NTAB], YTAB[NTAB], the data values. Output, double C[NTAB], the coefficients of the polynomial that passes through the data (XTAB,YTAB). C(0) is the constant term. */ { if ( !r8vec_distinct ( ntab, xtab ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DATA_TO_R8POLY - Fatal error!\n" ); fprintf ( stderr, " Two entries of XTAB are equal.\n" ); exit ( 1 ); } data_to_dif ( ntab, xtab, ytab, c ); dif_to_r8poly ( ntab, xtab, c, c ); return; } /******************************************************************************/ void dif_antideriv ( int ntab, double xtab[], double diftab[], int *ntab2, double xtab2[], double diftab2[] ) /******************************************************************************/ /* Purpose: DIF_ANTIDERIV integrates a polynomial in divided difference form. Discussion: This routine uses the divided difference representation (XTAB, DIFTAB) of a polynomial to compute the divided difference representation (XTAB, ANTTAB) of the antiderivative of the polynomial. The antiderivative of a polynomial P(X) is any polynomial Q(X) with the property that d/dX Q(X) = P(X). This routine chooses the antid