#! /usr/bin/env python # def w_mass_matrix ( n ): #*****************************************************************************80 # ## W_MASS_MATRIX computes the mass matrix for the Chebyshev W polynomial. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 20 July 2015 # # Author: # # John Burkardt # import numpy as np from w_polynomial import w_polynomial from w_quadrature_rule import w_quadrature_rule x, w = w_quadrature_rule ( n + 1 ) phi = w_polynomial ( n + 1, n, x ) phiw = np.zeros ( [ n + 1, n + 1 ] ) for i in range ( 0, n + 1 ): for j in range ( 0, n + 1 ): phiw[j,i] = w[i] * phi[i,j] a = np.dot ( phiw, phi ) return a def w_mass_matrix_test ( ): #*****************************************************************************80 # ## W_MASS_MATRIX_TEST tests W_MASS_MATRIX. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 20 July 2015 # # Author: # # John Burkardt # import platform from r8mat_print import r8mat_print print ( '' ) print ( 'W_MASS_MATRIX_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' W_MASS_MATRIX computes the mass matrix for the' ) print ( ' Chebyshev polynomials W(i,x).' ) print ( ' A(I,J) = integral ( -1 <=x <= +1 ) W(i,x) W(j,x) sqrt(1-x)/sqrt(1+x) dx' ) print ( ' 0 if i is not equal to j' ) print ( ' pi if i = j.' ) n = 3 a = w_mass_matrix ( n ) r8mat_print ( n + 1, n + 1, a, ' W mass matrix:' ) # # Terminate. # print ( '' ) print ( 'W_MASS_MATRIX_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) w_mass_matrix_test ( ) timestamp ( )