program main !*****************************************************************************80 ! !! JACOBI_POLYNOMIAL_TEST() tests JACOBI_POLYNOMIAL(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 April 2012 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'JACOBI_POLYNOMIAL_TEST' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test the JACOBI_POLYNOMIAL library.' call test01 ( ) call test02 ( ) call test03 ( ) call test04 ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'JACOBI_POLYNOMIAL_TEST' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine test01 ( ) !*****************************************************************************80 ! !! TEST01 tests J_POLYNOMIAL. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 April 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) e real ( kind = rk ) fx1 real ( kind = rk ) fx2 real ( kind = rk ), allocatable :: fx2_vec(:,:) integer m integer n integer n_data real ( kind = rk ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01:' write ( *, '(a)' ) ' J_POLYNOMIAL_VALUES stores values of' write ( *, '(a)' ) ' the Jacobi polynomials.' write ( *, '(a)' ) ' J_POLYNOMIAL evaluates the polynomial.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Tabulated Computed' write ( *, '(a)' ) ' N A B X J(N,A,B,X) J(N,A,B,X) Error' write ( *, '(a)' ) ' ' n_data = 0 do call j_polynomial_values ( n_data, n, a, b, x, fx1 ) if ( n_data == 0 ) then exit end if allocate ( fx2_vec(1,0:n) ) m = 1 call j_polynomial ( m, n, a, b, x, fx2_vec ) fx2 = fx2_vec(1,n) e = fx1 - fx2 write ( *, '(2x,i4,2x,f6.2,2x,f6.2,2x,f6.2,2x,g24.16,2x,g24.16,2x,g8.2)' ) & n, a, b, x, fx1, fx2, e deallocate ( fx2_vec ) end do return end subroutine test02 ( ) !*****************************************************************************80 ! !! TEST02 tests J_POLYNOMIAL_ZEROS. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 April 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: test_num = 3 real ( kind = rk ) a real ( kind = rk ), dimension ( test_num) :: a_test = (/ & 0.5D+00, 1.0D+00, 2.0D+00 /) real ( kind = rk ) b real ( kind = rk ), dimension ( test_num ) :: b_test = (/ & 0.5D+00, 1.5D+00, 0.5D+00 /) integer degree real ( kind = rk ), allocatable :: hz(:,:) integer test character ( len = 80 ) title real ( kind = rk ), allocatable :: z(:) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02:' write ( *, '(a)' ) ' J_POLYNOMIAL_ZEROS computes the zeros of J(n,a,b,x);' write ( *, '(a)' ) ' Check by calling J_POLYNOMIAL there.' do test = 1, 3 a = a_test(test); b = b_test(test); do degree = 1, 5 allocate ( z(1:degree) ) call j_polynomial_zeros ( degree, a, b, z ) write ( title, '(a,i1,a,f3.1,a,f3.1,a)' ) 'Zeros for J(', degree, ',', a, ',', b, ')' call r8vec_print ( degree, z, title ) allocate ( hz(degree,degree+1) ) call j_polynomial ( degree, degree, a, b, z, hz ) write ( title, '(a,i1,a,f3.1,a,f3.1,a)' ) 'Evaluate J(', degree, ',', a, ',', b, ')' call r8vec_print ( degree, hz(1:degree,degree+1), title ) deallocate ( hz ) deallocate ( z ) end do end do return end subroutine test03 ( ) !*****************************************************************************80 ! !! TEST03 tests J_QUADRATURE_RULE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 April 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) b integer i integer j real ( kind = rk ) j_double_product_integral real ( kind = rk ), allocatable :: ji(:,:) real ( kind = rk ), allocatable :: jj(:,:) integer k integer n real ( kind = rk ) q real ( kind = rk ) q_exact real ( kind = rk ), allocatable :: w(:) real ( kind = rk ), allocatable :: x(:) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03:' write ( *, '(a)' ) ' J_QUADRATURE_RULE computes the quadrature rule' write ( *, '(a)' ) ' associated with J(n,a,b,x);' n = 7 a = 1.0D+00 b = 2.5D+00 allocate ( x(1:n) ) allocate ( w(1:n) ) call j_quadrature_rule ( n, a, b, x, w ) call r8vec2_print ( n, x, w, ' X W' ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Use the quadrature rule to estimate:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Q = Integral (-1