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 Matrix;
      10             : template <size_t>
      11             : class Mesh;
      12             : namespace Spectral {
      13             : enum class Basis : uint8_t;
      14             : enum class Quadrature : uint8_t;
      15             : }  // namespace Spectral
      16             : /// \endcond
      17             : 
      18             : namespace Spectral {
      19             : /*!
      20             :  * \brief %Matrix used to perform an indefinite integral of a function over the
      21             :  * logical grid. The left boundary condition is such that the integral is 0 at
      22             :  * \f$\xi=-1\f$
      23             :  *
      24             :  * Currently only Legendre and Chebyshev polynomials are implemented, but we
      25             :  * provide a derivation for how to compute the indefinite integration matrix for
      26             :  * general Jacobi polynomials.
      27             :  *
      28             :  * #### Legendre Polynomials
      29             :  * The Legendre polynomials have the identity:
      30             :  *
      31             :  * \f{align*}{
      32             :  * P_n(x) = \frac{1}{2n+1}\frac{d}{dx}\left(P_{n+1}(x) - P_{n-1}(x)\right)
      33             :  * \f}
      34             :  *
      35             :  * The goal is to evaluate the integral of a function \f$u\f$ expanded in terms
      36             :  * of Legendre polynomials as
      37             :  *
      38             :  * \f{align*}{
      39             :  * u(x) = \sum_{i=0}^{N} c_i P_i(x)
      40             :  * \f}
      41             :  *
      42             :  * We similarly expand the indefinite integral of \f$u\f$ as
      43             :  *
      44             :  * \f{align*}{
      45             :  * \left.\int u(y) dy\right\rvert_{x}=&\sum_{i=0}^N \tilde{c}_i P_i(x) \\
      46             :  *   =&\left.\int\sum_{i=1}^{N}\frac{c_i}{2i+1}
      47             :  *      \left(P_{i+1}(y)-P_{i-1}(y)\right)dy\right\rvert_{x}
      48             :  *      + \tilde{c}_0 P_0(x) \\
      49             :  *   =&\sum_{i=1}^{N}\left(\frac{c_{i-1}}{2i-1} - \frac{c_{i+1}}{2i+3}\right)
      50             :  *     P_i(x) + \tilde{c}_0 P_0(x)
      51             :  * \f}
      52             :  *
      53             :  * Thus we get that for \f$i>0\f$
      54             :  *
      55             :  * \f{align*}{
      56             :  * \tilde{c}_i=\frac{c_{i-1}}{2i-1}-\frac{c_{i+1}}{2i+3}
      57             :  * \f}
      58             :  *
      59             :  * and \f$\tilde{c}_0\f$ is a constant of integration, which we choose such that
      60             :  * the integral is 0 at the left boundary of the domain (\f$x=-1\f$). The
      61             :  * condition for this is:
      62             :  *
      63             :  * \f{align*}{
      64             :  *   \tilde{c}_0=\sum_{i=1}^{N}(-1)^{i+1}\tilde{c}_i
      65             :  * \f}
      66             :  *
      67             :  * The matrix returned by this function is the product of the tridiagonal matrix
      68             :  * for the \f$\tilde{c}_i\f$ and the matrix for the boundary condition.
      69             :  *
      70             :  * #### Chebyshev Polynomials
      71             :  *
      72             :  * A similar derivation leads to the relations:
      73             :  *
      74             :  * \f{align*}{
      75             :  *  \tilde{c}_i=&\frac{c_{i-1}-c_{i+1}}{2i},&\mathrm{if}\;i>1 \\
      76             :  *  \tilde{c}_1=&c_0 - \frac{c_2}{2},&\mathrm{if}\;i=1 \\
      77             :  * \f}
      78             :  *
      79             :  * We again have:
      80             :  *
      81             :  * \f{align*}{
      82             :  * \tilde{c}_0=\sum_{i=1}^N(-1)^{i+1}\tilde{c}_i
      83             :  * \f}
      84             :  *
      85             :  * These are then used to define the indefinite integration matrix.
      86             :  *
      87             :  * #### Jacobi Polynomials
      88             :  *
      89             :  * For general Jacobi polynomials \f$P^{(\alpha,\beta)}_n(x)\f$ given by
      90             :  *
      91             :  * \f{align*}{
      92             :  *  (1-x)^\alpha(1+x)^\beta P^{(\alpha,\beta)}_n(x)=\frac{(-1)^n}{2^n n!}
      93             :  *  \frac{d^n}{dx^n}\left[(1-x)^{\alpha+n}(1+x)^{\beta+n}\right]
      94             :  * \f}
      95             :  *
      96             :  * we have that
      97             :  *
      98             :  * \f{align*}{
      99             :  * P^{(\alpha,\beta)}_n(x)=\frac{d}{dx}\left(
     100             :  * b^{(\alpha,\beta)}_{n-1,n}P^{(\alpha,\beta)}_{n-1}(x) +
     101             :  * b^{(\alpha,\beta)}_{n,n}P^{(\alpha,\beta)}_n(x) +
     102             :  * b^{(\alpha,\beta)}_{n+1,n}P^{(\alpha,\beta)}_{n+1}(x)
     103             :  * \right)
     104             :  * \f}
     105             :  *
     106             :  * where
     107             :  *
     108             :  * \f{align*}{
     109             :  * b^{(\alpha,\beta)}_{n-1,n}=&-\frac{1}{n+\alpha+\beta}
     110             :  *                            a^{(\alpha,\beta)}_{n-1,n} \\
     111             :  * b^{(\alpha,\beta)}_{n,n}=&-\frac{2}{\alpha+\beta}
     112             :  *                          a^{(\alpha,\beta)}_{n,n} \\
     113             :  * b^{(\alpha,\beta)}_{n+1,n}=&\frac{1}{n+1}
     114             :  *                            a^{(\alpha,\beta)}_{n+1,n} \\
     115             :  * a^{(\alpha,\beta)}_{n-1,n}=&\frac{2(n+\alpha)(n+\beta)}
     116             :  *            {(2n+\alpha+\beta+1)(2n+\alpha+\beta)} \\
     117             :  * a^{(\alpha,\beta)}_{n,n}=&-\frac{\alpha^2-\beta^2}
     118             :  *            {(2n+\alpha+\beta+2)(2n+\alpha+\beta)} \\
     119             :  * a^{(\alpha,\beta)}_{n-1,n}=&\frac{2(n+1)(n+\alpha+\beta+1)}
     120             :  *            {(2n+\alpha+\beta+2)(2n+\alpha+\beta+1)}
     121             :  * \f}
     122             :  *
     123             :  * Following the same derivation we get that
     124             :  *
     125             :  * \f{align*}{
     126             :  *   \tilde{c}_i=c_{i+1}b^{(\alpha,\beta)}_{i,i+1}
     127             :  *              +c_i b^{(\alpha,\beta)}_{i,i}
     128             :  *              +c_{i-1}b^{(\alpha,\beta)}_{i,i-1}
     129             :  * \f}
     130             :  *
     131             :  * and the boundary condition is
     132             :  *
     133             :  * \f{align*}{
     134             :  *  \tilde{c}_0=\sum_{i=1}^N(-1)^{i+1}
     135             :  *              \frac{\Gamma(i+\alpha+1)}{i!\Gamma(\alpha+1)} \tilde{c}_i
     136             :  * \f}
     137             :  *
     138             :  * where \f$\Gamma(x)\f$ is the Gamma function.
     139             :  */
     140             : template <Basis BasisType, Quadrature QuadratureType>
     141           1 : const Matrix& integration_matrix(size_t num_points);
     142             : 
     143             : /*!
     144             :  * \brief Indefinite integration matrix for a one-dimensional mesh.
     145             :  *
     146             :  * \see integration_matrix(size_t)
     147             :  */
     148           1 : const Matrix& integration_matrix(const Mesh<1>& mesh);
     149             : }  // namespace Spectral