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
|
