Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : /// \cond
       9             : class DataVector;
      10             : class Matrix;
      11             : namespace Spectral {
      12             : enum class Quadrature : uint8_t;
      13             : }  // namespace Spectral
      14             : /// \endcond
      15             : 
      16             : namespace Spectral {
      17             : 
      18             : /*!
      19             :  * \ingroup SpectralGroup
      20             :  *
      21             :  * \brief A collection of helper functions for Chebyshev polynomials
      22             :  *
      23             :  * \details The Chebyshev polynomials are given by:
      24             :  * \f[
      25             :  * T_n (x) = \cos(n \theta), \text{      } \theta = \arccos(x)
      26             :  * \f]
      27             :  *
      28             :  * The Chebyshev expansion of a function \f$f \in [-1,1] \f$ is given by:
      29             :  * \f[
      30             :  * f(x) = \sum_{n=0}^{\infty} f_n T_n(x)
      31             :  * \f]
      32             :  * where
      33             :  * \f[
      34             :  * f_n = \frac{1}{c_n} \int_{-1}^1 f(x) T_n(x) w(x) dx
      35             :  * \f]
      36             :  * where the weight function is \f$w(x) = (1-x^2)^{-1/2}\f$ and the basis
      37             :  * function normalization value is given by:
      38             :  * \f{align*}{
      39             :  *   c_n
      40             :  *   &=\begin{cases}
      41             :  *     \hfil \pi \hfil & \text{if } n = 0 \\
      42             :  *     \hfil \frac{\pi}{2} \hfil & \text{otherwise}
      43             :  *   \end{cases}
      44             :  * \f}
      45             :  *
      46             :  * If a function is discretized at \f$N+1\f$ collocation points, the modal
      47             :  * representation will have \f$N+1\f$ spectral coefficients consisting of
      48             :  * \f[
      49             :  * f_n \qquad \text{for } n = 0, \ldots, N
      50             :  * \f]
      51             :  *
      52             :  * For more details about using Chebyshev polynomials see e.g. \cite Boyd2001
      53             :  * and \cite Fornberg1996.
      54             :  */
      55           1 : class Chebyshev {
      56             :  public:
      57             :   /*!
      58             :    * \brief Value of the basis function \f$\Phi_k(x) = T_k(x)\f$
      59             :    */
      60             :   template <typename T>
      61           1 :   static T basis_function_value(size_t k, const T& x);
      62             : 
      63             :   /*!
      64             :    * \brief The normalization square \f$c_k\f$ of the basis function
      65             :    * \f$\Phi_k(x)\f$, i.e. the definite integral of its square.
      66             :    */
      67           1 :   static double basis_function_normalization_square(size_t k);
      68             : 
      69             :   /*!
      70             :    * \brief Collocation points \f$\{x_i\}\f$
      71             :    *
      72             :    * \details The collocation points on the interval \f$[-1, 1]\f$ are given by
      73             :    * \f{align*}{
      74             :    *   x_i
      75             :    *   &=\begin{cases}
      76             :    *     \hfil - \cos \frac{(2i+1)\pi}{2N+2} \hfil &
      77             :    *       \text{for Quadrature::Gauss} \\
      78             :    *     \hfil - \cos \frac{i \pi}{N} \hfil &
      79             :    *       \text{for Quadrature::GaussLobatto}
      80             :    *   \end{cases}
      81             :    * \f}
      82             :    */
      83             :   template <Quadrature quadrature>
      84           1 :   static DataVector collocation_points(size_t num_points);
      85             : 
      86             :   /*!
      87             :    * \brief Integration weights \f$\{w_i\}\f$
      88             :    *
      89             :    * \details The integration weights are used to approximate the weighted
      90             :    * integral \f$Q[f]=\int f(x)w(x)\mathrm{d}x\approx \sum_k f_k w_k\f$.
      91             :    * For Quadrature::Gauss, the weights are given by:
      92             :    * \f[
      93             :    * w_i = \frac{\pi}{N+1}
      94             :    * \f]
      95             :    *
      96             :    * For Quadrature::GaussLobatto, the weights are given by:
      97             :    * \f{align*}{
      98             :    *   w_i
      99             :    *   &=\begin{cases}
     100             :    *     \hfil \frac{\pi}{2N}  \hfil & \text{for } j = 0, N\\
     101             :    *     \hfil \frac{\pi}{N}   \hfil & \text{for } j = 1, \ldots, N-1
     102             :    *   \end{cases}
     103             :    * \f}
     104             :    *
     105             :    * \note These weights are used to compute the modes from the nodal
     106             :    * values.  They are not used to evaluate definite or indefinite
     107             :    * integrals directly from the nodal values.
     108             :    */
     109             :   template <Quadrature quadrature>
     110           1 :   static DataVector integration_weights(size_t num_points);
     111             : 
     112             :   /*!
     113             :    * \brief Matrix used to compute the modes of the indefinite
     114             :    * integral from the modes of the integrand such that the constant
     115             :    * of integration is determined by requiring the integral to be zero
     116             :    * at \f$x=-1\f$.
     117             :    *
     118             :    * \details Chebyshev polynomials satisfy the identity:
     119             :    * \f[
     120             :    * \int^x dy T_n (y) = \left\{ \frac{T_{n+1}(x)}{2(n+1)} +
     121             :    *   \frac{T_{n-1}(x)}{2(n-1)} \right\}, \qquad n \geq 2
     122             :    * \f]
     123             :    *
     124             :    * Thus the modes \f$\tilde{f}_j\f$ of the integral are given as:
     125             :    * \f{align*}{
     126             :    *   \tilde{f}_i
     127             :    *   &=\begin{cases}
     128             :    *     \hfil \frac{f_{i-1}-f_{i+1}}{2i}, \hfil & \text{for } i > 1\\
     129             :    *     \hfil f_0 - \frac{f_2}{2} \hfil & \text{for } i = 1
     130             :    *   \end{cases}
     131             :    * \f}
     132             :    * where \f$f_j\f$ are the modes of the integrand.
     133             :    *\f$\tilde{f}_0\f$ is a constant of integration, which we choose such that
     134             :    * the integral is 0 at the left boundary of the domain (\f$x=-1\f$). The
     135             :    * condition for this is:
     136             :    *
     137             :    * \f[
     138             :    *   \tilde{f}_0=\sum_{i=1}^{N}(-1)^{i+1}\tilde{f}_i
     139             :    * \f]
     140             :    */
     141           1 :   static Matrix indefinite_integral_matrix(size_t num_points);
     142             : 
     143             :   /*!
     144             :    * \brief Row-vector used to compute the definite integral from the modes of
     145             :    * the integrand
     146             :    *
     147             :    * \details Given the modes \f$\tilde{f}_j\f$ of the integral (see
     148             :    * indefinite_integral_matrix), the definite integral is given by evaluating
     149             :    * the series at \f$x=1\f$ and \f$x=-1\f$ and taking the difference.  Given
     150             :    * that \f$T_n(1) = 1\f$ and \f$T_n(-1) = (-1)^n\f$, this means multiplying
     151             :    * the indefinite_integral_matrix by the row-vector \f$\{0, 2, 0, 2, \ldots
     152             :    * \}\f$ which yields:
     153             :    * \f{align*}{
     154             :    *   \tilde{q}_i
     155             :    *   &=\begin{cases}
     156             :    *     \hfil -\frac{2}{n^2-1}, \hfil & \text{for } i \text{ even}\\
     157             :    *     \hfil 0 \hfil & \text{for } i \text{ odd}
     158             :    *   \end{cases}
     159             :    * \f}
     160             :    */
     161           1 :   static Matrix definite_integral_matrix(size_t num_points);
     162             : };
     163             : }  // namespace Spectral