ROL
Lagrange.hpp
Go to the documentation of this file.
1#include<vector>
2#include<iostream>
3#ifndef __LAGRANGE__
4#define __LAGRANGE__
5
6template<class Real>
7class Lagrange{
8 public:
9 Lagrange(const std::vector<Real> &xin, const std::vector<Real> &xev );
10 ~Lagrange();
11
12 // Interpolate from the interpolation to the evaluation points
13 void interp(const std::vector<Real> &f, std::vector<Real> &p);
14
15 void dinterp(const std::vector<Real> &f, std::vector<Real> &p);
16
17 // Evaluate the kth interpolant on the evaluation points
18 void interpolant(const int k, std::vector<Real> &l);
19
20 // Derivative of the interpolating polynomial
21 void derivative(const int k, std::vector<Real> &d);
22
23 /* Implement sum formulas as found in equation 4.2 of Trefethen
24 and Berrut SIAM Review, Vol. 46, No. 3, pp.501-517 */
25 void bi_sum(const std::vector<Real> &f, std::vector<Real> &y);
26
27
28 private:
30 const std::vector<Real> xin_;
31
32 // \param xwv_ Vector of evaluation points
33 const std::vector<Real> xev_;
34
35 // \param nin_ Number of interpolation points
36 const int nin_;
37
38 // \param nev_ Number of evaluation points
39 const int nev_;
40
41 // \param w_ Vector of interpolation weights
42 std::vector<Real> w_;
43
44 // \param ell_ Vector conatining barycentric multiplicative polynomial on evaluation points
45 std::vector<Real> ell_;
46
47 // Pseudospectral differentiation matrix on interpolation points (column stacked)
48 std::vector<Real> D_;
49};
50
51
52
56template<class Real>
57Lagrange<Real>::Lagrange(const std::vector<Real> &xin, const std::vector<Real> &xev):
58 xin_(xin), xev_(xev), nin_(xin.size()), nev_(xev.size()), w_(nin_,0), ell_(nev_,0), D_(nin_*nin_,0) {
59
60 // For storing displacements between interpolation points
61 double d;
62
63 /* Compute the weights using as slightly modified version of the
64 algorithm on page 504 in the Trefethen & Berrut paper */
65
66 w_[0] = 1;
67
68 for(int j=1;j<nin_;++j)
69 {
70 w_[j] = 1.0;
71
72 for(int k=0;k<j;++k)
73 {
74 d = xin_[k]-xin_[j];
75 w_[k] *= d;
76 w_[j] *= -d;
77 }
78 }
79
80 for(int j=0;j<nin_;++j)
81 {
82 w_[j] = 1.0/w_[j];
83 }
84
85 std::vector<Real> ones(nin_,1.0); // Create vector of ones
86
87 this->bi_sum(ones,ell_); // Compute the \f$\ell(x)\f$ polynomial
88
89 for(int j=0;j<nev_;++j)
90 {
91 ell_[j] = 1.0/ell_[j];
92 }
93
94 for(int k=0;k<nin_;++k) {
95 for(int i=0;i<nin_;++i) {
96 if(i!=k) {
97 D_[i+k*nin_] = (w_[k]/w_[i])/(xin_[i]-xin_[k]);
98 }
99 }
100
101 }
102 for(int k=0;k<nin_;++k){
103 for(int i=0;i<nin_;++i) {
104 if(i!=k){
105 D_[k+k*nin_] -= D_[k+i*nin_];
106 }
107 }
108 }
109}
110
111
112
113
114template<class Real>
116
121template<class Real>
122void Lagrange<Real>::bi_sum(const std::vector<Real> &f, std::vector<Real> &y){
123
124 for(int j=0;j<nev_;++j)
125 {
126 y[j] = 0;
127 for(int k=0;k<nin_;++k)
128 {
129 if(xev_[j] == xin_[k])
130 {
131 y[j] = f[j];
132 break;
133 }
134 else
135 {
136 y[j] += w_[k]*f[k]/(xev_[j]-xin_[k]);
137 }
138 }
139 }
140}
141
146template<class Real>
147void Lagrange<Real>::interp(const std::vector<Real> &f, std::vector<Real> &p){
148 this->bi_sum(f,p);
149
150 for(int j=0;j<nev_;++j)
151 {
152 p[j] *= ell_[j];
153 }
154}
155
156template<class Real>
157void Lagrange<Real>::dinterp(const std::vector<Real> &f, std::vector<Real> &p){
158 std::vector<Real> fb(nin_,0);
159
160 std::copy(f.begin(),f.end(),fb.begin()+1);
161
162 this->bi_sum(fb,p);
163
164 for(int j=0;j<nev_;++j)
165 {
166 p[j] *= ell_[j];
167 }
168}
169
170
171
175template<class Real>
176void Lagrange<Real>::interpolant(const int k, std::vector<Real> &l){
177 std::vector<Real> f(nin_,0);
178 std::fill(l.begin(),l.end(),0);
179 f[k] = 1.0;
180 this->bi_sum(f,l);
181
182 for(int j=0;j<nev_;++j)
183 {
184 l[j] = ell_[j]*w_[k]/(xev_[j]-xin_[k]);
185 }
186}
187
188
190template<class Real>
191void Lagrange<Real>::derivative(const int k, std::vector<Real> &d ) {
192 std::vector<Real> lp(nin_,0);
193 std::fill(d.begin(),d.end(),0);
194 std::copy(D_.begin()+k*nin_,D_.begin()+(k+1)*nin_,lp.begin());
195
196 // Interpolate the derivative
197 this->interp(lp,d);
198}
199
200#endif
ROL::Ptr< OP > D_
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:57
const std::vector< Real > xin_
Definition Lagrange.hpp:30
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:147
void interpolant(const int k, std::vector< Real > &l)
Evaluate the kth interpolant on the evaluation points.
Definition Lagrange.hpp:176
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:122
std::vector< Real > w_
Definition Lagrange.hpp:42
std::vector< Real > ell_
Definition Lagrange.hpp:45
void dinterp(const std::vector< Real > &f, std::vector< Real > &p)
Definition Lagrange.hpp:157
const int nin_
Definition Lagrange.hpp:36
void derivative(const int k, std::vector< Real > &d)
Derivative of the th interpolant on the interpolation points.
Definition Lagrange.hpp:191
const int nev_
Definition Lagrange.hpp:39
std::vector< Real > D_
Definition Lagrange.hpp:48
const std::vector< Real > xev_
Definition Lagrange.hpp:33