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
15template<class Real>
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
65template<class Real>
66Lagrange<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
123template<class Real>
125
130template<class Real>
131void 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
155template<class Real>
156void 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
165template<class Real>
166void 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
184template<class Real>
185void 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
199template<class Real>
200void 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
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:66
const std::vector< Real > xin_
Definition Lagrange.hpp:39
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
void interpolant(const int k, std::vector< Real > &l)
Evaluate the kth interpolant on the evaluation points.
Definition Lagrange.hpp:185
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
std::vector< Real > w_
Definition Lagrange.hpp:51
std::vector< Real > ell_
Definition Lagrange.hpp:54
void dinterp(const std::vector< Real > &f, std::vector< Real > &p)
Definition Lagrange.hpp:166
const int nin_
Definition Lagrange.hpp:45
void derivative(const int k, std::vector< Real > &d)
Derivative of the th interpolant on the interpolation points.
Definition Lagrange.hpp:200
const int nev_
Definition Lagrange.hpp:48
std::vector< Real > D_
Definition Lagrange.hpp:57
const std::vector< Real > xev_
Definition Lagrange.hpp:42