ROL
Lagrange.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include<vector>
11 #include<iostream>
12 #ifndef __LAGRANGE__
13 #define __LAGRANGE__
14 
15 template<class Real>
16 class Lagrange{
17  public:
18  Lagrange(const std::vector<Real> &xin, const std::vector<Real> &xev );
19  ~Lagrange();
20 
21  // Interpolate from the interpolation to the evaluation points
22  void interp(const std::vector<Real> &f, std::vector<Real> &p);
23 
24  void dinterp(const std::vector<Real> &f, std::vector<Real> &p);
25 
26  // Evaluate the kth interpolant on the evaluation points
27  void interpolant(const int k, std::vector<Real> &l);
28 
29  // Derivative of the interpolating polynomial
30  void derivative(const int k, std::vector<Real> &d);
31 
32  /* Implement sum formulas as found in equation 4.2 of Trefethen
33  and Berrut SIAM Review, Vol. 46, No. 3, pp.501-517 */
34  void bi_sum(const std::vector<Real> &f, std::vector<Real> &y);
35 
36 
37  private:
39  const std::vector<Real> xin_;
40 
41  // \param xwv_ Vector of evaluation points
42  const std::vector<Real> xev_;
43 
44  // \param nin_ Number of interpolation points
45  const int nin_;
46 
47  // \param nev_ Number of evaluation points
48  const int nev_;
49 
50  // \param w_ Vector of interpolation weights
51  std::vector<Real> w_;
52 
53  // \param ell_ Vector conatining barycentric multiplicative polynomial on evaluation points
54  std::vector<Real> ell_;
55 
56  // Pseudospectral differentiation matrix on interpolation points (column stacked)
57  std::vector<Real> D_;
58 };
59 
60 
61 
65 template<class Real>
66 Lagrange<Real>::Lagrange(const std::vector<Real> &xin, const std::vector<Real> &xev):
67  xin_(xin), xev_(xev), nin_(xin.size()), nev_(xev.size()), w_(nin_,0), ell_(nev_,0), D_(nin_*nin_,0) {
68 
69  // For storing displacements between interpolation points
70  double d;
71 
72  /* Compute the weights using as slightly modified version of the
73  algorithm on page 504 in the Trefethen & Berrut paper */
74 
75  w_[0] = 1;
76 
77  for(int j=1;j<nin_;++j)
78  {
79  w_[j] = 1.0;
80 
81  for(int k=0;k<j;++k)
82  {
83  d = xin_[k]-xin_[j];
84  w_[k] *= d;
85  w_[j] *= -d;
86  }
87  }
88 
89  for(int j=0;j<nin_;++j)
90  {
91  w_[j] = 1.0/w_[j];
92  }
93 
94  std::vector<Real> ones(nin_,1.0); // Create vector of ones
95 
96  this->bi_sum(ones,ell_); // Compute the \f$\ell(x)\f$ polynomial
97 
98  for(int j=0;j<nev_;++j)
99  {
100  ell_[j] = 1.0/ell_[j];
101  }
102 
103  for(int k=0;k<nin_;++k) {
104  for(int i=0;i<nin_;++i) {
105  if(i!=k) {
106  D_[i+k*nin_] = (w_[k]/w_[i])/(xin_[i]-xin_[k]);
107  }
108  }
109 
110  }
111  for(int k=0;k<nin_;++k){
112  for(int i=0;i<nin_;++i) {
113  if(i!=k){
114  D_[k+k*nin_] -= D_[k+i*nin_];
115  }
116  }
117  }
118 }
119 
120 
121 
122 
123 template<class Real>
125 
130 template<class Real>
131 void Lagrange<Real>::bi_sum(const std::vector<Real> &f, std::vector<Real> &y){
132 
133  for(int j=0;j<nev_;++j)
134  {
135  y[j] = 0;
136  for(int k=0;k<nin_;++k)
137  {
138  if(xev_[j] == xin_[k])
139  {
140  y[j] = f[j];
141  break;
142  }
143  else
144  {
145  y[j] += w_[k]*f[k]/(xev_[j]-xin_[k]);
146  }
147  }
148  }
149 }
150 
155 template<class Real>
156 void Lagrange<Real>::interp(const std::vector<Real> &f, std::vector<Real> &p){
157  this->bi_sum(f,p);
158 
159  for(int j=0;j<nev_;++j)
160  {
161  p[j] *= ell_[j];
162  }
163 }
164 
165 template<class Real>
166 void Lagrange<Real>::dinterp(const std::vector<Real> &f, std::vector<Real> &p){
167  std::vector<Real> fb(nin_,0);
168 
169  std::copy(f.begin(),f.end(),fb.begin()+1);
170 
171  this->bi_sum(fb,p);
172 
173  for(int j=0;j<nev_;++j)
174  {
175  p[j] *= ell_[j];
176  }
177 }
178 
179 
180 
184 template<class Real>
185 void Lagrange<Real>::interpolant(const int k, std::vector<Real> &l){
186  std::vector<Real> f(nin_,0);
187  std::fill(l.begin(),l.end(),0);
188  f[k] = 1.0;
189  this->bi_sum(f,l);
190 
191  for(int j=0;j<nev_;++j)
192  {
193  l[j] = ell_[j]*w_[k]/(xev_[j]-xin_[k]);
194  }
195 }
196 
197 
199 template<class Real>
200 void Lagrange<Real>::derivative(const int k, std::vector<Real> &d ) {
201  std::vector<Real> lp(nin_,0);
202  std::fill(d.begin(),d.end(),0);
203  std::copy(D_.begin()+k*nin_,D_.begin()+(k+1)*nin_,lp.begin());
204 
205  // Interpolate the derivative
206  this->interp(lp,d);
207 }
208 
209 #endif
void bi_sum(const std::vector< Real > &f, std::vector< Real > &y)
This routine evaluates sums of the form shown in equation (4.2) in the paper by J-P Berrut and L...
Definition: Lagrange.hpp:131
void interpolant(const int k, std::vector< Real > &l)
Evaluate the kth interpolant on the evaluation points.
Definition: Lagrange.hpp:185
ROL::Ptr< OP > D_
const int nev_
Definition: Lagrange.hpp:48
const std::vector< Real > xin_
Definition: Lagrange.hpp:39
void derivative(const int k, std::vector< Real > &d)
Derivative of the th interpolant on the interpolation points.
Definition: Lagrange.hpp:200
const int nin_
Definition: Lagrange.hpp:45
const std::vector< Real > xev_
Definition: Lagrange.hpp:42
std::vector< Real > D_
Definition: Lagrange.hpp:57
std::vector< Real > ell_
Definition: Lagrange.hpp:54
void interp(const std::vector< Real > &f, std::vector< Real > &p)
Given the values of a function on the interpolation points xin, stored in f, evaluate the interpolati...
Definition: Lagrange.hpp:156
Lagrange(const std::vector< Real > &xin, const std::vector< Real > &xev)
Interpolation object which interpolates from to the grid xin to xev.
Definition: Lagrange.hpp:66
void dinterp(const std::vector< Real > &f, std::vector< Real > &p)
Definition: Lagrange.hpp:166
std::vector< Real > w_
Definition: Lagrange.hpp:51