BarycentricRational.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <cstddef>
7 #include <unistd.h> // IWYU pragma: keep
8 #include <vector>
9
10 namespace PUP {
11 class er;
12 } // namespace PUP
13
14 namespace intrp {
15 /*!
16  * \ingroup NumericalAlgorithmsGroup
17  * \brief A barycentric rational interpolation class
18  *
19  * The class builds a barycentric rational interpolant of a specified order
20  * using the x_values and y_values passed into the constructor.
21  * Barycentric interpolation requires \f$3N\f$ storage, and costs
22  * \f$\mathcal{O}(Nd)\f$ to construct, where \f$N\f$ is the size of the x- and
23  * y-value vectors and \f$d\f$ is the order of the interpolant. The evaluation
24  * cost is \f$\mathcal{O}(N)\f$ compared to \f$\mathcal{O}(d)\f$ of a spline
25  * method, but constructing the barycentric interpolant does not require any
26  * derivatives of the function to be known.
27  *
28  * The interpolation function is
29  *
30  * \f[
31  * \mathcal{I}(x)=\frac{\sum_{i=0}^{N-1}w_i y_i /
32  * (x-x_i)}{\sum_{i=0}^{N-1}w_i/(x-x_i)}
33  * \f]
34  *
35  * where \f$w_i\f$ are the weights. The weights are computed using
36  *
37  * \f[
38  * w_k=\sum_{i=k-d\\0\le i < N-d}^{k}(-1)^{i}
39  * \prod_{j=i\\j\ne k}^{i+d}\frac{1}{x_k-x_j} \f]
40  *
41  * \requires x_values.size() == y_values.size() and
42  * x_values_.size() >= order
43  */
45  public:
46  BarycentricRational() noexcept = default;
48  std::vector<double> y_values, size_t order) noexcept;
49
50  double operator()(double x_to_interp_to) const noexcept;
51
52  size_t order() const noexcept;
53
54  // clang-tidy: no runtime references
55  void pup(PUP::er& p) noexcept; // NOLINT
56
57  private:
58  std::vector<double> x_values_;
59  std::vector<double> y_values_;
60  std::vector<double> weights_;
61  ssize_t order_{0};
62 };
63 } // namespace intrp
Definition: Strahlkorper.hpp:14