#! /usr/bin/env python # def gegenbauer_polynomial_value ( m, n, alpha, x ): #*****************************************************************************80 # ## GEGENBAUER_POLYNOMIAL_VALUE computes the Gegenbauer polynomials C(I,ALPHA)(X). # # Differential equation: # # (1-X*X) Y'' - (2 ALPHA + 1) X Y' + M (M + 2 ALPHA) Y = 0 # # Recursion: # # C(0,ALPHA,X) = 1, # C(1,ALPHA,X) = 2*ALPHA*X # C(M,ALPHA,X) = ( ( 2*M-2+2*ALPHA) * X * C(M-1,ALPHA,X) # + ( -M+2-2*ALPHA) * C(M-2,ALPHA,X) ) / M # # Restrictions: # # ALPHA must be greater than -0.5. # # Special values: # # If ALPHA = 1, the Gegenbauer polynomials reduce to the Chebyshev # polynomials of the second kind. # # Norm: # # Integral ( -1 <= X <= 1 ) ( 1 - X^2 )^( ALPHA - 0.5 ) * C(M,ALPHA,X)^2 dX # # = PI * 2^( 1 - 2 * ALPHA ) * Gamma ( M + 2 * ALPHA ) # / ( M! * ( M + ALPHA ) * ( Gamma ( ALPHA ) )^2 ) # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 23 November 2015 # # Author: # # John Burkardt # # Reference: # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Parameters: # # Input, integer M, the highest order polynomial to compute. # Note that polynomials 0 through N will be computed. # # Input, integer N, the number of evaluation points. # # Input, real ALPHA, a parameter which is part of the definition of # the Gegenbauer polynomials. It must be greater than -0.5. # # Input, real X(N), the evaluation points. # # Output, real C(1:M+1,N), the values of Gegenbauer polynomials 0 through M # at the N points X. # import numpy as np from gegenbauer_alpha_check import gegenbauer_alpha_check from sys import exit check = gegenbauer_alpha_check ( alpha ) if ( not check ): print ( '' ) print ( 'GEGENBAUER_POLYNOMIAL_VALUE - Fatal error!' ) print ( ' GEGENBAUER_ALPHA_CHECK failed.' ) exit ( 'GEGENBAUER_POLYNOMIAL_VALUE - Fatal error!' ) c = np.zeros ( [ m + 1, n ] ) if ( m < 0 ): return c for j in range ( 0, n ): c[0,j] = 1.0 if ( n == 0 ): return c if ( m < 1 ): return c for j in range ( 0, n ): c[1,j] = 2.0 * alpha * x[j] for i in range ( 2, m + 1 ): for j in range ( 0, n ): c[i,j] = ( ( 2 * i - 2 + 2.0 * alpha ) * x[j] * c[i-1,j] \ + ( - i + 2 - 2.0 * alpha ) * c[i-2,j] ) \ / i return c def gegenbauer_polynomial_value_test ( ): #*****************************************************************************80 # ## GEGENBAUER_POLYNOMIAL_VALUE_TEST tests GEGENBAUER_POLYNOMIAL_VALUE. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 23 November 2015 # # Author: # # John Burkardt # import numpy as np import platform from gegenbauer_polynomial_values import gegenbauer_polynomial_values print ( '' ) print ( 'GEGENBAUER_POLYNOMIAL_VALUE_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' GEGENBAUER_POLYNOMIAL_VALUE evaluates the Gegenbauer polynomial.' ) print ( '' ) print ( ' M ALPHA X GPV GEGENBAUER' ) print ( '' ) n_data = 0 x = np.zeros ( 1 ) while ( True ): n_data, m, alpha, xscalar, fx = gegenbauer_polynomial_values ( n_data ) if ( n_data == 0 ): break # # Since GEGENBAUER_POLYNOMIAL_VALUE expects a vector input X, we have to # do a little "rewrapping" of the input. # n = 1 x[0] = xscalar c = gegenbauer_polynomial_value ( m, n, alpha, x ) print ( ' %6d %8.2f %8.2f %12f %12f' % ( m, alpha, x[0], fx, c[m,0] ) ) # # Terminate. # print ( '' ) print ( 'GEGENBAUER_POLYNOMIAL_VALUE_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) gegenbauer_polynomial_value_test ( ) timestamp ( )