35 #include <boost/function.hpp>
36 #include <boost/numeric/ublas/vector.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 #include <boost/numeric/ublas/io.hpp>
45 #include <feel/feelcore/traits.hpp>
49 namespace ublas=boost::numeric::ublas;
69 template<
int N,
typename T =
double>
78 static const int order = N;
79 static const int nOrder = N;
101 Jacobi( value_type a = value_type( 0.0 ), value_type b = value_type( 0.0 ) )
122 self_type& operator=( self_type
const& p )
144 value_type
operator()( value_type
const& x )
const;
152 uint16_type degree()
const
181 value_type
value( value_type
const& x )
const
186 value_type derivate( value_type
const& x )
const;
198 template<
int N,
typename T>
199 typename Jacobi<N, T>::value_type
202 const value_type one = 1.0;
203 const value_type two = 2.0;
209 return 0.5 * ( M_a - M_b + ( M_a + M_b + two ) * x );
213 value_type apb = M_a + M_b;
214 value_type pn2 = one;
215 value_type pn1 = 0.5 * ( M_a - M_b + ( apb + two ) * x );
218 for (
int k = 2; k < N+1; ++k )
220 value_type kv = value_type( k );
221 value_type a1 = two * kv * ( kv + apb ) * ( two * kv + apb - two );
222 value_type a2 = ( two * kv + apb - one ) * ( M_a * M_a - M_b * M_b );
223 value_type a3 = ( ( two * kv + apb - two )
224 * ( two * kv + apb - one )
225 * ( two * kv + apb ) );
226 value_type a4 = ( two * ( kv + M_a - one ) * ( kv + M_b - one )
227 * ( two * kv + apb ) );
229 p = ( ( a2 + a3 * x ) * pn1 - a4 * pn2 )/a1;
237 template<
int N,
typename T>
238 typename Jacobi<N, T>::value_type
244 Jacobi<N-1, T> dp( M_a + 1.0, M_b + 1.0 );
245 value_type Nv = value_type( N );
246 return 0.5 * ( M_a + M_b + Nv + 1.0 ) * dp( x );
269 template<
typename T =
double>
284 typedef T value_type;
298 Jacobi( uint16_type N, value_type a = value_type( 0.0 ), value_type b = value_type( 0.0 ) )
307 M_degree( p.M_degree ),
321 self_type& operator=( self_type
const& p )
325 M_degree = p.M_degree;
344 value_type
operator()( value_type
const& x )
const;
352 uint16_type degree()
const
363 void setDegree( uint16_type N )
385 value_type
value( value_type
const& x )
const
390 value_type derivate( value_type
const& x )
const;
399 uint16_type M_degree;
404 typename Jacobi<T>::value_type
407 const uint16_type N = this->M_degree;
408 const value_type one = 1.0;
409 const value_type two = 2.0;
415 return 0.5 * ( M_a - M_b + ( M_a + M_b + two ) * x );
419 value_type apb = M_a + M_b;
420 value_type pn2 = one;
421 value_type pn1 = 0.5 * ( M_a - M_b + ( apb + two ) * x );
424 for ( uint16_type k = 2; k < N+1; ++k )
426 value_type kv = value_type( k );
427 value_type a1 = two * kv * ( kv + apb ) * ( two * kv + apb - two );
428 value_type a2 = ( two * kv + apb - one ) * ( M_a * M_a - M_b * M_b );
429 value_type a3 = ( ( two * kv + apb - two )
430 * ( two * kv + apb - one )
431 * ( two * kv + apb ) );
432 value_type a4 = ( two * ( kv + M_a - one ) * ( kv + M_b - one )
433 * ( two * kv + apb ) );
435 p = ( ( a2 + a3 * x ) * pn1 - a4 * pn2 )/a1;
444 typename Jacobi<T>::value_type
447 const uint16_type N = this->M_degree;
452 Jacobi<T> dp( N-1, M_a + 1.0, M_b + 1.0 );
453 value_type Nv = value_type( N );
454 return 0.5 * ( M_a + M_b + Nv + 1.0 ) * dp( x );
459 template<u
int16_type N,
typename T>
460 typename Eigen::Matrix<T,N+1,Eigen::Dynamic>
461 JacobiBatchEvaluation( T a, T b, Matrix<T,Eigen::Dynamic,1>
const& __pts )
463 typedef T value_type;
464 typedef Eigen::Matrix<T,Eigen::Dynamic,N+1> matrix_type;
465 typedef Eigen::Matrix<T,Eigen::Dynamic,1> vector_type;
466 matrix_type res( __pts.size(), N+1 );
467 res.col( 0 ).setOnes();
472 res.col( 1 ) = 0.5*( res.col( 1 ).Constant( res.rows(),a-b )+( a+b+2 )*__pts );
473 value_type apb = a + b;
475 for (
int k = 2; k < N+1; ++k )
477 value_type kv = value_type( k );
478 value_type a1 = 2.0 * kv * ( kv + apb ) * ( 2.0 * kv + apb - 2.0 );
479 value_type a2 = ( 2.0 * kv + apb - 1.0 ) * ( a * a - b * b );
480 value_type a3 = ( 2.0 * kv + apb - 2.0 ) * ( 2.0 * kv + apb - 1.0 ) * ( 2.0 * kv + apb );
481 value_type a4 = 2.0 * ( kv + a - 1.0 ) * ( kv + b - 1.0 ) * ( 2.0 * kv + apb );
487 ( a2* res.col( k-1 ) +
488 a3 * ( __pts.cwise()*res.col( k-1 ) ) -
489 a4 * res.col( k-2 ) );
493 return res.transpose();
495 template<
typename T,u
int16_type N, u
int16_type Np>
496 typename Eigen::Matrix<T,N+1,Np>
497 JacobiBatchEvaluation( T a, T b, Matrix<T,Np,1>
const& __pts )
499 typedef T value_type;
500 typedef Eigen::Matrix<T,Np,N+1> matrix_type;
501 typedef Eigen::Matrix<T,Np,1> vector_type;
502 matrix_type res( Np, N+1 );
503 res.col( 0 ).setOnes();
508 res.col( 1 ) = 0.5*( res.col( 1 ).Constant( res.rows(),a-b )+( a+b+2 )*__pts );
509 value_type apb = a + b;
511 for (
int k = 2; k < N+1; ++k )
513 value_type kv = value_type( k );
514 value_type a1 = 2.0 * kv * ( kv + apb ) * ( 2.0 * kv + apb - 2.0 );
515 value_type a2 = ( 2.0 * kv + apb - 1.0 ) * ( a * a - b * b );
516 value_type a3 = ( 2.0 * kv + apb - 2.0 ) * ( 2.0 * kv + apb - 1.0 ) * ( 2.0 * kv + apb );
517 value_type a4 = 2.0 * ( kv + a - 1.0 ) * ( kv + b - 1.0 ) * ( 2.0 * kv + apb );
523 ( a2* res.col( k-1 ) +
524 a3 * ( __pts.cwise()*res.col( k-1 ) ) -
525 a4 * res.col( k-2 ) );
529 return res.transpose();
533 template<u
int16_type N,
typename T>
535 JacobiBatchEvaluation( T a, T b, ublas::vector<T>
const& __pts )
537 typedef T value_type;
538 ublas::matrix<T> res( N+1, __pts.size() );
539 ublas::row( res, 0 ) = ublas::scalar_vector<value_type>( res.size2(), 1.0 );
543 ublas::row( res, 1 ) = 0.5 * ( ublas::scalar_vector<value_type>( res.size2(), a - b ) + ( a + b + 2.0 ) * __pts );
544 value_type apb = a + b;
546 for (
int k = 2; k < N+1; ++k )
548 value_type kv = value_type( k );
549 value_type a1 = 2.0 * kv * ( kv + apb ) * ( 2.0 * kv + apb - 2.0 );
550 value_type a2 = ( 2.0 * kv + apb - 1.0 ) * ( a * a - b * b );
551 value_type a3 = ( 2.0 * kv + apb - 2.0 ) * ( 2.0 * kv + apb - 1.0 ) * ( 2.0 * kv + apb );
552 value_type a4 = 2.0 * ( kv + a - 1.0 ) * ( kv + b - 1.0 ) * ( 2.0 * kv + apb );
557 ublas::row( res, k ) = a2* ublas::row( res, k-1 ) +
558 a3 * ublas::element_prod( __pts, ublas::row( res, k-1 ) ) -
559 a4 * ublas::row( res, k-2 );
566 template<u
int16_type N,
typename T>
568 JacobiBatchDerivation( ublas::matrix<T>& res, T a, T b, ublas::vector<T>
const& __pts, mpl::bool_<true> )
570 typedef T value_type;
571 ublas::subrange( res, 1, N+1, 0, __pts.size() ) = JacobiBatchEvaluation<N-1, T>( a+1.0, b+1.0, __pts );
573 for ( uint16_type i = 1; i < N+1; ++i )
574 ublas::row( res, i ) *= 0.5*( a+b+value_type( i )+1.0 );
577 template<u
int16_type N,
typename T>
579 JacobiBatchDerivation( ublas::matrix<T>& , T , T ,
580 ublas::vector<T>
const& , mpl::bool_<false> )
584 template<u
int16_type N,
typename T>
586 JacobiBatchDerivation( T a, T b, ublas::vector<T>
const& __pts )
588 typedef T value_type;
589 ublas::matrix<T> res( N+1, __pts.size() );
590 ublas::row( res, 0 ) = ublas::scalar_vector<value_type>( res.size2(), 0.0 );
591 static const bool cond = N>0;
592 JacobiBatchDerivation<N,T>( res, a, b, __pts, mpl::bool_<cond>() );
598 template<u
int16_type N,
typename T>
599 Eigen::Matrix<T,N+1,Eigen::Dynamic>
600 JacobiBatchDerivation( T a, T b, Eigen::Matrix<T,Eigen::Dynamic,1>
const& __pts )
602 typedef T value_type;
603 Eigen::Matrix<T,N+1,Eigen::Dynamic> res( N+1, __pts.size() );
604 res.row( 0 ).setZero();
608 res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<N-1, T>( a+1.0, b+1.0, __pts );
611 for ( uint16_type i = 1; i < N+1; ++i )
612 res.row( i ) *= 0.5*( a+b+value_type( i )+1.0 );
620 Eigen::Matrix<T,N+1,Ncols>
621 JacobiBatchDerivation( T a, T b, Eigen::Matrix<T,Ncols,1>
const& __pts )
623 typedef T value_type;
624 Eigen::Matrix<T,N+1,Ncols> res;
625 res.row( 0 ).setZero();
629 res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<T,N-1,Ncols>( a+1.0, b+1.0, __pts );
632 for ( uint16_type i = 1; i < N+1; ++i )
633 res.row( i ) *= 0.5*( a+b+value_type( i )+1.0 );
639 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
640 JacobiBatchDerivation( uint16_type N,
641 T a, T b, Eigen::Matrix<T,Eigen::Dynamic,1>
const& __pts )
643 typedef T value_type;
644 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> res( N+1, __pts.size() );
645 res.row( 0 ).setZero();
649 res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<T>( N-1,
654 for ( uint16_type i = 1; i < N+1; ++i )
655 res.row( i ) *= 0.5*( a+b+value_type( i )+1.0 );
666 JacobiBatchEvaluation( uint16_type N, T a, T b, ublas::vector<T>
const& __pts )
668 typedef T value_type;
669 ublas::matrix<T> res( N+1, __pts.size() );
670 ublas::row( res, 0 ) = ublas::scalar_vector<value_type>( res.size2(), 1.0 );
674 ublas::row( res, 1 ) = 0.5 * ( ublas::scalar_vector<value_type>( res.size2(), a - b ) + ( a + b + 2.0 ) * __pts );
675 value_type apb = a + b;
677 for ( uint16_type k = 2; k < N+1; ++k )
679 value_type kv = value_type( k );
680 value_type a1 = 2.0 * kv * ( kv + apb ) * ( 2.0 * kv + apb - 2.0 );
681 value_type a2 = ( 2.0 * kv + apb - 1.0 ) * ( a * a - b * b );
682 value_type a3 = ( 2.0 * kv + apb - 2.0 ) * ( 2.0 * kv + apb - 1.0 ) * ( 2.0 * kv + apb );
683 value_type a4 = 2.0 * ( kv + a - 1.0 ) * ( kv + b - 1.0 ) * ( 2.0 * kv + apb );
688 ublas::row( res, k ) = a2* ublas::row( res, k-1 ) +
689 a3 * ublas::element_prod( __pts, ublas::row( res, k-1 ) ) -
690 a4 * ublas::row( res, k-2 );
699 JacobiBatchDerivation( uint16_type N, T a, T b, ublas::vector<T>
const& __pts )
701 typedef T value_type;
702 ublas::matrix<T> res( N+1, __pts.size() );
703 ublas::row( res, 0 ) = ublas::scalar_vector<value_type>( res.size2(), 0.0 );
707 ublas::subrange( res, 1, N+1, 0, __pts.size() ) = JacobiBatchEvaluation<T>( N-1, a+1.0, b+1.0, __pts );
709 for ( uint16_type i = 1; i < N+1; ++i )
710 ublas::row( res, i ) *= 0.5*( a+b+value_type( i )+1.0 );
726 template<
typename JacobiP,
typename Vector>
728 roots( JacobiP
const& p, Vector& xr )
730 const uint16_type N = p.degree();
731 typedef typename JacobiP::value_type value_type;
735 value_type eps = type_traits<value_type>::epsilon();
739 for (
int k = 0; k < N; ++k )
741 value_type pi = 4.0*math::atan( value_type( 1.0 ) );
743 r = -math::cos( ( 2.0*value_type( k ) + 1.0 ) * pi / ( 2.0 * value_type( N ) ) );
747 r = 0.5 * ( r + xr[k-1] );
750 value_type jf = 2.0*eps;
751 value_type delta = 2.0*eps;
758 for (
int i = 0; i < k; ++i )
760 s += value_type( 1.0 ) / ( r - xr[i] );
764 value_type jdf = p.derivate( r );
765 delta = jf / ( jdf - jf * s );
772 while ( math::abs( jf ) > eps && j < max_iter );
786 return n*fact( n-1.0 );
788 template<
int N,
typename T,
typename VectorW,
typename VectorN>
790 gaussjacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
792 typedef T value_type;
794 Jacobi<N, T> p( a, b );
798 const value_type two = 2.0;
799 const value_type power = a+b+1.0;
800 value_type a1 = math::pow( two,power );
801 value_type a2 = fact( a+value_type( N ) );
802 value_type a3 = fact( b+value_type( N ) );
803 value_type a4 = fact( a+b+value_type( N ) );
804 value_type a5 = fact( value_type( N ) );
805 value_type a6 = a1 * a2 * a3;
807 for (
int k = 0; k < N; ++k )
809 value_type fp = p.derivate( xr[k] );
810 value_type dn = fp*fp*( 1.0 - xr[k]*xr[k] )*a4*a5;
823 template<
int N,
typename T,
typename VectorW,
typename VectorN>
825 gausslobattojacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
827 typedef T value_type;
829 Jacobi<N-2, T> p( a+ 1.0, b+ 1.0 );
831 Jacobi<N-1, T> q( a, b );
833 VectorN prexr( N-2 );
839 for (
int i = 1; i < N-1; ++i )
844 const value_type two = 2.0;
846 value_type a1 = math::pow( two,
int( a+b+1.0 ) );
847 value_type a2 = fact( a+value_type( N ) -1.0 );
848 value_type a3 = fact( b+value_type( N ) -1.0 );
850 value_type a4 = fact( a+b+value_type( N ) );
851 value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )-1.0 );
853 value_type a6 = a1 * a2 * a3;
855 for (
int k = 0; k < N; ++k )
857 value_type fq = q.value( xr[k] );
859 value_type dn = fq*fq*a4*a5;
864 wr[0] = wr[0]*( b+1.0 );
865 wr[N-1] = wr[N-1]*( a+1.0 );
875 template<
typename T,
typename VectorW,
typename VectorN>
877 gausslobattojacobi(
size_type N, VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ),
bool interior =
false )
879 wr.resize( N - 2*interior );
880 xr.resize( N - 2*interior );
882 typedef T value_type;
888 VectorN prexr( N-2 );
898 for (
size_type i = 1-interior; i < N-( 1+interior ); ++i )
899 xr[i] = prexr[i-( 1-interior )];
902 const value_type two = 2.0;
904 value_type a1 = math::pow( two,
int( a+b+1.0 ) );
905 value_type a2 = fact( a+value_type( N ) -1.0 );
906 value_type a3 = fact( b+value_type( N ) -1.0 );
908 value_type a4 = fact( a+b+value_type( N ) );
909 value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )-1.0 );
911 value_type a6 = a1 * a2 * a3;
913 for (
int k = 1-interior; k < int( N-1-interior ); ++k )
915 value_type fq = q.value( xr[k] );
917 value_type dn = fq*fq*a4*a5;
924 wr[0] = wr[0]*( b+1.0 );
925 wr[N-1] = wr[N-1]*( a+1.0 );
937 template<
int N,
typename T,
typename VectorW,
typename VectorN>
939 left_gaussradaujacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
941 typedef T value_type;
943 Jacobi<N-1, T> p( a, b+1.0 );
945 Jacobi<N-1, T> q( a, b );
947 VectorN prexr( N-1 );
953 for (
int i = 1; i < N; ++i )
956 const value_type two = 2.0;
958 value_type a1 = math::pow( two,
int( a+b ) );
959 value_type a2 = fact( a+value_type( N ) -1.0 );
960 value_type a3 = fact( b+value_type( N ) -1.0 );
962 value_type a4 = fact( a+b+value_type( N ) );
963 value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )+b );
965 value_type a6 = a1 * a2 * a3;
967 for (
int k = 0; k < N; ++k )
969 value_type fq = q.value( xr[k] );
971 value_type dn = fq*fq*a4*a5;
973 wr[k] =a6*( 1.0-xr[k] ) /dn ;
976 wr[0] = wr[0]*( b+1.0 );
980 template<
int N,
typename T,
typename VectorW,
typename VectorN>
982 right_gaussradaujacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
984 typedef T value_type;
986 Jacobi<N-1, T> p( a+1.0, b );
988 Jacobi<N-1, T> q( a, b );
990 VectorN prexr( N-1 );
994 for (
int i = 0; i < N-1; ++i )
999 const value_type two = 2.0;
1001 value_type a1 = math::pow( two,
int( a+b ) );
1002 value_type a2 = fact( a+value_type( N ) -1.0 );
1003 value_type a3 = fact( b+value_type( N ) -1.0 );
1005 value_type a4 = fact( a+b+value_type( N ) );
1006 value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )+a );
1008 value_type a6 = a1 * a2 * a3;
1010 for (
int k = 0; k < N; ++k )
1012 value_type fq = q.value( xr[k] );
1014 value_type dn = fq*fq*a4*a5;
1016 wr[k] =a6*( 1.0+xr[k] ) /dn ;
1019 wr[N-1] = wr[N-1]*( a+1.0 );
1023 template<
int N,
typename T>
1024 T
integrate( boost::function<T( T
const& )>
const& f )
1026 typedef T value_type;
1027 ublas::vector<T> xr( Jacobi<N, T>::nOrder );
1028 ublas::vector<T> wr( Jacobi<N, T>::nOrder );
1031 details::gaussjacobi<N, T, ublas::vector<T> >( wr, xr );
1033 value_type res = 0.0;
1035 for (
int k = 0; k < Jacobi<N, T>::nOrder; ++k )
1036 res += wr[k]*f( xr[k] );