Linear model reduction using SPOD modes

Space-only POD aims to reconstruct snapshots of the state by adding together prominent modes weighted by expansion coefficients. The first mode is defined to maximize $\lambda[\bm{\phi}({\bf{x}})]$, the expected value of the energy captured by it,

	$$\lambda[\bm{\phi}({\bf{x}})]=\frac{\mathbb{E}\big{[}|\langle\bm{q}({{\bf x}}),\bm{\phi}({{\bf x}})\rangle_{{\bf x}}|^{2}\big{]}}{\|\bm{\phi}({{\bf x}})\|^{2}}\text{,}$$		(2.1a)
	$$\bm{\phi}_{1}({{\bf x}})=\operatorname*{arg\,max}\lambda[\bm{\phi}({\bf{x}})]\text{.}$$		(2.1b)

The subsequent modes are defined to maximize the energy captured under the constraint that they are orthogonal to all previous ones,

$$\bm{\phi}_{j}({\bf x})=\operatorname*{arg\,max}_{\langle\bm{\phi}({{\bf x}}),\bm{\phi}_{k<j}({{\bf x}})\rangle_{{\bf x}}=0}\lambda[\bm{\phi}({\bf{x}})]\text{.}$$

(2.2)

The space-only inner product $\langle\cdot,\cdot\rangle_{{\bf x}}$, which defines the energy captured, is defined as an integral over the spatial domain $\Omega$,

$$\langle\bm{q}({{\bf x}}),\bm{\phi}({{\bf x}})\rangle_{{\bf x}}=\int_{\Omega}\bm{\phi}^{*}({{\bf x}})\mathsfbi{W}({\bf x})\bm{q}({\bf x})d{\bf x}\text{,}$$

(2.3)

where $\mathsfbi{W}({\bf x})$ is a weight matrix used to account for inter-variable importance or possibly to preference certain regions of the domain. One can show (Towne2018spectral; Lumley67) that the solution to the optimization problem (2.1b,2.2) is modes which are eigenfunctions of the space-only correlation tensor,

$$\int_{\Omega}\mathsfbi{C}({\bf x}_{1},{\bf x}_{2})\mathsfbi{W}({\bf x}_{2})\bm{\phi}_{j}({\bf x}_{2})d{\bf x}_{2}=\lambda_{j}\bm{\phi}_{j}({\bf x}_{1})\text{,}$$

(2.4)

where the eigenvalue is equal to the energy of the mode, i.e., $\lambda_{j}=\lambda[\bm{\phi}_{j}]$, and the correlation tensor is

$$\mathsfbi{C}({\bf x}_{1},{\bf x}_{2})=\mathbb{E}[\bm{q}({\bf x}_{1})\bm{q}^{*}({\bf x}_{2})]\text{.}$$

(2.5)

Whereas space-only POD modes optimally represent snapshots, space-time POD modes optimally represent trajectories over the time window $[0,T]$. The formulation is much the same as in space-only POD; the modes optimize the expected energy (2.1b,2.2), but in space-time POD, the inner product is over time as well as space,

$$\langle\bm{q}({{\bf x}},t),\bm{\phi}({{\bf x}},t)\rangle_{{\bf x},t}=\int_{0}^{T}\int_{\Omega}\bm{\phi}^{*}({{\bf x}},t)\mathsfbi{W}({\bf x})\bm{q}({\bf x},t)d{\bf x}\ dt\text{.}$$

(2.6)

The space-time POD modes that solve this optimization are eigenfunctions of the space-time correlation $\mathsfbi{C}({\bf x}_{1},t_{1},{\bf x}_{2},t_{2})=\mathbb{E}[\bm{q}({\bf x}_{1},t_{1})\bm{q}^{*}({\bf x}_{2},t_{2})]$, and the eigenvalues represent the energy (importance) of each mode. The most important property of space-time POD modes for the purpose of this paper is that the space-time POD reconstruction of a trajectory achieves lower error, on average, than the reconstruction with the same number of modes in any other space-time basis. More concretely, using the first $r$ space-time POD modes to reconstruct the trajectory,

$$\bm{q}^{r}({\bf x},t)=\sum_{j=1}^{r}\bm{\phi}_{j}({{\bf x}},t)\langle\bm{q}({\bf x},t),\bm{\phi}_{j}({{\bf x}},t)\rangle_{{\bf x},t}\text{,}$$

(2.7)

yields lower expected error

$$\mathbb{E}[\|\bm{q}^{r}({\bf x},t)-\bm{q}({\bf x},t)\|_{{\bf x},t}^{2}]$$

(2.8)

than would any other space-time basis. Above, expected error is measured over space and time using the norm $\|\cdot\|_{{\bf x},t}$ induced by the inner product.

Spectral POD is most easily understood as the frequency domain variant of space-only POD for statistically stationary systems. In other words, SPOD modes at a particular frequency optimally reconstruct (in the same sense as above) the state at that frequency, on average. The property that makes them attractive for model reduction, however, is that SPOD modes are also the long-time limit of space-time POD modes for statistically stationary systems. These ideas are made precise below, but for a more complete discussion, see (Towne2018spectral).

Spectral POD modes at frequency $k$ maximize

$$\lambda_{k}[\bm{\psi}({\bf{x}})]=\frac{\mathbb{E}\big{[}|\langle\hat{\bm{q}}_{k}({{\bf x}}),\bm{\psi}({{\bf x}})\rangle_{{\bf x}}|^{2}\big{]}}{\|\bm{\psi}({{\bf x}})\|^{2}}\text{,}$$

(2.9)

again, subject to the constraint that each mode $\bm{\psi}_{k,j}({\bf x})$ is orthogonal to the previous ones at that frequency $\bm{\psi}_{k,i<j}({\bf x})$. The Fourier-transformed state is defined as

$$\hat{\bm{q}}_{k}=\int_{-\infty}^{\infty}e^{-i\omega_{k}t}\bm{q}({\bf x},t)dt\text{.}$$

(2.10)

The solution to the optimization (2.9) is that the modes are eigenvectors of the cross-spectral density $\mathsfbi{S}_{k}({\bf x}_{1},{\bf x}_{2})=\mathbb{E}[\hat{\bm{q}}_{k}({\bf x}_{1})\hat{\bm{q}}_{k}^{*}({\bf x}_{2})]$

$$\int\limits_{\Omega}\mathsfbi{S}_{k}({{\bf x}}_{1},{\bm{x}}_{2})\mathsfbi{W}({{\bf x}}_{2})\bm{\psi}_{k,j}({{\bf x}}_{2})d{{\bf x}}_{2}=\lambda_{k,j}\bm{\psi}_{k,j}({{\bf x}}_{1})\text{.}$$

(2.11)

The eigenvalue is again equal to the energy of the mode, i.e., $\lambda_{k,j}=\lambda_{k}[\bm{\psi}_{k,j}]$. Modes at frequency $k$ have an implicit time dependence of $e^{i\omega_{k}t}$.

SPOD modes and their energies become identical to space-time POD modes as the time interval on which the latter are defined becomes long (Lumley67; Towne2018spectral; Frame23). Thus, the SPOD modes with the largest energies among all frequencies are the dominant space-time POD modes (for long times) and are most efficient for reconstructing long-time trajectories. We denote by $\tilde{\lambda}_{j}$ the $j$-th largest SPOD eigenvalue among all frequencies, which may be compared to the space-time POD eigenvalues. The convergence in the energy of the trajectory they capture is relatively fast, so for time intervals beyond a few correlation times, the SPOD modes capture nearly as much energy as the space-time modes (Frame23).

If the simulation time of a reduced-order model is long enough for this convergence to be met, the ability of the SPOD modes to capture structures is not diminished relative to that of space-time POD modes. SPOD modes also have two properties that make them more suitable for model reduction than space-time modes: they have analytic time dependence, and they are separable in space and time. The former makes some analytic progress possible in writing the equations that govern the modes and enables Fourier theory to be applied. The latter means that storing the modes requires $N_{t}$ times less memory, where $N_{t}$ is the number of times steps in the simulation.

Figure 1 shows the convergence in the representational ability of space-time POD and SPOD as the time interval becomes long. Specifically, to represent trajectories with some level of accuracy, $98\%$, say, one needs the same number of SPOD coefficients as space-time POD coefficients if the interval is long compared to the correlation time in the system. This convergence in representation ability occurs because the SPOD modes themselves converge to space-time POD modes in the limit of a long time interval. With space-only POD, one must specify the coefficients for every time step, which leads to a far less efficient encoding of the data because the coefficients are highly correlated from one time step to the next. That SPOD modes are near-optimal in representing trajectories, and that they are substantially more efficient than space-only POD modes motivate this work. If one can efficiently solve for some number of the SPOD coefficients of a trajectory, then these coefficients will lead to substantially lower error than solving for the same number of space-only POD coefficients.

Upon discretization, the continuous tensors become discrete, and the modes become $N_{x}$-component vectors in $\mathbb{C}^{N_{x}}$. Integration over space becomes matrix-vector multiplication. For example, the discrete SPOD modes at frequency $\omega_{k}$ are defined as the eigenvectors of the (weighted) cross-spectral density matrix,

$${\bf S}_{k}{\bf W}{\bf\Psi}_{k}={\bf\Psi}_{k}{\bf\Lambda}_{k}\text{.}$$

(2.12)

${\bf S}_{k}$ is the cross-spectral density, ${\bf\Psi}_{k}$ is the matrix with the discrete SPOD modes as its columns, and ${\bf\Lambda}_{k}$ is the diagonal matrix of eigenvalues which are the SPOD mode energies. It is important to note that in practice, the matrix ${\bf S}_{k}$ is not formed; the SPOD modes are calculated by the method of snapshots (Sirovich87; Towne2018spectral). In particular, given $r_{d}$ trajectories from which to obtain SPOD modes, each $N_{\omega}$ time steps in length, the discrete Fourier transform (DFT) of each trajectory is taken. This yields $r_{d}$ realizations of the $k$-th frequency, for every frequency $k$. These can be formed into a data matrix

$${\bf Q}_{k}=[\hat{\bm{q}}_{k}^{1},\hat{\bm{q}}_{k}^{2},\dots,\hat{\bm{q}}_{k}^{r_{d}}]\text{,}$$

(2.13)

where $\hat{\bm{q}}_{k}^{i}\in\mathbb{C}^{N_{x}}$ is the $k$-th frequency of the DFT of the $i$-th trajectory. The SPOD modes at frequency $\omega_{k}$ and the associated energies may then be obtained by first taking the singular value decomposition ${\bf U}{\bf\Sigma}{\bf V}^{*}=1/\sqrt{r_{d}}{\bf W}^{1/2}{\bf Q}_{k}$. The SPOD modes are then given by ${\bf W}^{-1/2}{\bf U}$ and the energies by ${\bf\Sigma}^{2}$ (Towne2018spectral).

A finite number of evenly spaced frequencies are retained. The lowest one, $\omega_{1}$, induces a time $T=2\pi/\omega_{1}$, which determines the interval $[0,T]$ on which the modes are periodic. The trajectories themselves are, of course, not periodic on this interval, so $T$ is the longest we may use SPOD modes for prediction, though, if a longer prediction is needed, the method may be repeated. Using the SPOD modes, the trajectory may be written

$$\bm{q}(t)=\frac{1}{N_{\omega}}\sum_{j=0}^{N_{\omega}-1}{\bf\Psi}_{j}\bm{a}_{j}e^{i\omega_{j}t}\text{,}$$

(2.14)

where $\bm{a}_{k}={\bf\Psi}_{k}^{*}{\bf W}\hat{\bm{q}}_{k}\in\mathbb{C}^{N_{x}}$ is the vector of expansion coefficients at the $k$-th frequency. The factor of $\frac{1}{N_{\omega}}$ makes (2.14) an (interpolated) inverse discrete Fourier transform with ${\bf\Psi}_{k}\bm{a}_{k}$ as the $k$-th coefficient in the DFT of $\bm{q}(t)$. The fully reduced representation of the trajectory is

$$\bm{q}^{r}(t)=\frac{1}{N_{\omega}}\sum_{j=0}^{N_{\omega}-1}{\bf\Psi}^{r}_{j}\bm{a}_{j}^{r}e^{i\omega_{j}t}\text{.}$$

(2.15)

In this paper, an $r$ superscript indicates a reduced quantity. The mean number of modes retained at each frequency is $r$, and $N_{\omega}r$ modes are retained in total. The same number is not retained at all frequencies because the most energetic space-time POD modes are the most energetic SPOD modes over all frequencies; the number of modes retained at the $k$-th frequency $r_{k}$ is the number of modes at this frequency that are among the $N_{\omega}r$ most energetic overall,

$$r_{k}=|\{l:\lambda_{k,l}\geq\tilde{\lambda}_{N_{\omega}r}\}|\text{.}$$

(2.16)

With these notations established, $\bm{a}_{k}^{r}\in\mathbb{C}^{r_{k}}$ and ${\bf\Psi}_{k}^{r}\in\mathbb{C}^{N_{x}\times r_{k}}$. Finally, the Fourier transformed trajectory is defined using the DFT as

$$\hat{\bm{q}}_{k}=\sum_{j=0}^{N_{\omega}-1}\bm{q}(j\Delta t)e^{-i\omega_{k}j\Delta t}\text{.}$$

(2.17)

The starting point is the following equation for $\bm{a}^{r}_{k}$ in terms of $\hat{\bm{q}}_{k}$:

$$\bm{a}^{r}_{k}={\bf\Psi}^{r*}_{k}{\bf W}\hat{\bm{q}}_{k}\text{.}$$

(3.2)

This equation is exact and follows from the fact that $\hat{\bm{q}}_{k}={\bf\Psi}_{k}\bm{a}_{k}$. The Fourier transformed state $\hat{\bm{q}}_{k}$ must be obtained from the known forcing and initial condition, and there is some subtlety here. For the sake of clarity in describing the model reduction approach, we ignore the subtlety in this section. The correct relation between the Fourier-transformed forcing and initial condition and the Fourier-transformed state is derived in Section 4 and is given in (4.14), and the model reduction method is adjusted for this correction in Section 5.

To write the equations in the frequency domain, we take the Fourier transform, replacing $\bm{q}(t)$ and $\bm{f}(t)$ with $\hat{\bm{q}}_{k}$ and $\hat{\bm{f}}_{k}$ and (naively, it turns out) replacing $\dot{\bm{q}}(t)$ with $i\omega_{k}\hat{\bm{q}}_{k}$, yielding

$$i\omega_{k}\hat{\bm{q}}_{k}={\bf A}\hat{\bm{q}}_{k}+{\bf B}\hat{\bm{f}}_{k}\text{.}$$

(3.3)

Then, assuming the stability of ${\bf A}$, the state is related to the forcing by

$$\hat{\bm{q}}_{k}={\bf R}_{k}{\bf B}\hat{\bm{f}}_{k}\text{,}$$

(3.4)

where the resolvent operator is defined as

$${\bf R}_{k}=(i\omega_{k}{\bf I}-{\bf A})^{-1}\text{.}$$

(3.5)

Plugging (3.4) into (3.2), we have

$$\bm{a}^{r}_{k}={\bf\Psi}^{r*}_{k}{\bf W}{\bf R}_{k}{\bf B}\hat{\bm{f}}_{k}\text{.}$$

(3.6)

This method may be derived as a Petrov-Galerkin method in the following way. Starting from the equations in the frequency domain, ${\bf L}_{k}\hat{\bm{q}}_{k}={\bf B}\hat{\bm{f}}_{k}$, where ${\bf L}_{k}=i\omega_{k}{\bf I}-{\bf A}$, and plugging in the representation for $\hat{\bm{q}}_{k}$, we have

$${\bf L}_{k}{\bf\Psi}_{k}^{r}\bm{a}_{k}^{r}={\bf B}\hat{\bm{f}}_{k}\text{.}$$

(3.7)

By left-multiplying the equations by ${\bf\Psi}_{k}^{r*}{\bf W}$ and inverting the matrices, the method would be Galerkin (Lin19; Towne21). Instead, (3.6) may be recovered by multiplying the equations on the left by the test basis ${\bf\Psi}_{k}^{r*}{\bf R}_{k}{\bf W}$, yielding

$${\bf\Psi}_{k}^{r*}{\bf W}{\bf R}_{k}{\bf L}_{k}{\bf\Psi}_{k}^{r}\bm{a}_{k}^{r}={\bf\Psi}_{k}^{r*}{\bf W}{\bf R}_{k}{\bf B}\hat{\bm{f}}_{k}\text{.}$$

(3.8)

The operator ${\bf\Psi}_{k}^{r*}{\bf W}{\bf R}_{k}{\bf L}_{k}{\bf\Psi}_{k}^{r}$ is the identity, so (3.6) is recovered and the method is Petrov-Galerkin.

We denote the matrix that maps the forcing at frequency $k$ to the SPOD mode coefficients at that frequency as ${\bf M}_{k}={\bf\Psi}^{r*}_{k}{\bf W}{\bf R}_{k}{\bf B}\in\mathbb{C}^{r_{k}\times N_{f}}$. For small problems (i.e., where $\mathcal{O}(N_{x}^{3})$ is not prohibitive) the resolvent, and thus ${\bf M}_{k}$, may be computed directly. More care is needed in larger problems, and, indeed, significant research has been devoted to approximating the resolvent operator in fluid mechanics (Farghadan23; Martini21), where it is often used in the context of stability theory (Trefethen93; Mckeon10; McKeon17; Schmidt18). These techniques may be used here to approximate ${\bf M}_{k}$ using some number of the dominant resolvent modes. Implementing these methods, however, may require substantial effort, and an alternative is possible due to the availability of data (the same data used to calculate the SPOD modes).

An accurate and simple-to-implement approximation of ${\bf M}_{k}$ can be obtained by leveraging the availability of data as follows. Defining

$$\bm{g}_{k}^{i}={\bf L}_{k}\hat{\bm{q}}_{k}^{i}\text{,}$$

(3.9)

where, again, $\hat{\bm{q}}_{k}^{i}$ is the $i$-th realization of $\hat{\bm{q}}_{k}$ in the training data, we have

$$\hat{\bm{q}}_{k}^{i}={\bf R}_{k}\bm{g}_{k}^{i}\text{.}$$

(3.10)

Using the many realizations of the training data (the same ones used to generate the SPOD modes), we have

$${\bf Q}_{k}={\bf R}_{k}{\bf G}_{k}\text{,}$$

(3.11)

where ${\bf Q}_{k}=[\hat{\bm{q}}_{k}^{1},\hat{\bm{q}}_{k}^{2},\dots]$ and ${\bf G}_{k}=[\bm{g}_{k}^{1},\bm{g}_{k}^{2},\dots]$. Multiplying this equation by its Hermitian transpose gives

$${\bf Q}_{k}{\bf Q}_{k}^{*}={\bf R}_{k}{\bf G}_{k}{\bf G}_{k}^{*}{\bf R}_{k}^{*}\text{.}$$

(3.12)

Both ${\bf Q}_{k}{\bf Q}_{k}^{*}$ and ${\bf G}_{k}{\bf G}_{k}^{*}$ may be represented by their eigendecompositions, i.e., their POD decompositions, giving

$${\bf\Psi}_{k}^{r}{\bf\Lambda}_{k}^{r}{\bf\Psi}_{k}^{r*}={\bf R}_{k}{\bf\Psi}_{k}^{gr}{\bf\Lambda}_{k}^{gr}{\bf\Psi}_{k}^{gr*}{\bf R}_{k}^{*}\text{,}$$

(3.13)

where ${\bf\Psi}_{k}^{gr}$ and ${\bf\Lambda}_{k}^{gr}$ come from the POD of ${\bf G}_{k}$ (the $r$ superscript is a reminder that ${\bf\Psi}_{k}^{gr}$ is not full-rank). By left-multiplying by ${\bf\Psi}^{r*}_{k}{\bf W}$ and right-multiplying by ${\bf L}^{*}_{k}{\bf\Psi}_{k}^{gr}{\bf\Lambda}_{k}^{gr-1}{\bf\Psi}_{k}^{gr*}{\bf B}$, we have

$${\bf\Lambda}_{k}^{r}({\bf L}_{k}{\bf\Psi}_{k}^{r})^{*}{\bf\Psi}^{gr}_{k}{\bf\Lambda}_{k}^{gr-1}{\bf\Psi}^{gr*}_{k}{\bf B}={\bf\Psi}_{k}^{r*}{\bf W}{\bf R}_{k}{\bf P}_{k}^{gr}{\bf B}\text{,}$$

(3.14)

where ${\bf P}_{k}^{g}={\bf\Psi}_{k}^{gr}{\bf\Psi}_{k}^{gr*}$ is the projection onto the modes of $\bm{g}_{k}$. The right-hand-side of (3.14) is an approximation of ${\bf M}_{k}$, which we label ${\bf E}_{k}$

$${\bf E}_{k}={\bf\Psi}_{k}^{r*}{\bf W}{\bf R}_{k}{\bf P}_{k}^{g}{\bf B}\approx{\bf M}_{k}\text{.}$$

(3.15)

If ${\bf P}_{k}^{g}$ were the identity, the approximation would be exact. The operator ${\bf M}_{k}$ is used to operate on forcing vectors, so the approximation is accurate to the extent that the forcing is within the span of the realizations of $\bm{g}_{k}$. The forcing is likely to be near this span because it is closely related to $\bm{g}$, which will become clear in Section 4.2.

The $k$-th Fourier coefficient of the time derivative $\hat{\dot{\bm{q}}}_{k}$ is not $i\omega_{k}\hat{\bm{q}}_{k}$. It is shown in appendix LABEL:App:Derivative_coefficient that in fact,

$$\hat{\dot{\bm{q}}}_{k}=i\omega_{k}\hat{\bm{q}}_{k}+\frac{\Delta\bm{q}}{T}\text{,}$$

(4.1)

where $\Delta\bm{q}=\bm{q}(T)-\bm{q}(0)$ is the change of the state on the time interval (Martini19; Nogueira20). If the state were periodic, this term would be zero, and the usual equation would hold, but the solution to the linear ODE has no reason to be periodic. How, then, does one solve in the frequency domain? Modifying (3.4) with $\frac{\Delta\bm{q}}{T}$ to give

$$\hat{\bm{q}}_{k}={\bf R}_{k}\left(\hat{\bm{f}}_{k}-\frac{\Delta\bm{q}}{T}\right)$$

(4.2)

is not useful as $\Delta\bm{q}$ is now known a priori (and, in some sense, is what we are solving for)!

To circumvent the problem described above, we use the analytic solution in time of (3.1). We then take the discrete Fourier transform (DFT) of this solution analytically to obtain $\hat{\bm{q}}$ in terms of $\hat{\bm{f}}$. We use the DFT rather than the Fourier series because, from the exact DFT, one can easily obtain the exact solution at a set of points in the time domain by taking the inverse DFT. On the other hand, with the first $N_{\omega}$ Fourier series coefficients, point values cannot be recovered exactly.

The analytic solution of (3.1) is (Grigoriu21)

$$\bm{q}(t)=e^{{\bf A}t}{\bf q}(0)+\int_{0}^{t}e^{{\bf A}(t-t^{\prime})}{\bf B}\bm{f}(t^{\prime})\ dt^{\prime}\text{.}$$

(4.3)

We want the (analytic) DFT of (4.3) where the $k$-th component of the DFT $\hat{\bm{q}}_{k}$ of $\bm{q}(t)$ is defined

$$\hat{\bm{q}}_{k}=\sum_{j=0}^{N_{\omega}-1}\bm{q}_{j}e^{-\frac{2\pi i}{N_{\omega}}jk}\text{,}$$

(4.4)

where $\bm{q}_{j}=\bm{q}(j\Delta t)$. Plugging (4.3) into (4.4), we have

$$\hat{\bm{q}}_{k}=\sum_{j=0}^{N_{\omega}-1}\left(e^{{\bf A}j\Delta t}\bm{q}_{0}+\int_{0}^{j\Delta t}e^{{\bf A}(j\Delta t-t^{\prime})}{\bf B}\bm{f}(t^{\prime})\ dt^{\prime}\right)e^{-\frac{2\pi i}{N_{\omega}}jk}\text{.}$$

(4.5)

For later convenience, we refer to the initial condition and forcing terms in (4.5) as $\hat{\bm{q}}_{k,ic}$ and $\hat{\bm{q}}_{k,force}$, respectively.

First, we evaluate $\hat{\bm{q}}_{k,ic}$. It may be rewritten as

$$\hat{\bm{q}}_{k,ic}=\sum_{j=0}^{N_{\omega}-1}e^{({\bf A}-i\omega_{k})j\Delta t}\bm{q}_{0}\text{,}$$

(4.6)

where we have defined $\Delta t=\frac{T}{N_{\omega}}$ and substituted $\omega_{k}=\frac{2\pi k}{T}$. The key to evaluating (4.6) is to notice that it is a matrix geometric sum, i.e., each term is the previous one multiplied by $e^{({\bf A}-i\omega_{k})\Delta t}$. The solution, entirely analogous to the scalar geometric sum, is

$$\hat{\bm{q}}_{k,ic}=({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t})^{-1}({\bf I}-e^{({\bf A}-i\omega_{k})N_{\omega}\Delta t})\bm{q}_{0}\text{.}$$

(4.7)

Because $e^{i\omega_{k}N_{\omega}\Delta t}=1$, this simplifies to

$$\hat{\bm{q}}_{k,ic}=({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t})^{-1}({\bf I}-e^{{\bf A}T})\bm{q}_{0}\text{.}$$

(4.8)

We view the first factor $({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t})^{-1}$ as a time-discrete resolvent operator; expanding the matrix exponential to first order in $\Delta t$, one recovers a multiple of the resolvent. Note, however, that truncating the expansion to first order is always invalid for the higher frequencies because, for these frequencies, $\omega_{k}\Delta t$ is order unity.

To evaluate $\hat{\bm{q}}_{k,force}$, we assume that $\bm{f}(t)$ can be written as

$$\bm{f}(t)=\frac{1}{N_{\omega}}\sum_{l=0}^{N_{\omega}-1}\hat{\bm{f}}_{l}e^{i\omega_{l}t}\text{.}$$

(4.9)

Though this is likely not exactly true in practice, it introduces minimal error (as we show in our numerical examples), and is necessary for evaluating the integral in (4.5). After pulling out the constant matrix exponential term and including the Fourier interpolation of $\bm{f}(t)$, the integral in (4.5) is

$$\hat{\bm{q}}_{k,force}=\frac{1}{N_{\omega}}\sum_{j=0}^{N_{\omega}-1}e^{-i\omega_{k}j\Delta t}e^{{\bf A}j\Delta t}\int_{0}^{j\Delta t}\sum_{l=0}^{N_{\omega}-1}e^{(i\omega_{l}-{\bf A})t^{\prime}}{\bf B}\hat{\bm{f}}_{l}\ dt^{\prime}\text{.}$$

(4.10)

Integrating (and hence inverting the matrix exponential within the integral) brings out a resolvent, and the integral evaluates to

$$\hat{\bm{q}}_{k,force}=\frac{1}{N_{\omega}}\sum_{j=0}^{N_{\omega}-1}e^{-i\omega_{k}j\Delta t}\sum_{l=0}^{N_{\omega}-1}{\bf R}_{l}\left(e^{i\omega_{l}j\Delta t}-e^{{\bf A}j\Delta t}\right){\bf B}\hat{\bm{f}}_{l}\text{,}$$

(4.11)

where, again, ${\bf R}_{l}=(i\omega_{l}{\bf I}-{\bf A})^{-1}$ is the resolvent operator at the $l$-th frequency. Now, using this result, the forcing term in (4.5) is

$$\hat{\bm{q}}_{k,force}=\frac{1}{N_{\omega}}\sum_{j=0}^{N_{\omega}-1}\sum_{l=0}^{N_{\omega}-1}{\bf R}_{l}\left(e^{i(\omega_{l}-\omega_{k})j\Delta t}-e^{({\bf A}-i\omega_{k}{\bf I})j\Delta t}\right){\bf B}\hat{\bm{f}}_{l}$$

(4.12)

The frequency difference term evaluates to zero for $\omega_{l}\neq\omega_{k}$, and the other term may be evaluated by the same geometric sum argument as before. The entire forcing term in (4.5) becomes

$$\hat{\bm{q}}_{k,force}={\bf R}_{k}{\bf B}\hat{\bm{f}}_{k}-\frac{1}{N_{\omega}}({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t})^{-1}({\bf I}-e^{{\bf A}T})\sum_{l=0}^{N_{\omega}-1}{\bf R}_{l}{\bf B}\hat{\bm{f}}_{l}\text{.}$$

(4.13)

Adding $\hat{\bm{q}}_{k,ic}$ and $\hat{\bm{q}}_{k,force}$, we obtain the following equation for the $k$-th component of the DFT of the state

$$\hat{\bm{q}}_{k}={\bf R}_{k}{\bf B}\hat{\bm{f}}_{k}+\left({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t}\right)^{-1}\left({\bf I}-e^{{\bf A}T}\right)\left(\bm{q}_{0}-\frac{1}{N_{\omega}}\sum_{l=0}^{N_{\omega}-1}{\bf R}_{l}{\bf B}\hat{\bm{f}}_{l}\right)\text{.}$$

(4.14)

The first term is the naive one. For a stable system, it may be shown that the correction term is inversely proportional to the length of the time interval. Therefore, while the error from ignoring it does decrease as the time interval increases in this case, it is not negligible unless the interval is extremely long. Intuitively, the correction term accounts for the transient in the solution, which decays exponentially at the rate prescribed by the slowest decaying eigenvector of ${\bf A}$. The decay time is (of course) independent of the time interval $T$ on which we take the transform, so the fraction of the time interval occupied by this decay is inversely proportional to $T$. The term $\frac{1}{N_{\omega}}\sum_{l=0}^{N_{\omega}-1}{\bf R}_{l}\hat{\bm{f}}_{l}$ is the initial condition that is ‘in sync’ with the forcing. If $\bm{q}_{0}$ were equal to this, then there would be no transient, $\bm{q}(t)$ would be exactly periodic on the interval, and the naive equation (3.4) would be correct. For long time intervals relative to the slowest decaying eigenvector of ${\bf A}$, the term $({\bf I}-e^{{\bf A}T})$ may be ignored. We also note that summing the operator multiplying the initial condition and forcing sum over all frequencies can be shown to give $N_{\omega}{\bf I}$

$$\sum_{k=0}^{N_{\omega}-1}\left({\bf I}-e^{({\bf A}-i\omega_{k})\Delta t}\right)^{-1}\left({\bf I}-e^{{\bf A}T}\right)=N_{\omega}{\bf I}\text{.}$$

(4.15)