Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
jacobi.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-07-26
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2012 Universite Joseph Fourier Grenoble 1
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __Jacobi_H
31 #define __Jacobi_H 1
32 
33 #include <cmath>
34 
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>
40 
41 #if 0
42 #include <Eigen/Core>
43 #endif
44 
45 #include <feel/feelcore/traits.hpp>
46 
47 namespace Feel
48 {
49 namespace ublas=boost::numeric::ublas;
69 template<int N, typename T = double>
70 class Jacobi
71 {
72 public:
73 
77 
78  static const int order = N;
79  static const int nOrder = N;
80 
82 
86 
87  typedef T value_type;
88  typedef Jacobi<N, T> self_type;
89 
91 
95 
96 
101  Jacobi( value_type a = value_type( 0.0 ), value_type b = value_type( 0.0 ) )
102  :
103  M_a( a ),
104  M_b( b )
105  {}
106 
107  Jacobi( Jacobi const & p )
108  :
109  M_a( p.M_a ),
110  M_b( p.M_b )
111  {}
112 
113  ~Jacobi()
114  {}
115 
117 
121 
122  self_type& operator=( self_type const& p )
123  {
124  if ( this != &p )
125  {
126  M_a = p.M_a;
127  M_b = p.M_b;
128  }
129 
130  return *this;
131  }
132 
144  value_type operator()( value_type const& x ) const;
145 
147 
151 
152  uint16_type degree() const
153  {
154  return N;
155  }
156 
158 
162 
163 
165 
169 
181  value_type value( value_type const& x ) const
182  {
183  return this->operator()( x );
184  }
185 
186  value_type derivate( value_type const& x ) const;
187 
189 
190 
191 
192 protected:
193 
194 private:
195  value_type M_a;
196  value_type M_b;
197 };
198 template<int N, typename T>
199 typename Jacobi<N, T>::value_type
200 Jacobi<N, T>::operator()( value_type const& x ) const
201 {
202  const value_type one = 1.0;
203  const value_type two = 2.0;
204 
205  if ( N == 0 )
206  return one;
207 
208  else if ( N == 1 )
209  return 0.5 * ( M_a - M_b + ( M_a + M_b + two ) * x );
210 
211  else // N >= 2
212  {
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 );
216  value_type p = 0.0;
217 
218  for ( int k = 2; k < N+1; ++k )
219  {
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 ) );
228 
229  p = ( ( a2 + a3 * x ) * pn1 - a4 * pn2 )/a1;
230  pn2 = pn1;
231  pn1 = p;
232  }
233 
234  return p;
235  }
236 }
237 template<int N, typename T>
238 typename Jacobi<N, T>::value_type
239 Jacobi<N, T>::derivate( value_type const& x ) const
240 {
241  if ( N == 0 )
242  return 0.0;
243 
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 );
247 }
248 namespace dyna
249 {
269 template<typename T = double>
270 class Jacobi
271 {
272 public:
273 
277 
279 
283 
284  typedef T value_type;
285  typedef Jacobi<T> self_type;
286 
288 
292 
293 
298  Jacobi( uint16_type N, value_type a = value_type( 0.0 ), value_type b = value_type( 0.0 ) )
299  :
300  M_degree( N ),
301  M_a( a ),
302  M_b( b )
303  {}
304 
305  Jacobi( Jacobi const & p )
306  :
307  M_degree( p.M_degree ),
308  M_a( p.M_a ),
309  M_b( p.M_b )
310  {}
311 
312  ~Jacobi()
313  {}
314 
316 
320 
321  self_type& operator=( self_type const& p )
322  {
323  if ( this != &p )
324  {
325  M_degree = p.M_degree;
326  M_a = p.M_a;
327  M_b = p.M_b;
328  }
329 
330  return *this;
331  }
332 
344  value_type operator()( value_type const& x ) const;
345 
347 
351 
352  uint16_type degree() const
353  {
354  return M_degree;
355  }
356 
358 
362 
363  void setDegree( uint16_type N )
364  {
365  M_degree = N;
366  }
367 
369 
373 
385  value_type value( value_type const& x ) const
386  {
387  return this->operator()( x );
388  }
389 
390  value_type derivate( value_type const& x ) const;
391 
393 
394 
395 
396 protected:
397 
398 private:
399  uint16_type M_degree;
400  value_type M_a;
401  value_type M_b;
402 };
403 template<typename T>
404 typename Jacobi<T>::value_type
405 Jacobi<T>::operator()( value_type const& x ) const
406 {
407  const uint16_type N = this->M_degree;
408  const value_type one = 1.0;
409  const value_type two = 2.0;
410 
411  if ( N == 0 )
412  return one;
413 
414  else if ( N == 1 )
415  return 0.5 * ( M_a - M_b + ( M_a + M_b + two ) * x );
416 
417  else // N >= 2
418  {
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 );
422  value_type p = 0.0;
423 
424  for ( uint16_type k = 2; k < N+1; ++k )
425  {
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 ) );
434 
435  p = ( ( a2 + a3 * x ) * pn1 - a4 * pn2 )/a1;
436  pn2 = pn1;
437  pn1 = p;
438  }
439 
440  return p;
441  }
442 }
443 template<typename T>
444 typename Jacobi<T>::value_type
445 Jacobi<T>::derivate( value_type const& x ) const
446 {
447  const uint16_type N = this->M_degree;
448 
449  if ( N == 0 )
450  return 0.0;
451 
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 );
455 }
456 }
457 
458 #if 0
459 template<uint16_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 )
462 {
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();
468 
469  if ( N > 0 )
470  {
471  //ublas::col( res, 1 ) = 0.5 * ( ublas::scalar_vector<value_type>( res.size2(), a - b) + ( a + b + 2.0 ) * __pts );
472  res.col( 1 ) = 0.5*( res.col( 1 ).Constant( res.rows(),a-b )+( a+b+2 )*__pts );
473  value_type apb = a + b;
474 
475  for ( int k = 2; k < N+1; ++k )
476  {
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 );
482  a2 = a2 / a1;
483  a3 = a3 / a1;
484  a4 = a4 / a1;
485 
486  res.col( k ) =
487  ( a2* res.col( k-1 ) +
488  a3 * ( __pts.cwise()*res.col( k-1 ) ) -
489  a4 * res.col( k-2 ) );
490  }
491  }
492 
493  return res.transpose();
494 }
495 template<typename T,uint16_type N, uint16_type Np>
496 typename Eigen::Matrix<T,N+1,Np>
497 JacobiBatchEvaluation( T a, T b, Matrix<T,Np,1> const& __pts )
498 {
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();
504 
505  if ( N > 0 )
506  {
507  //ublas::col( res, 1 ) = 0.5 * ( ublas::scalar_vector<value_type>( res.size2(), a - b) + ( a + b + 2.0 ) * __pts );
508  res.col( 1 ) = 0.5*( res.col( 1 ).Constant( res.rows(),a-b )+( a+b+2 )*__pts );
509  value_type apb = a + b;
510 
511  for ( int k = 2; k < N+1; ++k )
512  {
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 );
518  a2 = a2 / a1;
519  a3 = a3 / a1;
520  a4 = a4 / a1;
521 
522  res.col( k ) =
523  ( a2* res.col( k-1 ) +
524  a3 * ( __pts.cwise()*res.col( k-1 ) ) -
525  a4 * res.col( k-2 ) );
526  }
527  }
528 
529  return res.transpose();
530 }
531 #endif // 0 - Eigen
532 
533 template<uint16_type N, typename T>
534 ublas::matrix<T>
535 JacobiBatchEvaluation( T a, T b, ublas::vector<T> const& __pts )
536 {
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 );
540 
541  if ( N > 0 )
542  {
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;
545 
546  for ( int k = 2; k < N+1; ++k )
547  {
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 );
553  a2 = a2 / a1;
554  a3 = a3 / a1;
555  a4 = a4 / a1;
556 
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 );
560  }
561  }
562 
563  return res;
564 }
565 
566 template<uint16_type N, typename T>
567 void
568 JacobiBatchDerivation( ublas::matrix<T>& res, T a, T b, ublas::vector<T> const& __pts, mpl::bool_<true> )
569 {
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 );
572 
573  for ( uint16_type i = 1; i < N+1; ++i )
574  ublas::row( res, i ) *= 0.5*( a+b+value_type( i )+1.0 );
575 }
576 
577 template<uint16_type N, typename T>
578 void
579 JacobiBatchDerivation( ublas::matrix<T>& /*res*/, T /*a*/, T /*b*/,
580  ublas::vector<T> const& /*__pts*/, mpl::bool_<false> )
581 {
582 }
583 
584 template<uint16_type N, typename T>
585 ublas::matrix<T>
586 JacobiBatchDerivation( T a, T b, ublas::vector<T> const& __pts )
587 {
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>() );
593  return res;
594 }
595 
596 
597 #if 0
598 template<uint16_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 )
601 {
602  typedef T value_type;
603  Eigen::Matrix<T,N+1,Eigen::Dynamic> res( N+1, __pts.size() );
604  res.row( 0 ).setZero();
605 
606  if ( N > 0 )
607  {
608  res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<N-1, T>( a+1.0, b+1.0, __pts );
609 
610  //ublas::subrange( res, 1, N+1, 0, __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 );
613  }
614 
615  return res;
616 }
617 template<typename T,
618  uint16_type N,
619  uint16_type Ncols>
620 Eigen::Matrix<T,N+1,Ncols>
621 JacobiBatchDerivation( T a, T b, Eigen::Matrix<T,Ncols,1> const& __pts )
622 {
623  typedef T value_type;
624  Eigen::Matrix<T,N+1,Ncols> res;
625  res.row( 0 ).setZero();
626 
627  if ( N > 0 )
628  {
629  res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<T,N-1,Ncols>( a+1.0, b+1.0, __pts );
630 
631  //ublas::subrange( res, 1, N+1, 0, __pts.size() ) = JacobiBatchEvaluation<N-1, T>( 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 );
634  }
635 
636  return res;
637 }
638 template<typename T>
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 )
642 {
643  typedef T value_type;
644  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> res( N+1, __pts.size() );
645  res.row( 0 ).setZero();
646 
647  if ( N > 0 )
648  {
649  res.block( 1, 0, N, __pts.size() ) = JacobiBatchEvaluation<T>( N-1,
650  a+1.0,
651  b+1.0,
652  __pts );
653 
654  for ( uint16_type i = 1; i < N+1; ++i )
655  res.row( i ) *= 0.5*( a+b+value_type( i )+1.0 );
656  }
657 
658  return res;
659 }
660 #endif // 0 - Eigen
661 
662 namespace dyna
663 {
664 template<typename T>
665 ublas::matrix<T>
666 JacobiBatchEvaluation( uint16_type N, T a, T b, ublas::vector<T> const& __pts )
667 {
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 );
671 
672  if ( N > 0 )
673  {
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;
676 
677  for ( uint16_type k = 2; k < N+1; ++k )
678  {
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 );
684  a2 = a2 / a1;
685  a3 = a3 / a1;
686  a4 = a4 / a1;
687 
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 );
691  }
692  }
693 
694  return res;
695 }
696 
697 template<typename T>
698 ublas::matrix<T>
699 JacobiBatchDerivation( uint16_type N, T a, T b, ublas::vector<T> const& __pts )
700 {
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 );
704 
705  if ( N > 0 )
706  {
707  ublas::subrange( res, 1, N+1, 0, __pts.size() ) = JacobiBatchEvaluation<T>( N-1, a+1.0, b+1.0, __pts );
708 
709  for ( uint16_type i = 1; i < N+1; ++i )
710  ublas::row( res, i ) *= 0.5*( a+b+value_type( i )+1.0 );
711  }
712 
713  return res;
714 }
715 } // dyna
716 
718 namespace details
719 {
720 
726 template<typename JacobiP, typename Vector>
727 void
728 roots( JacobiP const& p, Vector& xr )
729 {
730  const uint16_type N = p.degree();
731  typedef typename JacobiP::value_type value_type;
732 
733  if ( N != 0 )
734  {
735  value_type eps = type_traits<value_type>::epsilon();
736  value_type r;
737  int max_iter = 30;
738 
739  for ( int k = 0; k < N; ++k )
740  {
741  value_type pi = 4.0*math::atan( value_type( 1.0 ) );
742  // use k-th checbychev point to initiliaze newton
743  r = -math::cos( ( 2.0*value_type( k ) + 1.0 ) * pi / ( 2.0 * value_type( N ) ) );
744 
745  // use average of r and xr[k-1] as starting point (see KS)
746  if ( k > 0 )
747  r = 0.5 * ( r + xr[k-1] );
748 
749  int j = 0;
750  value_type jf = 2.0*eps;
751  value_type delta = 2.0*eps;
752 
753  do
754  {
755  // use deflation as proposed in KS
756  value_type s = 0.0;
757 
758  for ( int i = 0; i < k; ++i )
759  {
760  s += value_type( 1.0 ) / ( r - xr[i] );
761  }
762 
763  jf = p( r );
764  value_type jdf = p.derivate( r );
765  delta = jf / ( jdf - jf * s );
766 
767  // newton step done
768  r = r - delta;
769 
770  ++j;
771  }
772  while ( math::abs( jf ) > eps && j < max_iter );
773 
774  // store k-th root
775  xr[k] = r;
776  }
777  }
778 }
779 
780 template<typename T>
781 T fact( T n )
782 {
783  if ( n == 0.0 )
784  return 1;
785 
786  return n*fact( n-1.0 );
787 }
788 template<int N, typename T, typename VectorW, typename VectorN>
789 void
790 gaussjacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
791 {
792  typedef T value_type;
793 
794  Jacobi<N, T> p( a, b );
795 
796  roots( p, xr );
797 
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 ) );//gamma(a + m + 1);
802  value_type a3 = fact( b+value_type( N ) );//gamma(b + m + 1);
803  value_type a4 = fact( a+b+value_type( N ) );//gamma(a + b + m + 1);
804  value_type a5 = fact( value_type( N ) );
805  value_type a6 = a1 * a2 * a3;// / a4 / a5;
806 
807  for ( int k = 0; k < N; ++k )
808  {
809  value_type fp = p.derivate( xr[k] );
810  value_type dn = fp*fp*( 1.0 - xr[k]*xr[k] )*a4*a5;
811 
812  wr[k] =a6 / dn;
813  //wr[k] = a6 / ( 1.0 - xr[k]*xr[k]) / ( fp*fp );
814  }
815 }
816 
817 
823 template<int N, typename T, typename VectorW, typename VectorN>
824 void
825 gausslobattojacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
826 {
827  typedef T value_type;
828 
829  Jacobi<N-2, T> p( a+ 1.0, b+ 1.0 );
830 
831  Jacobi<N-1, T> q( a, b );
832 
833  VectorN prexr( N-2 );
834 
835  roots( p, prexr );
836 
837  xr[0]= -1.0;
838 
839  for ( int i = 1; i < N-1; ++i )
840  xr[i] = prexr[i-1];
841 
842  xr[N-1]= 1.0;
843 
844  const value_type two = 2.0;
845 
846  value_type a1 = math::pow( two,int( a+b+1.0 ) );
847  value_type a2 = fact( a+value_type( N ) -1.0 );//gamma(a + Q);
848  value_type a3 = fact( b+value_type( N ) -1.0 );//gamma(b + Q);
849 
850  value_type a4 = fact( a+b+value_type( N ) );//gamma(a + b + m + 1);
851  value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )-1.0 ); // (Q-1)!(Q-1)
852 
853  value_type a6 = a1 * a2 * a3;// Numerateur
854 
855  for ( int k = 0; k < N; ++k )
856  {
857  value_type fq = q.value( xr[k] );
858 
859  value_type dn = fq*fq*a4*a5;
860 
861  wr[k] =a6 /dn ;
862  }
863 
864  wr[0] = wr[0]*( b+1.0 );
865  wr[N-1] = wr[N-1]*( a+1.0 );
866 }
867 
868 namespace dyna
869 {
875 template<typename T, typename VectorW, typename VectorN>
876 void
877 gausslobattojacobi( size_type N, VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ), bool interior = false )
878 {
879  wr.resize( N - 2*interior );
880  xr.resize( N - 2*interior );
881 
882  typedef T value_type;
883 
884  Feel::dyna::Jacobi<T> p( N-2, a+ 1.0, b+ 1.0 );
885 
886  Feel::dyna::Jacobi<T> q( N-1, a, b );
887 
888  VectorN prexr( N-2 );
889 
890  roots( p, prexr );
891 
892  if ( !interior )
893  {
894  xr[0]= -1.0;
895  xr[N-1]= 1.0;
896  }
897 
898  for ( size_type i = 1-interior; i < N-( 1+interior ); ++i )
899  xr[i] = prexr[i-( 1-interior )];
900 
901 
902  const value_type two = 2.0;
903 
904  value_type a1 = math::pow( two,int( a+b+1.0 ) );
905  value_type a2 = fact( a+value_type( N ) -1.0 );//gamma(a + Q);
906  value_type a3 = fact( b+value_type( N ) -1.0 );//gamma(b + Q);
907 
908  value_type a4 = fact( a+b+value_type( N ) );//gamma(a + b + m + 1);
909  value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )-1.0 ); // (Q-1)!(Q-1)
910 
911  value_type a6 = a1 * a2 * a3;// Numerateur
912 
913  for ( int k = 1-interior; k < int( N-1-interior ); ++k )
914  {
915  value_type fq = q.value( xr[k] );
916 
917  value_type dn = fq*fq*a4*a5;
918 
919  wr[k] =a6 /dn ;
920  }
921 
922  if ( !interior )
923  {
924  wr[0] = wr[0]*( b+1.0 );
925  wr[N-1] = wr[N-1]*( a+1.0 );
926  }
927 }
928 
929 }
930 
937 template<int N, typename T, typename VectorW, typename VectorN>
938 void
939 left_gaussradaujacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
940 {
941  typedef T value_type;
942 
943  Jacobi<N-1, T> p( a, b+1.0 );
944 
945  Jacobi<N-1, T> q( a, b );
946 
947  VectorN prexr( N-1 );
948 
949  roots( p, prexr );
950 
951  xr[0]= -1.0;
952 
953  for ( int i = 1; i < N; ++i )
954  xr[i] = prexr[i-1];
955 
956  const value_type two = 2.0;
957 
958  value_type a1 = math::pow( two,int( a+b ) );
959  value_type a2 = fact( a+value_type( N ) -1.0 );//gamma(a + Q);
960  value_type a3 = fact( b+value_type( N ) -1.0 );//gamma(b + Q);
961 
962  value_type a4 = fact( a+b+value_type( N ) );//gamma(a + b + m + 1);
963  value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )+b ); // (Q-1)!(Q+b)
964 
965  value_type a6 = a1 * a2 * a3;// Numerateur
966 
967  for ( int k = 0; k < N; ++k )
968  {
969  value_type fq = q.value( xr[k] );
970 
971  value_type dn = fq*fq*a4*a5;
972 
973  wr[k] =a6*( 1.0-xr[k] ) /dn ;
974  }
975 
976  wr[0] = wr[0]*( b+1.0 );
977 }
978 
979 
980 template<int N, typename T, typename VectorW, typename VectorN>
981 void
982 right_gaussradaujacobi( VectorW& wr, VectorN& xr, T a = T( 0.0 ), T b = T( 0.0 ) )
983 {
984  typedef T value_type;
985 
986  Jacobi<N-1, T> p( a+1.0, b );
987 
988  Jacobi<N-1, T> q( a, b );
989 
990  VectorN prexr( N-1 );
991 
992  roots( p, prexr );
993 
994  for ( int i = 0; i < N-1; ++i )
995  xr[i] = prexr[i];
996 
997  xr[N-1]=1.0;
998 
999  const value_type two = 2.0;
1000 
1001  value_type a1 = math::pow( two,int( a+b ) );
1002  value_type a2 = fact( a+value_type( N ) -1.0 );//gamma(a + Q);
1003  value_type a3 = fact( b+value_type( N ) -1.0 );//gamma(b + Q);
1004 
1005  value_type a4 = fact( a+b+value_type( N ) );//gamma(a + b + m + 1);
1006  value_type a5 = fact( value_type( N )-1.0 )*( value_type( N )+a ); // (Q-1)!(Q+a)
1007 
1008  value_type a6 = a1 * a2 * a3;// numerator
1009 
1010  for ( int k = 0; k < N; ++k )
1011  {
1012  value_type fq = q.value( xr[k] );
1013 
1014  value_type dn = fq*fq*a4*a5;
1015 
1016  wr[k] =a6*( 1.0+xr[k] ) /dn ;
1017  }
1018 
1019  wr[N-1] = wr[N-1]*( a+1.0 );
1020 }
1021 
1022 
1023 template<int N, typename T>
1024 T integrate( boost::function<T( T const& )> const& f )
1025 {
1026  typedef T value_type;
1027  ublas::vector<T> xr( Jacobi<N, T>::nOrder );
1028  ublas::vector<T> wr( Jacobi<N, T>::nOrder );
1029 
1030  // get weights and nodes for Legendre polynomials
1031  details::gaussjacobi<N, T, ublas::vector<T> >( wr, xr );
1032 
1033  value_type res = 0.0;
1034 
1035  for ( int k = 0; k < Jacobi<N, T>::nOrder; ++k )
1036  res += wr[k]*f( xr[k] );
1037 
1038  return res;
1039 }
1040 } // details
1042 
1043 
1044 } // Feel
1045 #endif /* __Jacobi_H */

Generated on Sun Oct 20 2013 08:25:00 for Feel++ by doxygen 1.8.4