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 : /// \endcond
12 :
13 : namespace Spectral {
14 :
15 : /*!
16 : * \ingroup SpectralGroup
17 : *
18 : * \brief A collection of helper functions for using a Fourier series
19 : *
20 : * \details The general real-valued Fourier series is given by:
21 : * \f[
22 : * f(x) = a_0 + \sum_{n=1}^{\infty} a_n \cos(nx) + b_n \sin(nx)
23 : * \f]
24 : * where
25 : * \f{align*}
26 : * a_0 &= \frac{1}{2\pi} \int_0^{2\pi} f(x) dx \\
27 : * a_n &= \frac{1}{\pi} \int_0^{2\pi} f(x) \cos(nx) dx \\
28 : * b_n &= \frac{1}{\pi} \int_0^{2\pi} f(x) \sin(nx) dx
29 : * \f}
30 : *
31 : * If a function is discretized at \f$N\f$ collocation points, the modal
32 : * representation will have \f$N\f$ spectral coefficients consisting of
33 : * \f{align*}
34 : * a_n & \qquad \text{for } n = 0, \ldots, \frac{N}{2} \\
35 : * b_n & \qquad \text{for } n = 1, \ldots, \frac{N-1}{2}
36 : * \f}
37 : * Thus to represent all terms through the harmonic \f$M\f$ requires \f$N = 2M
38 : * +1\f$ collocation points.
39 : *
40 : * For more details about using Fourier series see e.g. \cite Boyd2001 and
41 : * \cite Fornberg1996.
42 : */
43 1 : class Fourier {
44 : public:
45 : /*!
46 : * \brief Value of the basis function \f$\Phi_k(x)\f$
47 : *
48 : * \details We define the basis functions as
49 : * \f{align*}{
50 : * \Phi_k(x)
51 : * &=\begin{cases}
52 : * \hfil \cos(kx) \hfil & \text{if } k \geq 0 \\
53 : * \hfil \sin(-kx) \hfil & \text{if } k < 0
54 : * \end{cases}
55 : * \f}
56 : */
57 : template <typename T>
58 1 : static T basis_function_value(int k, const T& x);
59 :
60 : /*!
61 : * \brief The normalization square \f$c_k\f$ of the basis function
62 : * \f$\Phi_k(x)\f$, i.e. the definite integral of its square.
63 : *
64 : * \details
65 : * In particular,
66 : * \f{align*}{
67 : * c_k
68 : * &=\begin{cases}
69 : * \hfil 2 \pi \hfil & \text{if } k = 0 \\
70 : * \hfil \pi \hfil & \text{otherwise}
71 : * \end{cases}
72 : * \f}
73 : */
74 1 : static double basis_function_normalization_square(int k);
75 :
76 : /*!
77 : * \brief Collocation points \f$\{x_i\}\f$
78 : *
79 : * \details The collocation points on the interval \f$0 \leq x < 2 \pi\f$ are
80 : * given by
81 : * \f[
82 : * x_i = \frac{2 \pi i}{N}
83 : * \f]
84 : */
85 1 : static DataVector collocation_points(size_t num_points);
86 :
87 : /*!
88 : * \brief Quadrature weights \f$\{w_i\}\f$
89 : *
90 : * \details The quadrature weights are given by
91 : * \f[
92 : * w_i = \frac{2 \pi}{N}
93 : * \f]
94 : */
95 1 : static DataVector quadrature_weights(size_t num_points);
96 :
97 : /*!
98 : * \brief Storage index (offset) \f$i\f$ into a ModalVector (representing the
99 : * coefficients) of the given mode \f$k\f$
100 : *
101 : * \details The modal coefficients are stored in a ModalVector as
102 : * \f$\{u_0, u_1, u_{-1}, u_2, u_{-2}, \ldots, u_M, u_{-M}\}\f$,
103 : * where \f$u_{-M}\f$ is omitted if the number of coefficients is even.
104 : * Therefore the storage index \f$i\f$ for the mode \f$u_k\f$ is:
105 : * \f{align*}{
106 : * i
107 : * &=\begin{cases}
108 : * \hfil 0 \hfil & \text{if } k = 0 \\
109 : * \hfil 2k-1 \hfil & \text{if } k > 0 \\
110 : * \hfil -2k \hfil & \text{if } k < 0
111 : * \end{cases}
112 : * \f}
113 : */
114 1 : static size_t modal_storage_index(int k);
115 :
116 : /*!
117 : * \brief Mode \f$k\f$ corresponding to given storage index (offset) \f$i\f$
118 : * into a ModalVector (representing the coefficients)
119 : *
120 : * \details The modal coefficients are stored in a ModalVector as
121 : * \f$\{u_0, u_1, u_{-1}, u_2, u_{-2}, \ldots, u_M, u_{-M}\}\f$,
122 : * where \f$u_{-M}\f$ is omitted if the number of coefficients is even.
123 : * Therefore the storage index \f$i\f$ for the mode \f$u_k\f$ is:
124 : * \f{align*}{
125 : * k
126 : * &=\begin{cases}
127 : * \hfil 0 \hfil & \text{if } i = 0 \\
128 : * \hfil 1 + \frac{i}{2} \hfil & \text{if } i \text{ is odd} \\
129 : * \hfil -\frac{i}{2} \hfil & \text{if } i \text{ is even}
130 : * \end{cases}
131 : * \f}
132 : */
133 1 : static int mode_at_storage_index(size_t storage_index);
134 :
135 : /*!
136 : * \brief Matrix \f$D_{i,j}\f$ used to obtain the first derivative
137 : *
138 : * \details The differentiation matrix is given by:
139 : * \f{align*}{
140 : * D_{i,j}
141 : * &=\begin{cases}
142 : * \hfil 0 \hfil & \text{if } i = j \\
143 : * \hfil \frac{1}{2} (-1)^{i-j} \csc \left[0.5 (x_i - x_j)\right] \hfil &
144 : * \text{if } i \neq j,\; N \text{ is odd} \\
145 : * \hfil \frac{1}{2} (-1)^{i-j} \cot \left[0.5 (x_i - x_j)\right] \hfil &
146 : * \text{if } i \neq j,\; N \text{ is even}
147 : * \end{cases}
148 : * \f}
149 : */
150 1 : static Matrix differentiation_matrix(size_t num_points);
151 :
152 : /*!
153 : * \brief %Matrix used to interpolate to the \p target_points.
154 : *
155 : * \details Each row of the matrix is given by the interpolation weights for
156 : * interpolating to a particular target point \f$x\f$. At each target point,
157 : * \f$x\f$, the interpolation weights are given by:
158 : * \f{align*}{
159 : * C_j(x)
160 : * &=\begin{cases}
161 : * \hfil \frac{1}{N} \sin \left[0.5 N(x - x_j)\right]
162 : * \csc \left[0.5 (x - x_j)\right] \hfil & \text{if } N \text{ is odd} \\
163 : * \hfil \frac{1}{N} \sin \left[0.5 N (x - x_j)\right]
164 : * \cot \left[0.5 (x - x_j)\right] \hfil & \text{if } N \text{ is even}
165 : * \end{cases}
166 : * \f}
167 : * where \f$x_j\f$ are the collocation_points.
168 : */
169 : template <typename T>
170 1 : static Matrix interpolation_matrix(size_t num_points, const T& target_points);
171 : };
172 : } // namespace Spectral
|
