NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
integration.hpp
Go to the documentation of this file.
1 /*********************************************************
2  * Martin Eckstein, 2010
3  * Gregory integration rules,
4  * and backward differntiation formulae
5  *********************************************************/
6 #ifndef INTEGRATION
7 #define INTEGRATION
8 
9 #include <cmath>
10 #include <cassert>
11 #include <iostream>
12 #include <complex>
13 #include <vector>
14 #include <stdlib.h>
15 
16 namespace integration{
17 
18 #define GREGORY_KMAX 8
19  void read_poly_interpolation(int k,long double *P); // P must be of size k1*k1
20  void read_poly_differentiation(int k,long double *w); // w must be of size k1*k1
21  void read_poly_integration(int k,long double *w); // w must be of size k1*k1*k1
22  void read_poly_interpolation(int k,double *P); // P must be of size k1*k1
23  void read_poly_differentiation(int k,double *P); // P must be of size k1*k1
24  void read_poly_integration(int k,double *P); // P must be of size k1*k1*k1
25  void read_poly_interpolation(int k,float *P); // P must be of size k1*k1
26  void read_poly_differentiation(int k,float *P); // P must be of size k1*k1
27  void read_poly_integration(int k,float *P); // P must be of size k1*k1*k1
28  void read_gregory_weights(int k,double *w);
29  void read_bd_weights(int k,double *w);
30  void read_rcorr(int k,double *w);
31 
32  /* #######################################################################################
33  #
34  # Integrator class
35  #
36  # The integrator contains all kinds of weights for integration and differntiation
37  # of a function at equally spaced abscissae, up to a certain order k of accuray
38  #
39  # !!!! OUT OF RANGE ERRORS ARE NOT CATCHED FOR THE ARRAYS BELOW !!!!
40  #
41  # (1) GREGORY WEIGHTS: I.gregory_weight and I.gregory_omega
42  #
43  # The Gregory integration formula of order k is defined by
44  #
45  # int_0^m f(x) dx = sum_{l=0}^{ max(k,m) } f_l I.gregory_weights(m,l)
46  #
47  # for m>2*k, the weights have a simpler form:
48  # gregory_weights(m,j)= I.gregory_omega(j) j<=k
49  # gregory_weights(m,j)= 1 k < j < m-k
50  # gregory_weights(m,j)= I.gregory_omega(m-j) j>=m-k
51  # I.gregory_omega returns for 0 <= j <= k
52  #
53  #
54  # (2) RCORR: I.rcorr
55  #
56  # These are weights for the special boundary term, which is needed for
57  # the convolution of matsubara Greenfunctions, which cannot be
58  # contuinuously be continued to x<0 and x>beta:
59  #
60  # For 0 < m < k:
61  # R=int_0^m dx a(m-x) b(x) = sum_{j,l=0}^k I.rcorr(m,j,l) a(j)*b(l)
62  #
63  # The formula should be used if both a(x) and b(x) are known only for
64  # x>0. The integral is computed by approximating a(m-x) in the
65  # interval [m-k,m] and b(x) in the interval [0,k], and integrating the
66  # product of the approxiating polynomials.
67  #
68  #
69  # (3) BACKWARD DIFFERENTIATION: I.bd_weight
70  #
71  # Backward differentiation formula !!!! OF ORDER K+1 !!!!
72  #
73  # d/dx f(x=x0) = sum_{l=0}^{k+1} I.bd_weights(l)*f(x0-l)
74  #
75  # The weights are related to the poly_differentiation weights (see 4a),
76  # I.bd_weights(l)=I._poly_differentiation(k+1,k+1-l)
77  # but (...k+1 are not stored)
78  #
79  # (4a) POLYNOMIAL INTERPOLATION: I.poly_interpolation
80  #
81  # Provided for completeness only:
82  # The k-th order approximation polynomial trough function
83  # values {(x_l,f(x_l)),l=0...k} is given by
84  #
85  # P(x) = sum_{alpha,l=0}^{k} x^alpha*I.interpolation(alpha,l)*f_l
86  #
87  # (4b) POLYNOMIAL DIFFERENTIATION: I.poly_differentiation
88  #
89  # To compute the derivative of the kth order polynomial
90  # through points (x_l,f_l), l=0...k:
91  #
92  # d/dx f(x=x_i) = sum_{l=0}^k I.poly_differentiation(i,l)*f(l)
93  #
94  #
95  # (4c) POLYNOMIAL INTEGRATION: I.poly_integration
96  #
97  # These are weights to compute the integral
98  #
99  # For 0 <= i,j <= k
100  # int_i^j dx f(x) = sum_{l=0}^k I.poly_integration(i,j,l)*f(l)
101  #
102  # This is needed for starting formulas of the volterra Integration.
103  # The integral is ontained by computing the approximating polynomial of
104  # f(x) in [0,k] and integrating it.
105  #
106  ####################################################################################### */
107 
108 
109 
128  template <typename T=double> class Integrator {
129 public:
147  Integrator(int k=0){
148  int k1=k+1;
149  if(k<0 || k > GREGORY_KMAX ){ std::cout << "Integrator: k out of range " << std::endl; abort();}
150  poly_interpolation_= new T [k1*k1];
151  poly_differentiation_= new T [k1*k1];
152  poly_integration_= new T [k1*k1*k1];
153  bd_weights_= new T [k1+1];
154  gregory_weights_= new T [ 4*k1*k1 ];
155  gregory_omega_=gregory_weights_ + ((2*k+1)*2*k1); /*pointer to last line*/
156  if(k>1) rcorr_ = new T [ (k-1)*k1*k1] ; else rcorr_=0;
157  k_=k;
158  integration::read_poly_interpolation(k_,poly_interpolation_);
159  integration::read_poly_differentiation(k_,poly_differentiation_);
160  integration::read_poly_integration(k_,poly_integration_);
161  integration::read_bd_weights(k+1,bd_weights_);
162  integration::read_gregory_weights(k,gregory_weights_);
163  integration::read_rcorr(k_,rcorr_);
164 }
166  int k1;
167  if( this->k_==I.k_ ) return *this;
168  delete [] poly_interpolation_;
169  delete [] poly_differentiation_;
170  delete [] poly_integration_;
171  delete [] bd_weights_;
172  delete [] gregory_weights_;
173  delete [] rcorr_;
174  k_=I.k_;
175  k1=k_+1;
176  poly_interpolation_= new T [k1*k1];
177  poly_differentiation_= new T [k1*k1];
178  poly_integration_= new T [k1*k1*k1];
179  bd_weights_= new T [k1+1];
180  gregory_weights_= new T [ 4*k1*k1 ];
181  gregory_omega_=gregory_weights_ + ((2*k_+1)*2*k1); /*pointer to last line*/
182  if(k_>1) rcorr_ = new T [ (k_-1)*k1*k1] ; else rcorr_=0;
183  integration::read_poly_interpolation(k_,poly_interpolation_);
184  integration::read_poly_differentiation(k_,poly_differentiation_);
185  integration::read_poly_integration(k_,poly_integration_);
186  integration::read_bd_weights(k_+1,bd_weights_);
187  integration::read_gregory_weights(k_,gregory_weights_);
188  integration::read_rcorr(k_,rcorr_);
189  return *this;
190 }
192  delete [] poly_interpolation_;
193  delete [] poly_differentiation_;
194  delete [] poly_integration_;
195  delete [] bd_weights_;
196  delete [] gregory_weights_;
197  delete [] rcorr_;
198 }
200  int k1;
201  k_=I.k_;
202  k1=k_+1;
203  poly_interpolation_= new T [k1*k1];
204  poly_differentiation_= new T [k1*k1];
205  poly_integration_= new T [k1*k1*k1];
206  bd_weights_= new T [k1+1];
207  gregory_weights_= new T [ 4*k1*k1 ];
208  gregory_omega_=gregory_weights_ + ((2*k_+1)*2*k1); /*pointer to last line*/
209  if(k_>1) rcorr_ = new T [ (k_-1)*k1*k1] ; else rcorr_=0;
210  integration::read_poly_interpolation(k_,poly_interpolation_);
211  integration::read_poly_differentiation(k_,poly_differentiation_);
212  integration::read_poly_integration(k_,poly_integration_);
213  integration::read_bd_weights(k_+1,bd_weights_);
214  integration::read_gregory_weights(k_,gregory_weights_);
215  integration::read_rcorr(k_,rcorr_);
216 }
226  int get_k(void) {return k_;}
236  int k(void) {return k_;}
237 
262  T poly_interpolation(int alpha,int l) {return poly_interpolation_[alpha*(k_+1)+l];}
263 
264 
290  T poly_differentiation(int i,int l) {return poly_differentiation_[i*(k_+1)+l];}
291 
321  T poly_integration(int i,int j,int l) {return poly_integration_[i*(k_+1)*(k_+1) + j*(k_+1) + l];}
322 
323 
346  T bd_weights(int l) {return bd_weights_[l];}
347 
382  T gregory_weights(int n,int j){
383  int k1=k_+1,k2=2*k1;
384  assert( j>=0 && ( j<=k_ || j<=n) );
385  if( n>=k2 ){
386  if( j > k_ && j < n-k_ ) { return 1;}
387  else if ( j<k1 ) { return gregory_omega_[j];}
388  else { return gregory_omega_[n-j]; }
389 }else{ return gregory_weights_[n*k2 + j];}
390 }
391 
405  T gregory_omega(int j) {return gregory_omega_[j];}
406 
407 
439  T rcorr(int m,int j,int l){
440  int k1=k_+1;
441  if(k_<2) return 0;
442  assert(m<k_ && j<=k_ && l<=k_);
443  return *(rcorr_ + (m-1)*k1*k1 + j*k1 +l);
444 }
445 
446  // The following is a particulary useless function,
447  // more or less provided as an example
448  // compute the integral int_j^n dx of a function f(x) -
449  // in the class M, M=0 must make sense :
450 
474  template <class M> M integrate(const std::vector<M> &f,int n){
475  int k1=k_+1,i,k2=2*k1,n1=n-k_;
476  T *w;
477  M sum=f[0];
478  sum=0;
479 
480  n1=(n<=k_ ? k_ : n);
481  assert( int(f.size()) >= n1);
482 
483  if( n>=k2 ){
484  n1=n-k_;
485  for(i=0;i<=k_;i++) sum += f[i]*gregory_omega_[i];
486  for(i=k1;i<n1;i++) sum += f[i];
487  for(i=n1;i<=n;i++) sum += f[i]*gregory_omega_[n-i];
488 }else if(n>0){
489  w=gregory_weights_ + n*k2;
490  for(i=0;i<=n1;i++) sum += f[i]*w[i];
491 }
492  return sum;
493 }
494 
495 private:
496  int k_;
497  T *poly_interpolation_;
498  T *poly_differentiation_;
499  T *poly_integration_;
500  T *bd_weights_;
501  T *gregory_weights_;
502  T *gregory_omega_;
503  T *rcorr_;
504 };
505 
506  template<typename T> Integrator<T> &I(int k){
507  static Integrator<T> Isave[6];
508  static int init=0;
509  if(!init){
510  init=1;
511  for(int k=1;k<=5;k++) Isave[k]=Integrator<T>(k);
512 }
513  return Isave[k%6];
514 }
515 
516 #ifndef CNTR_NO_EXTERN_TEMPLATES
517  extern template class Integrator<double>;
518  extern template Integrator<double> &I<double>(int k);
519 #endif
520 
521 } //namespace
522 
523 #endif
524 
525 
Class Integrator contains all kinds of weights for integration and differentiation of a function at ...
T rcorr(int m, int j, int l)
Returns the special quadrature weights for computing integrals on the Matsubara axis.
int get_k(void)
Returns the order of the integrator class.
Integrator(int k=0)
Initializes the Integrator class for a given order k.
Integrator< T > & I(int k)
T gregory_omega(int j)
Returns the Gregory weights at the integral boundaries (see gregory_weights).
T poly_integration(int i, int j, int l)
Returns the the weight needed for polynomial integration.
T gregory_weights(int n, int j)
Returns the Gregory weights for integration.
T bd_weights(int l)
Returns the the backwards differencing coefficients.
T poly_interpolation(int alpha, int l)
Returns the the weight needed for polynomial interpolation.
T poly_differentiation(int i, int l)
Returns the the weight needed for polynomial differentiation.
int k(void)
Returns the order of the integrator class.
template Integrator< double > & I< double >(int k)
Integrator & operator=(const Integrator &I)
M integrate(const std::vector< M > &f, int n)
Integrates a function using Gregory quadrature.
Integrator(const Integrator &I)