time-periodic problems, parallel-in-time integration, projection-based PP-PC, Fourier-aware projection spaces

1 Introduction↩︎

Time-periodic problems arise naturally in computational models with periodic forcing or cyclic operation, where the quantity of interest is a periodic steady state rather than a transient trajectory. Examples include eddy-current simulations, cyclic chemical processes, and fluid–structure interaction models [1]–[4]. A typical finite-dimensional formulation seeks a function \(\mathbf{u}:\left[0,T\right]\to\mathbb{R}^d\) such that \[\label{eq:intro-periodic-problem} \mathbf{u}'\left(t\right) = \mathbf{f}\left(\mathbf{u}\left(t\right),t\right), \qquad t\in \left[0,T\right], \qquad \mathbf{u}\left(0\right)=\mathbf{u}\left(T\right),\tag{1}\] where \(\mathbf{f}:\mathbb{R}^d\times\mathbb{R}\to\mathbb{R}^d\) satisfies \(\mathbf{f}(\mathbf{x},t+T)=\mathbf{f}(\mathbf{x},t)\) for all \((\mathbf{x},t)\in\mathbb{R}^d\times\mathbb{R}\).

The periodic condition in 1 has been treated by several classes of numerical methods. Frequency-domain formulations exploit periodicity directly, from the classical Galerkin analysis of [5] to harmonic-balance methods for periodic flows [6]. Multigrid methods solve the periodicity-constrained space-time system for time-periodic parabolic problems [7], and waveform-relaxation methods provide an iterative alternative [8]. Shooting and Newton–Picard methods compute periodic states as fixed points of the period map [9], [10].

In large-scale time-domain simulations, spatial parallelism eventually saturates, leaving sequential time propagation as the bottleneck. This has motivated parallel-in-time algorithms that introduce concurrency across time subintervals. Representative parallel-in-time approaches include block boundary-value methods, parareal, PFASST, multigrid-reduction-in-time, and all-at-once time-domain preconditioning methods [11]–[15]; see [16] for a survey. The parareal algorithm realizes this idea through an iterative fine–coarse correction. At each iteration, the fine propagator is applied independently on the time subintervals, and the coarse propagator updates the solution values at the coarse time points. The original parareal method was introduced by Lions, Maday, and Turinici for initial-value problems [17]; see also [18] for its convergence analysis. For time-periodic problems, the periodic condition must be incorporated into the parareal iteration. Gander, Jiang, Song, and Zhang introduced and analyzed a periodic variant [19]: the periodic parareal algorithm with a periodic coarse problem (PP-PC). PP-PC preserves periodicity on the coarse time points at every iteration. However, the PP-PC iteration may converge slowly for oscillatory or wave-like dynamics.

Projection-based correction has been widely used to improve the convergence of parallel-in-time methods. An early time-decomposed framework for implicit parallel-in-time integration was introduced in [20]. Projection has also been used to stabilize parareal iterations for first- and second-order hyperbolic systems [21]. For initial-value problems, the Krylov-enhanced parareal method of [22] improves the iteration through projection. The Krylov-enhanced PP-PC method of [23] extends the Krylov subspace correction to PP-PC, the periodic form of parareal, for time-periodic problems. In these projection-based parareal methods, the projection space determines how fine and coarse propagation are combined, so its choice affects the convergence behavior of the iteration.

In this paper, we introduce a Fourier-aware projection-based PP-PC method, which combines a new correction scheme with a new construction of the projection space. The scheme is discrepancy-based and differs from that of Krylov-enhanced PP-PC. For the projection space, Krylov-enhanced PP-PC uses the full solution-snapshot space, built from the computed solution values over one period. However, Fourier techniques are widely used to exploit temporal structure in time-dependent and time-periodic solvers [14], [24]–[26]. We develop a Fourier-aware construction of the projection space that exploits this structure. The construction represents coarse-time histories in the Fourier basis associated with the coarse-time index. A mode-selection rule specifies which temporal modes to retain. The corresponding modal components are then combined into a mixed Fourier space, which serves as the projection space.

We analyze the convergence of projection-based PP-PC with the discrepancy-based correction scheme for general nonlinear time-periodic problems. The starting point is a local one-step error estimate that holds for a discrepancy-form propagator and an arbitrary projection space, under local smoothness, nondegeneracy, and smallness assumptions. The estimate shows that the error after one iteration is bounded in terms of the unresolved error and two explicit nonlinear contributions. When this unresolved part is represented in temporal Fourier modes, it is bounded by a tail from modes not selected together with a leak from selected modes not fully captured by the projection space. The tail and the leak are both small when the selected set contains the error’s dominant temporal modes and the projection space captures their modal coefficients. Our main result (1) is a Fourier-aware one-step estimate in which the tail–leak form enters the unresolved-error term and the nonlinear remainder, while the base-state term remains unchanged. For linear time-periodic problems, the nonlinear terms vanish; under weaker assumptions, the estimates then reduce exactly to global bounds with the same constants: a general-projection estimate and a tail–leak estimate.

The rest of this paper is organized as follows. 2 recalls PP-PC and presents its projection-based framework. 3 introduces Fourier-aware projection-based PP-PC, its discrepancy-based correction scheme, and its mixed Fourier space. 4 develops the local nonlinear estimate, the Fourier-aware results, and the linear corollaries. 5 reports experiments on two linear and two nonlinear time-periodic problems, and 6 concludes the paper.

2 PP-PC and a framework for its projection-based variants↩︎

In this section we first recall PP-PC and present projection-based PP-PC as a general framework. Then we show that Krylov-enhanced PP-PC is an instance of this framework; our own instance is developed in 3.

2.1 Periodic parareal with a periodic coarse problem↩︎

We begin by recalling the original PP-PC method. Let \[0=T_0<T_1<\cdots<T_N=T, \qquad T_n=n\Delta T,\qquad \Delta T=T/N .\] The points \(T_n\) and the intervals \(\left[T_n,T_{n+1}\right]\) are called the coarse time points and coarse time intervals, respectively. The solution values at the coarse time points are called the interface values. They are vectors in a real finite-dimensional state space, denoted by \(\mathbb{R}^d\). At iteration \(k\) the interface values are \[\mathbf{U}^k = \left(\mathbf{U}_0^k,\ldots,\mathbf{U}_{N-1}^k\right) \in\left(\mathbb{R}^d\right)^N .\] The periodic condition identifies the endpoint value with the initial value: \[\mathbf{U}_N^k=\mathbf{U}_0^k .\]

The PP-PC method is built from a fine propagator and a coarse propagator on each coarse time interval. For \(0\le n<N\) and \(\mathbf{x}\in\mathbb{R}^d\), let \[\mathcal{F}\left(T_{n+1},T_n,\mathbf{x}\right), \qquad \mathcal{G}\left(T_{n+1},T_n,\mathbf{x}\right)\] denote the corresponding fine and coarse propagator values over \(\left[T_n,T_{n+1}\right]\), respectively, starting from the state \(\mathbf{x}\) at \(T_n\). More precisely, if \(\mathbf{u}\left(\cdot;T_n,\mathbf{x}\right)\) denotes the solution of the underlying evolution problem satisfying \(\mathbf{u}\left(T_n;T_n,\mathbf{x}\right)=\mathbf{x}\), then \(\mathcal{F}\left(T_{n+1},T_n,\mathbf{x}\right)\) gives an accurate approximation of \(\mathbf{u}\left(T_{n+1};T_n,\mathbf{x}\right)\). The coarse propagator \(\mathcal{G}\left(T_{n+1},T_n,\mathbf{x}\right)\) gives a less accurate but less expensive approximation of the same state. When no ambiguity arises, we write \[\mathcal{F}_n\left(\mathbf{x}\right):=\mathcal{F}\left(T_{n+1},T_n,\mathbf{x}\right), \qquad \mathcal{G}_n\left(\mathbf{x}\right):=\mathcal{G}\left(T_{n+1},T_n,\mathbf{x}\right).\]

The PP-PC iteration [19] computes \(\mathbf{U}^{k+1}\) from \(\mathbf{U}^k\) by \[\label{eq:pppc-expanded} \begin{align} \mathbf{U}_0^{k+1} &= \mathcal{F}_{N-1}\left(\mathbf{U}_{N-1}^k\right) + \mathcal{G}_{N-1}\left(\mathbf{U}_{N-1}^{k+1}\right) - \mathcal{G}_{N-1}\left(\mathbf{U}_{N-1}^k\right),\\ \mathbf{U}_{n+1}^{k+1} &= \mathcal{F}_n\left(\mathbf{U}_n^k\right) + \mathcal{G}_n\left(\mathbf{U}_n^{k+1}\right) - \mathcal{G}_n\left(\mathbf{U}_n^k\right), \qquad n=0,\ldots,N-2 . \end{align}\tag{2}\] Once \(\mathbf{U}^k\) is known, the fine and coarse propagator evaluations at \(\mathbf{U}_n^k\) are available. These fine propagator evaluations are independent across the coarse time intervals and can therefore be performed in parallel. The new values \(\mathbf{U}_n^{k+1}\) remain coupled through the terms \(\mathcal{G}_n\left(\mathbf{U}_n^{k+1}\right)\). The first equation closes the cycle by coupling \(\mathbf{U}_0^{k+1}\) to \(\mathbf{U}_{N-1}^{k+1}\). This coupled solve is the periodic coarse problem in PP-PC.

2.2 A framework for projection-based PP-PC↩︎

In PP-PC the correction in 2 is carried by the coarse propagator \(\mathcal{G}_n\). Projection-based PP-PC generalizes this step. At iteration \(k\) it constructs a projection space \(\mathcal{V}^k\subset\mathbb{R}^d\), with orthogonal projection \(P^k\) onto it, and replaces \(\mathcal{G}_n\) in the correction by a projection-based propagator \(\mathcal{X}_{P^k,n}\). It computes \(\mathbf{U}^{k+1}\) from the coupled periodic problem \[\label{eq:projection-corrected-pppc} \begin{align} \mathbf{U}_0^{k+1} &= \mathcal{F}_{N-1}\left(\mathbf{U}_{N-1}^k\right) + \mathcal{X}_{P^k,N-1}\left(\mathbf{U}_{N-1}^{k+1}\right) - \mathcal{X}_{P^k,N-1}\left(\mathbf{U}_{N-1}^k\right),\\ \mathbf{U}_{n+1}^{k+1} &= \mathcal{F}_n\left(\mathbf{U}_n^k\right) + \mathcal{X}_{P^k,n}\left(\mathbf{U}_n^{k+1}\right) - \mathcal{X}_{P^k,n}\left(\mathbf{U}_n^k\right), \qquad n=0,\ldots,N-2 . \end{align}\tag{3}\] Taking \(\mathcal{X}_{P,n}=\mathcal{G}_n\) recovers PP-PC 2 . The framework leaves two ingredients to be specified: the projection space \(\mathcal{V}^k\) and the projection-based propagator \(\mathcal{X}_{P^k,n}\). Different choices give different methods, while the structure of the iteration is unchanged. Once \(P^k\) and \(\mathcal{X}_{P^k,n}\) are fixed, the values \(\mathcal{F}_n(\mathbf{U}_n^k)\) and \(\mathcal{X}_{P^k,n}(\mathbf{U}_n^k)\) at the known iterate are independent across the coarse intervals and are computed in parallel; the new values \(\mathbf{U}_n^{k+1}\) remain coupled only through \(\mathcal{X}_{P^k,n}(\mathbf{U}_n^{k+1})\), and the first equation closes the cycle. 1 summarizes the framework.

Krylov-enhanced PP-PC [23] is an instance of this framework. For initial-value problems, the Krylov-enhanced parareal method [22] accelerates parareal by a projection built from solution snapshots; the periodic version applies the same idea to PP-PC. Its projection space is the full solution-snapshot space \[\label{eq:krylov-snapshot-space} \mathcal{V}_{\rm K}^k := \mathop{\mathrm{span}}\left\{\,\mathbf{U}_n^j:\;0\le j\le k,\;0\le n<N\,\right\},\tag{4}\] which in the present periodic notation collects all interface values available up to iteration \(k\). Let \(P_{\rm K}^k\) be the orthogonal projection onto \(\mathcal{V}_{\rm K}^k\) and set \(Q_{\rm K}^k:=I-P_{\rm K}^k\), where \(I\) is the identity matrix on \(\mathbb{R}^d\). Its projection-based propagator is the split form \[\label{eq:krylov-enhanced-propagator} \mathcal{K}_{{\rm K},n}^k\left(\mathbf{x}\right) := \mathcal{F}_n\left(P_{\rm K}^k\mathbf{x}\right) + \mathcal{G}_n\left(Q_{\rm K}^k\mathbf{x}\right),\tag{5}\] which applies the fine propagator on the snapshot space and the coarse propagator on its complement. Krylov-enhanced PP-PC is therefore 1 with the projection space of step [line:gen-proj] set to the snapshot space \(\mathcal{V}_{\rm K}^k\) and the projection-based propagator of step [line:gen-solve] to the split form \(\mathcal{K}_{{\rm K},n}^k\); equivalently, the iteration 3 with \(\mathcal{X}_{P^k,n}=\mathcal{K}_{{\rm K},n}^k\).

3 Fourier-aware projection-based PP-PC↩︎

Under the framework of 2.2, our method makes two new choices: a discrepancy-based correction scheme and a Fourier-aware projection space. This correction scheme is realized by the discrepancy-form propagator, the precise object we define first.

The propagator \(\mathcal{H}_{P,n}\) is not a coarse propagator. It adds to the coarse value \(\mathcal{G}_n(\mathbf{x})\) the fine–coarse discrepancy evaluated at the projected state \(P\mathbf{x}\); projection-based PP-PC with this propagator is 3 with \(\mathcal{X}_{P^k,n}=\mathcal{H}_{P^k,n}\).

The remaining choice is the projection space, the central algorithmic issue in projection-based PP-PC. Because the discrepancy-form propagator applies fine–coarse information only to the projected component, a projection space is useful to the extent that it captures the components that the iteration still has to correct. For time-periodic problems these components have temporal structure that the discrete Fourier transform (DFT) in the coarse-time index makes explicit.

At iteration \(k\), the interface values \[\mathbf{U}^k = \left(\mathbf{U}_0^k,\ldots,\mathbf{U}_{N-1}^k\right)\] form the solution history over one period. Its temporal Fourier coefficients are \[\label{eq:temporal-fourier-coefficients} \widehat{\mathbf{U}}_\ell^k = \frac{1}{\sqrt N} \sum_{n=0}^{N-1} \mathbf{U}_n^k \exp\left(-2\pi \imath\ell n/N\right), \qquad \ell=0,\ldots,N-1.\tag{6}\] The same transform applies to any coarse-time history. Here \(\imath\) is the imaginary unit. The index \(\ell\) represents a temporal Fourier mode, with phase increment \(2\pi\ell/N\) on the coarse time grid. For low positive indices the corresponding physical angular frequency is \(2\pi\ell/T\); indices near \(N\) represent the negative frequencies \(2\pi(\ell-N)/T\). The coefficient \(\widehat{\mathbf{U}}_\ell^k\) is the content of the solution history in that mode. For notational simplicity we identify a mode with its index \(\ell\), so that a set of temporal Fourier modes is a subset of \(\{0,\ldots,N-1\}\).

The mode \(\ell=0\) is the mean component over the period. Modes with small positive indices describe slowly varying components on the cyclic coarse time grid. Thus the coefficients \(\widehat{\mathbf{U}}_\ell^k\) describe how the solution history is distributed among temporal frequencies of the periodic iteration. Truncated Fourier representations have proved effective in several time-periodic applications [1], [6]. The construction exploits temporal concentration when it is present: a mode-selection rule first chooses the temporally relevant modes (3.1), and the projection space is then built from the modal components of the solution and the fine–coarse discrepancy histories at the selected modes (3.2).

3.1 Selection of temporal modes↩︎

The Fourier-aware construction first chooses a selected temporal mode set. At iteration \(k\), the selected temporal mode set is denoted by \[\mathcal{I}^k\subset\left\{0,\ldots,N-1\right\}.\] When only a few temporal modes carry significant content, the construction aims to retain those modes. Convergence may degrade if the selected set omits a mode carrying a significant part of the current error. The two choices below correspond to two practical regimes:

• The fixed choice applies when the dominant temporal frequencies of the response are known in advance. It then uses the same prescribed mode set at every iteration, \[\mathcal{I}^k=\mathcal{I}_{\mathrm{fix}}\quad\text{for all } k, \qquad \mathcal{I}_{\mathrm{fix}}\subset\left\{0,\ldots,N-1\right\},\] with \(\mathcal{I}_{\mathrm{fix}}\) the modes at those frequencies. It is appropriate when the periodic response is dominated by low temporal frequencies, for which we take \(\mathcal{I}_{\mathrm{fix}}\) to be the symmetric low-frequency band \[\label{eq:low-frequency-band} \mathcal{I}_L = \left\{0,1,\ldots,L\right\}\cup\left\{N-L,\ldots,N-1\right\}, \qquad 1\le L<N/2 .\tag{7}\] This band retains the mean mode and the first \(L\) positive temporal frequencies together with their negative-frequency counterparts.

• The adaptive increment-energy choice applies when the dominant temporal frequencies are not known in advance, or change during the iteration. It uses the current iteration increment. For \(k\ge 1\), set \[\Delta \mathbf{U}_n^k = \mathbf{U}_n^k-\mathbf{U}_n^{k-1}.\] Let \(\widehat{\Delta\mathbf{U}}_\ell^k\) be its temporal Fourier coefficients and define \[\label{eq:increment-energy} a_\ell^k = \left\|\widehat{\Delta\mathbf{U}}_\ell^k\right\|_2^2, \qquad \ell=0,\ldots,N-1.\tag{8}\] Since the histories are real, \(a_\ell^k=a_{(N-\ell)\bmod N}^k\), so the energies are properties of the conjugate groups \(\{\ell,(N-\ell)\bmod N\}\). The adaptive rule ranks these groups by their combined energy and retains the \(M\) highest-ranked groups, where \(M\) is a prescribed number. The mode \(\ell=0\) is self-conjugate, as is \(\ell=N/2\) when \(N\) is even, so the selected set \(\mathcal{I}^k\) contains at most \(2M\) modes. For \(k=0\) no increment is available; we set \(\mathcal{I}^0=\emptyset\), so the first update reduces to the PP-PC step. The increment estimates the iteration error: for a contracting iteration, \(\Delta\mathbf{U}^k\) is dominated by the error still to be removed, so its leading temporal modes follow those of the error. This link is heuristic: for a noncontracting or oscillatory iteration the temporal spectrum of \(\Delta\mathbf{U}^k\) need not match that of the error, and the adequacy of the selection is assessed a posteriori in 5.

3.2 Construction of the Fourier-aware projection space↩︎

A selected mode affects the correction only through its modal components in the projection space. The construction of the projection space is guided by the role of \(P\) in the projection-based PP-PC update 3 and in the discrepancy-form propagator ?? . For an orthogonal projection \(P\) and the current solution history \(\mathbf{U}^k\), define the correction source \[\mathbf{c}_{n+1}^k\left(P\right) := \mathcal{F}_n\left(\mathbf{U}_n^k\right) - \mathcal{H}_{P,n}\left(\mathbf{U}_n^k\right), \qquad n=0,\ldots,N-1 .\] By the definition of \(\mathcal{H}_{P,n}\), this source is \[\mathbf{c}_{n+1}^k\left(P\right) = \mathcal{D}_n\left(\mathbf{U}_n^k\right) - \mathcal{D}_n\left(P\mathbf{U}_n^k\right).\] With the discrepancy-form propagator, the projection-based PP-PC step is the cyclic system—its operator matrix understood rowwise, as made precise below— \[\label{eq:cyclic-system} \begin{bmatrix} I& 0 & \cdots & -\mathcal{H}_{P,N-1}\left(\cdot\right)\\ -\mathcal{H}_{P,0}\left(\cdot\right) & I& & 0\\ \vdots & \ddots & \ddots & \vdots\\ 0 & \cdots & -\mathcal{H}_{P,N-2}\left(\cdot\right) & I \end{bmatrix} \begin{bmatrix} \mathbf{U}_0^{k+1}\\ \mathbf{U}_1^{k+1}\\ \vdots\\ \mathbf{U}_{N-1}^{k+1} \end{bmatrix} = \begin{bmatrix} \mathbf{c}_N^k\left(P\right)\\ \mathbf{c}_1^k\left(P\right)\\ \vdots\\ \mathbf{c}_{N-1}^k\left(P\right) \end{bmatrix}.\tag{9}\] We refer to the operator on the left-hand side as the cyclic operator associated with \(\mathcal{H}_{P,n}\), and denote it by \(\mathcal{C}_P\). Throughout, indices of histories are read modulo \(N\), so that \(\mathcal{C}_P\) acts rowwise as \(\left(\mathcal{C}_P\mathbf{w}\right)_{n+1}=\mathbf{w}_{n+1}-\mathcal{H}_{P,n}(\mathbf{w}_n)\), \(n=0,\ldots,N-1\). The projection acts at the states where the discrepancy-form propagator is evaluated and in the correction source \(\mathbf{c}_{n+1}^k\left(P\right)\); the projection space should therefore be informed by both the solution and the fine–coarse discrepancy histories:

• The solution histories contain those evaluation states \(\mathbf{U}_n^m\); their selected modal components identify state-space directions associated with the retained temporal content of the periodic iterates.

• The discrepancy histories contain the correction directions \(\mathbf{d}_n^m=\mathcal{D}_n\left(\mathbf{U}_n^m\right)\); since \(\mathbf{c}_{n+1}^k\left(P\right)\) compares \(\mathcal{D}_n\left(\mathbf{U}_n^k\right)\) with \(\mathcal{D}_n\left(P\mathbf{U}_n^k\right)\), they identify the discrepancy components the projection should capture.

The projection space is therefore built from both histories. The one-step estimate in 4 contains a projection-dependent source contribution obtained by applying a local fine–coarse discrepancy bound to the unresolved part of the current error. A cyclic stability factor also enters the estimate, together with two additional contributions in the nonlinear case.

For each \(0\le m\le k\) and \(\ell=0,\ldots,N-1\), the solution history \(\mathbf{U}^m\) and the fine–coarse discrepancy history \(\mathbf{d}^m\) have temporal Fourier coefficients \[\label{eq:solution-discrepancy-dft} \widehat{\mathbf{U}}_\ell^m = \frac{1}{\sqrt N}\sum_{n=0}^{N-1}\mathbf{U}_n^m \exp\left(-2\pi\imath\ell n/N\right), \quad \widehat{\mathbf{d}}_\ell^m = \frac{1}{\sqrt N}\sum_{n=0}^{N-1}\mathbf{d}_n^m \exp\left(-2\pi\imath\ell n/N\right).\tag{10}\] These coefficients are complex in general, whereas the projection space is a subspace of \(\mathbb{R}^d\); we therefore take their real and imaginary parts, omitting any that vanish. The selected modal components of the solution history span the solution-generated Fourier space \[\label{eq:solution-generated-fourier-space} \mathcal{V}_{\rm sol}^k\left(\mathcal{I}^k\right) := \mathop{\mathrm{span}}\left\{ \operatorname{Re}\widehat{\mathbf{U}}_\ell^m, \operatorname{Im}\widehat{\mathbf{U}}_\ell^m: \ell\in\mathcal{I}^k,\;0\le m\le k \right\}\tag{11}\] and those of the fine–coarse discrepancy history span the discrepancy-generated Fourier space \[\label{eq:discrepancy-generated-fourier-space} \mathcal{V}_{\rm disc}^k\left(\mathcal{I}^k\right) := \mathop{\mathrm{span}}\left\{ \operatorname{Re}\widehat{\mathbf{d}}_\ell^m, \operatorname{Im}\widehat{\mathbf{d}}_\ell^m: \ell\in\mathcal{I}^k,\;0\le m\le k \right\}.\tag{12}\] The Fourier-aware projection space is their sum, the mixed Fourier space \[\label{eq:mixed-fourier-projection-space} \mathcal{V}_{\rm mix}^k\left(\mathcal{I}^k\right) := \mathcal{V}_{\rm sol}^k\left(\mathcal{I}^k\right) + \mathcal{V}_{\rm disc}^k\left(\mathcal{I}^k\right).\tag{13}\] It is the smallest subspace containing both, so a single projection \(P\) onto \(\mathcal{V}_{\rm mix}^k\left(\mathcal{I}^k\right)\) captures the selected modal components of the evaluation states and of the discrepancy directions entering the correction source.

3.3 The Fourier-aware projection-based PP-PC algorithm↩︎

We now combine these ingredients into Fourier-aware projection-based PP-PC. At each iteration \(k\), choose the selected temporal mode set \(\mathcal{I}^k\) by one of the rules in 3.1. The projection space is the mixed Fourier space \(\mathcal{V}_{\rm mix}^k\left(\mathcal{I}^k\right)\) of 13 , and \(P_{\mathcal{I}^k}\) is the orthogonal projection onto it. The Fourier-aware projection-based PP-PC update is the periodic solve 3 with the discrepancy-form propagator and \(P^k=P_{\mathcal{I}^k}\), namely \[\label{eq:fourier-aware-pppc} \begin{align} \mathbf{U}_0^{k+1} &= \mathcal{F}_{N-1}\left(\mathbf{U}_{N-1}^k\right) + \mathcal{H}_{P_{\mathcal{I}^k},N-1}\left(\mathbf{U}_{N-1}^{k+1}\right) - \mathcal{H}_{P_{\mathcal{I}^k},N-1}\left(\mathbf{U}_{N-1}^k\right),\\ \mathbf{U}_{n+1}^{k+1} &= \mathcal{F}_n\left(\mathbf{U}_n^k\right) + \mathcal{H}_{P_{\mathcal{I}^k},n}\left(\mathbf{U}_n^{k+1}\right) - \mathcal{H}_{P_{\mathcal{I}^k},n}\left(\mathbf{U}_n^k\right), \qquad n=0,\ldots,N-2 . \end{align}\tag{14}\] 2 summarizes the procedure. Its main parameter is the band width \(L\) of the fixed rule or the mode count \(M\) of the adaptive rule; the experimental values are given in 5. In an implementation, past Fourier coefficients are stored, so only those of the newly generated solution and fine–coarse discrepancy histories are computed at iteration \(k\).

The discrepancy-form propagator \(\mathcal{H}_{P,n}(\mathbf{x})=\mathcal{G}_n(\mathbf{x})+\mathcal{D}_n(P\mathbf{x})\) evaluates the fine propagator at the projected state \(P\mathbf{x}\). In the corrected update 14 this evaluation falls on the next iterate \(\mathbf{U}^{k+1}\), the unknown of the cyclic system 9 , so the cost of the correction depends on how that system is solved. For affine \(\mathcal{F}_n\) and \(\mathcal{G}_n\) the cyclic system is linear. Its operator is assembled from the fine propagator applied to the projection basis, and the correction is a single linear solve with no further fine propagations; the fine propagator is otherwise evaluated only at the current iterate \(\mathbf{U}^k\), in parallel across the coarse intervals as in the original PP-PC.

4 Error and convergence analysis↩︎

This section analyzes the error and convergence of projection-based PP-PC with our discrepancy-based correction scheme, realized by the discrepancy-form propagator ?? . The analysis considers the general nonlinear time-periodic problem \[\label{eq:nonlinear-periodic-problem} \mathbf{u}'(t)=\mathbf{f}(\mathbf{u}(t),t), \qquad t\in[0,T], \qquad \mathbf{u}(0)=\mathbf{u}(T),\tag{15}\] where \(\mathbf{f}(\mathbf{x},t+T)=\mathbf{f}(\mathbf{x},t)\) for all \(\mathbf{x}\) in a domain containing the states considered below.

Throughout, \(\mathrm D\) denotes the Jacobian of a map with respect to the state, the subscript \(\mathbf{u}\) marking the partial Jacobian of \(\mathbf{f}(\mathbf{u},t)\) in its state argument; \(\mathbf{f}\) and \(\mathrm D_{\mathbf{u}}\mathbf{f}\) are assumed continuous in \((\mathbf{u},t)\). Let \(\mathbf{u}^\star_{\mathrm c}\) denote a \(T\)-periodic solution, assumed to exist; the subscript marks the continuous-time solution. Linearizing 15 about \(\mathbf{u}^\star_{\mathrm c}\) gives the variational equation \[\label{eq:floquet-variational} \dot{\mathbf{v}}=A_\star(t)\,\mathbf{v}, \qquad A_\star(t):=\mathrm D_{\mathbf{u}}\mathbf{f}\!\left(\mathbf{u}^\star_{\mathrm c}(t),t\right),\tag{16}\] a linear system with \(T\)-periodic coefficients. Let its fundamental matrix \(\Phi\) satisfy \(\dot{\Phi}=A_\star\Phi\) and \(\Phi(0)=I\). Then \(\mathbf{v}(T)=\Phi(T)\mathbf{v}(0)\), where \(\Phi(T)\) is the monodromy matrix, whose eigenvalues are the Floquet multipliers [10], [27], [28].

Define the period map \(\varphi_T(\mathbf{x}):=\mathbf{u}(T;\mathbf{x})\), where \(\mathbf{u}(0;\mathbf{x})=\mathbf{x}\). Since \(\mathrm D\varphi_T(\mathbf{u}^\star_{\mathrm c}(0))=\Phi(T)\) and \(I-\Phi(T)\) is invertible by 1, the implicit function theorem implies that \(\mathbf{u}^\star_{\mathrm c}(0)\) is locally unique among the fixed points of \(\varphi_T\). If the fine period map \(\Phi_{\mathcal{F}}:=\mathcal{F}_{N-1}\circ\cdots\circ\mathcal{F}_0\) converges to \(\varphi_T\) in \(C^1\) near \(\mathbf{u}^\star_{\mathrm c}(0)\), then for sufficiently small fine steps \(\Phi_{\mathcal{F}}\) has a locally unique fixed point \(\mathbf{u}_0^\star\) near \(\mathbf{u}^\star_{\mathrm c}(0)\). This fixed point is nondegenerate: \(I-\mathrm D\Phi_{\mathcal{F}}(\mathbf{u}_0^\star)\) is invertible. We thus define the fine periodic reference \[\label{eq:fine-periodic-reference} \mathbf{u}^\star = \left(\mathbf{u}_0^\star,\ldots,\mathbf{u}_{N-1}^\star\right) \in \left(\mathbb{R}^d\right)^N,\tag{17}\] by the fine-propagator interface equations \[\label{eq:fine-periodic-reference-interface} \mathbf{u}_{n+1}^\star=\mathcal{F}_n\left(\mathbf{u}_n^\star\right), \qquad n=0,\ldots,N-1, \qquad \mathbf{u}_N^\star=\mathbf{u}_0^\star .\tag{18}\] The preceding argument gives a sufficient condition for the existence of this discrete reference. The estimates below assume its existence and measure the iteration error of \(\mathbf{U}^k\) relative to it. The error at the coarse time points is \[\mathbf{e}_n^k:=\mathbf{u}_n^\star-\mathbf{U}_n^k, \qquad n=0,\ldots,N-1,\] its maximum-in-time size is \[E_\infty^k:=\max_{0\le n<N}\left\lVert \mathbf{e}_n^k \right\rVert_2 ,\] and the error history is \[\mathbf{e}^k:=\left(\mathbf{e}_0^k,\ldots,\mathbf{e}_{N-1}^k\right) \in\left(\mathbb{R}^d\right)^N .\] On histories \(\mathbf{w}=(\mathbf{w}_0,\ldots,\mathbf{w}_{N-1})\in(\mathbb{R}^d)^N\) we use the maximum-in-time norm and the operator norm it induces, \[\left\lVert \mathbf{w} \right\rVert_\infty:=\max_{0\le n<N}\left\lVert \mathbf{w}_n \right\rVert_2 , \qquad \left\lVert S \right\rVert_{\infty\to\infty}:=\sup_{\mathbf{w}\neq\mathbf{0}}\frac{\left\lVert S\mathbf{w} \right\rVert_\infty}{\left\lVert \mathbf{w} \right\rVert_\infty} ,\] so that \(E_\infty^k=\left\lVert \mathbf{e}^k \right\rVert_\infty\). History indices are taken modulo \(N\); in particular, \(\mathbf{U}_N^{k+1}=\mathbf{U}_0^{k+1}\), \(\mathbf{u}_N^\star=\mathbf{u}_0^\star\), and \(\mathbf{e}_N^{k+1}=\mathbf{e}_0^{k+1}\). For one projection-based PP-PC update, fix an arbitrary projection space \(\mathcal{V}\subset\mathbb{R}^d\), let \(P\) be its orthogonal projection, and set \(Q:=I-P\). We call \(Q\mathbf{w}=(Q\mathbf{w}_0,\ldots,Q\mathbf{w}_{N-1})\in(\mathbb{R}^d)^N\) the unresolved part of \(\mathbf{w}\). The local one-step estimate uses three assumptions.

The following additional assumption is used only in the contractive special case discussed in 5.

4.1 Projection-dependent one-step analysis↩︎

For one update of 3 using the discrepancy-form propagator with \(P^k=P\), define the unresolved fine–coarse discrepancy \[\label{eq:nonlinear-unresolved-correction} \mathcal{R}_{P,n}(\mathbf{x}) := \mathcal{F}_n(\mathbf{x})-\mathcal{H}_{P,n}(\mathbf{x}) = \mathcal{D}_n(\mathbf{x})-\mathcal{D}_n(P\mathbf{x}).\tag{19}\] If \(E_\infty^k>0\), define the unresolved-error ratio \(\theta_Q^k(P)\) and the projected-error ratio \(\theta_P^k(P)\) by \[\theta_Q^k(P) := \frac{\max_{0\le n<N}\left\lVert Q\mathbf{e}_n^k \right\rVert_2}{E_\infty^k}, \qquad \theta_P^k(P) := \frac{\max_{0\le n<N}\left\lVert P\mathbf{e}_n^k \right\rVert_2}{E_\infty^k},\] and define the base-state weight by \[U_k^\perp(P) := \max_{0\le n<N}\left\lVert Q\mathbf{U}_n^k \right\rVert_2 .\] The next lemma combines these quantities in the projection-dependent factor \(\rho_k(P)\).

The linear time-periodic problem \[\label{eq:linear-periodic-problem} \mathbf{u}'(t)=A(t)\mathbf{u}(t)+\mathbf{g}(t), \qquad t\in[0,T], \qquad \mathbf{u}(0)=\mathbf{u}(T),\tag{23}\] where \(A(t)\in\mathbb{R}^{d\times d}\) and \(\mathbf{g}(t)\in\mathbb{R}^d\) are \(T\)-periodic, is the special case of 15 for which the fine and coarse propagators are affine maps of the initial state on each coarse time interval \([T_n,T_{n+1}]\), as holds for standard fixed-step integrators.

For this special case, the local smoothness conditions become global and their constants simplify. The discrepancy \(\mathcal{D}_n\) and the projection-based propagator \(\mathcal{H}_{P,n}\) have constant Jacobians, and one may take \(\mathcal{B}_k=\mathbb{R}^d\). Writing \(F_n\) and \(G_n\) for the linear parts of \(\mathcal{F}_n\) and \(\mathcal{G}_n\), the projection-based propagator has the linear part \(F_nP+G_nQ\); the constants of 4 reduce to \(\mu_k=\max_{0\le n<N}\left\lVert F_n-G_n \right\rVert_2\) and \(\nu_k=\kappa_k=0\), and the smallness condition ?? is vacuous. The cyclic operator \(\mathcal{C}_P\) is now affine, and its linearization \(\widehat{\mathcal{C}}_P\) no longer depends on the base state: it is the block-cyclic matrix with diagonal blocks \(I\) and cyclic subdiagonal blocks \(-(F_nP+G_nQ)\). The affine shifts cancel when two histories are subtracted, so the error history solves the linear system \[\widehat{\mathcal{C}}_P\,\mathbf{e}^{k+1} = \bigl((F_n-G_n)\,Q\mathbf{e}_n^k\bigr)_{n+1},\] with indices modulo \(N\). The update exists uniquely whenever \(\widehat{\mathcal{C}}_P\) is invertible. In this affine setting, the following assumption on the fine–coarse discrepancy \(\mathcal{D}_n\) replaces both the local nonlinear region and local smoothness assumptions ([ass:local-nonlinear-neighborhood,ass:local-discrepancy-smoothness]):

Under 6, \(\mu_k\le C_D\left(\Delta T\right)^{p+1}\), and the one-step factor of 1 reduces to \(\rho_k(P)=\widehat\Gamma_k\mu_k\theta_Q^k(P)\); the admissible radius \(R\) may be taken as any bound on \(E_\infty^{k+1}\) (such a bound exists under the assumed invertibility) and drops out of the estimate since \(\nu_k=\kappa_k=0\). For an arbitrary orthogonal projection, 1 records the resulting global estimate.

The projection enters the next-step error bound ?? through two quantities: the cyclic stability factor \(\Gamma\) and the unresolved-error norm \(\max_{0\le n<N}\left\lVert Q\mathbf{e}_n^k \right\rVert_2\).

4.2 Fourier-aware convergence↩︎

The unresolved-error quantity \(\left\lVert Q\mathbf{e}^k \right\rVert_\infty\) enters the projection-dependent factor \(\rho_k(P)\) in ?? through \(\theta_Q^k(P)\) and appears directly in the linear estimate ?? . To bound this quantity, consider a coarse-time history \(\mathbf{w}=(\mathbf{w}_0,\ldots,\mathbf{w}_{N-1})\in(\mathbb{R}^d)^N\) with temporal Fourier coefficients \(\widehat{\mathbf{w}}_\ell\) defined by 6 . For an orthogonal projection \(P\) we use the same symbol for its complex-linear extension to \(\mathbb{C}^d\), \(P(\mathbf{a}+\imath\mathbf{b}):=P\mathbf{a}+\imath P\mathbf{b}\), and the Euclidean norm on \(\mathbb{C}^d\) is \(\left\lVert \mathbf{a}+\imath\mathbf{b} \right\rVert_2^2=\left\lVert \mathbf{a} \right\rVert_2^2+\left\lVert \mathbf{b} \right\rVert_2^2\). For a selected temporal mode set \(\mathcal{I}\), define the tail and the selected-mode leak by \[\label{eq:linear-tail-definition} \operatorname{Tail}_{\mathcal{I}}(\mathbf{w}):=\sum_{\ell\notin\mathcal{I}}\left\lVert \widehat{\mathbf{w}}_\ell \right\rVert_2, \qquad \operatorname{Leak}_{\mathcal{I}}^{P}(\mathbf{w}):=\sum_{\ell\in\mathcal{I}}\left\lVert (I-P)\widehat{\mathbf{w}}_\ell \right\rVert_2 .\tag{25}\] The corresponding maximum-in-time estimate for the unresolved part \(Q\mathbf{w}\) of a history is stated next.

For the Fourier-aware projection, let \(\mathcal{I}^k\) be the selected temporal mode set and let \(P_{\mathcal{I}^k}\) be the orthogonal projection onto the mixed Fourier space \(\mathcal{V}_{\rm mix}^k(\mathcal{I}^k)\). If \(E_\infty^k>0\), define the Fourier tail–leak ratio \(\Theta_Q^k\) and the projected-error ratio \(\Theta_P^k\) by \[\label{eq:fourier-aware-ratios} \Theta_Q^k := \frac{ \operatorname{Tail}_{\mathcal{I}^k}(\mathbf{e}^k) + \operatorname{Leak}_{\mathcal{I}^k}^{P_{\mathcal{I}^k}}(\mathbf{e}^k) }{ \sqrt N\,E_\infty^k }, \qquad \Theta_P^k:=\theta_P^k(P_{\mathcal{I}^k}).\tag{26}\] 2 gives \(\theta_Q^k(P_{\mathcal{I}^k})\le\Theta_Q^k\).

The next result combines 1 with 2.

The factor \(\rho_{k,\mathrm F}\) in ?? contains three bracketed contributions. Through \(\Theta_Q^k\), the Fourier tail and leak control the unresolved-error term \(\mu_k\Theta_Q^k\); the base-state term \(\nu_kU_k^\perp(P_{\mathcal{I}^k})\Theta_P^k\) and discrepancy-curvature remainder \(\tfrac{\nu_k}{2}\Theta_Q^k\Theta_P^kE_\infty^k\) are the two nonlinear bracketed contributions. When each \(\mathcal{D}_n\) is affine, one may take \(\nu_k=0\), so only \(\mu_k\Theta_Q^k\) remains among these three bracketed contributions, as in the linear tail–leak estimate in 2. Since \(\Theta_Q^k\) and \(\Theta_P^k\) depend on the unknown error, they are a posteriori diagnostics. The factor \(\rho_{k,\mathrm F}\) may vary with \(k\) and need not be less than one; hence 1 gives a one-step estimate and does not by itself imply a uniform contraction.

This estimate separates temporal-mode selection from projection-space choice. Selecting all temporal modes gives \(\operatorname{Tail}_{\mathcal{I}^k}(\mathbf{e}^k)=0\). The equality \(\operatorname{Leak}_{\mathcal{I}^k}^{P_{\mathcal{I}^k}}(\mathbf{e}^k)=0\) holds if and only if the projection space contains both \(\operatorname{Re}\widehat{\mathbf{e}}_\ell^k\) and \(\operatorname{Im}\widehat{\mathbf{e}}_\ell^k\) for every \(\ell\in\mathcal{I}^k\). The solution-generated space \(\mathcal{V}_{\rm sol}^k(\mathcal{I}^k)\) alone need not contain these error coefficients. The source \(\mathcal{R}_{P,n}(\mathbf{u}_n^\star)-\mathcal{R}_{P,n}(\mathbf{U}_n^k)\) in 20 motivates enriching the projection space by \(\mathcal{V}_{\rm disc}^k(\mathcal{I}^k)\).

5 Numerical experiments↩︎

In this section, we present numerical experiments for Fourier-aware projection-based PP-PC (2; Fourier-aware PP-PC for short). The experiments serve three purposes:

1. to compare the outer-iteration convergence and projection-space dimensions of the proposed iteration with PP-PC [19] and Krylov-enhanced PP-PC [23] on two linear and two nonlinear time-periodic problems;

2. to evaluate the linear tail–leak estimate of 2 and illustrate, with frozen projection spaces, the roles of the unresolved-error norm and cyclic stability factor in 1;

3. to diagnose the three unweighted projection-dependent factors that enter the nonlinear one-step numerator in 1 and to evaluate the nonlinear tail–leak estimate of 1 using diagnostic estimates of the local constants.

Relative to the fine periodic reference \(\mathbf{u}^\star\) of 4, the reported error is \(E_\infty^k=\max_{0\le n<N} \left\lVert \mathbf{u}_n^\star-\mathbf{U}_n^k \right\rVert_2\). Each run starts from the periodic coarse solution \(\mathbf{U}^0\) and uses the problem-specific outer-iteration budget \(K\). All periodic solves are direct in the linear tests and use damped Newton shooting in the nonlinear tests, with a full-period residual tolerance of \(10^{-13}\). Newton Jacobians are assembled by forward differences with componentwise steps \(10^{-7}\max\{1,|x_j|\}\). All experiments were carried out in MATLAB R2024a on a laptop with an Intel Core i5-13500H CPU at 3.19 GHz and 32 GB RAM.

5.1 Linear problems↩︎

Both problems have period \(T=1\) and homogeneous Dirichlet conditions on \(x\in(0,1)\). We use second-order centered differences with \(\Delta x=1/32\); the first-order wave and realified Schrödinger systems both have dimension \(d=62\). On each coarse time interval, \(\mathcal{F}\) uses \(32\) Crank–Nicolson steps and \(\mathcal{G}\) one backward Euler step.

The wave problem is \[\label{eq:numerics-wave} \partial_{tt}u=\partial_{xx}u+f(x,t), \qquad f(x,t)=-\left(4\pi^2x(x-1)+2\right)\sin(2\pi t),\tag{27}\] for which a time-periodic solution is \(u(x,t)=x(x-1)\sin(2\pi t)\). We use \(N=16\), \(K=8\), and the fixed mode set \(\mathcal{I}_1=\{0,1,N-1\}\), which contains the known frequency pair. At \(T=1\), periodic even spatial modes make the continuous problem degenerate, while the odd-mode forcing admits the displayed solution. The discrete fine period map is nevertheless nondegenerate: the semidiscrete frequency closest to \(2\pi\) is \(64\sin(\pi/32)\approx6.273\ne2\pi\), and its nearest multiplier is about \(10^{-2}\) from \(1\).

The forced Schrödinger problem is \[\label{eq:numerics-schrodinger} \imath\,\partial_t\psi=-\partial_{xx}\psi+f(x,t),\tag{28}\] with periodic solution \[\psi(x,t) = \phi_1(x)e^{2\pi\imath t} +0.35\,\phi_2(x)e^{-4\pi\imath t},\] where \(\phi_1(x)=x(1-x)\) and \(\phi_2(x)=x(1-x)(1+0.25\sin\pi x)\). Writing \(\psi_h\) for its grid samples and \(L_h\) for the centered-difference Laplacian, we set \(f_h=\imath\,\dot{\psi}_h+L_h\psi_h\), so that \(\psi_h\) solves the semi-discrete system exactly. We use \(N=64\), \(K=16\), and, treating the frequencies as unknown, apply the adaptive increment-energy rule of 3.1 with a budget of five conjugate groups (\(\lvert\mathcal{I}^k\rvert\le10\)).

5.1.0.1 Convergence comparison

[fig:linear-convergence,tab:linear-summary] show that Fourier-aware PP-PC reaches \(E_\infty^k\le10^{-10}\) first: at \(k=8\) for the wave problem and \(k=3\) for the Schrödinger problem, versus \(k=15\) and \(k=8\) for Krylov-enhanced PP-PC; PP-PC reaches neither threshold within the respective budgets. The wave panel extends to \(k=15\) for the later Krylov crossing; the Schrödinger panel stops at \(k=7\). The table gives each nominal run’s final error, threshold iteration, maximal dimension, and mode count. A dash marks no threshold hit within \(K\) or no Fourier selection, as appropriate. Both projection spaces use the same relative singular-value tolerance \(10^{-12}\).

5.1.0.2 Evaluation of the tail–leak estimate

We evaluate 2 along the Fourier-aware iterations. For its cyclic factor \(\Gamma_k=\left\lVert \widehat{\mathcal{C}}_{P_{\mathcal{I}^k}}^{-1} \right\rVert_{\infty\to\infty}\), we use the computable block-row majorant \[\label{eq:block-row-majorant} \Gamma_k\le\overline{\Gamma}_k:=\max_{0\le n<N}\sum_{m=0}^{N-1} \left\lVert \bigl(\widehat{\mathcal{C}}_{P_{\mathcal{I}^k}}^{-1}\bigr)_{nm} \right\rVert_2.\tag{29}\] At the fixed discretization, the affine discrepancy has uniform Lipschitz constant \(\delta_{FG}:=\max_n\left\lVert F_n-G_n \right\rVert_2\), where \(F_n\) and \(G_n\) are the linear parts of \(\mathcal{F}_n\) and \(\mathcal{G}_n\). Applying the estimate with these computable quantities gives \[\label{eq:numerics-evaluated-bound} B_{k+1} := \overline{\Gamma}_k\,\delta_{FG}\, \frac{\operatorname{Tail}_{\mathcal{I}^k}\left(\mathbf{e}^k\right) +\operatorname{Leak}_{\mathcal{I}^k}^{P_{\mathcal{I}^k}}\left(\mathbf{e}^k\right)}{\sqrt N}.\tag{30}\] Thus \(E_\infty^{k+1}\le B_{k+1}\) in exact arithmetic. The tail and leak use the measured error history \(\mathbf{e}^k\), so \(B_{k+1}\) is a posteriori. The estimate therefore applies to the Crank–Nicolson wave propagator without a fine-propagator contraction assumption.

[tab:linear-tail-leak,fig:linear-bound] show that \(B_k\) lies above \(E_\infty^k\) at every computed step and follows the initial error decay. The bound is not sharp, remaining roughly two to five orders above the measured error. The prefactor \(\overline{\Gamma}_k\delta_{FG}\) contributes to this slack. The block-row majorant \(\overline{\Gamma}_k\) ranges over \(8\times10^1\)–\(4\times10^3\) for the wave problem and \(3\times10^1\)–\(3\times10^2\) for the Schrödinger problem, while \(\delta_{FG}\approx63\) and \(1.2\), respectively. The larger wave values of \(\overline{\Gamma}_k\) are consistent with its near-resonant discretization.

The tail and leak columns separate the two contributions to the unresolved-error bound. For the fixed-band wave problem, the Fourier-aware tail stays at roundoff level and the leak decays. For the adaptive Schrödinger problem, \(\mathcal{I}^0\) is empty, so all initial error modes contribute to the tail. After the first selection, the tail falls to roundoff level and the leak decays. For comparison, the table also reports Krylov-enhanced diagnostics in the all-mode representation of 1; its diagnostic tail is zero, but its leak need not be. Thus a small tail alone does not make the tail–leak bound small.

5.1.0.3 A posteriori test of projection dependence

For a fixed \(P\), 1 identifies two projection-dependent quantities in the next-step error bound: the unresolved-error norm \(\max_{0\le n<N}\left\lVert Q\mathbf{e}_n^k \right\rVert_2\), with \(Q=I-P\), and the cyclic stability factor \(\Gamma\) associated with \(\widehat{\mathcal{C}}_P\). To examine these quantities a posteriori, we use two frozen projection spaces. We construct \(P_F\) from the mixed Fourier space and \(P_K\) from the solution-snapshot space, using their respective histories through the same iteration \(k_\ast\). For each problem, \(k_\ast\) is the first iteration at which Fourier-aware PP-PC satisfies \(E_\infty^{k_\ast}\le10^{-10}\). We restart the same discrepancy-form update 3 from \(\mathbf{U}^0\) with each projection fixed, so only \(P\) differs. 3 reports \(\max_n\left\lVert Q\mathbf{e}_n^0 \right\rVert_2\), the corresponding block-row majorant \(\overline{\Gamma}(P)\) obtained from 29 , and \(E_\infty^1\); 5 shows the subsequent fixed-projection histories.

Within each problem, the computed cyclic majorants are comparable, whereas the unresolved initial error for \(P_F\) is more than five orders of magnitude smaller than for \(P_K\). Consistently, \(P_F\) produces one-step errors of order \(10^{-12}\), while \(P_K\) produces errors of order \(10^{-5}\) and \(10^{-7}\) for the wave and Schrödinger problems, respectively. The subsequent fixed-projection histories preserve this separation. Within the evaluated estimate, the unresolved factor therefore provides the dominant separation between the two projections. Thus, in these tests, the reported Fourier-aware construction yields a more refined projection space for the error-relevant components.

5.2 Nonlinear problems↩︎

Both nonlinear problems have period \(T=1\) and periodic boundary conditions on \([0,1)\). We use second-order centered differences on \(32\) grid points (\(d=64\)), \(N=16\), and \(K=8\). The fine propagator is classical fourth-order Runge–Kutta with \(\delta t=1/512\). Writing the semi-discrete system as \(\mathbf{w}'=A\mathbf{w}+\mathbf{g}(\mathbf{w},t)\), the coarse propagator is an implicit–explicit (IMEX) Euler step [29], [30] \[\label{eq:numerics-semi-implicit} \mathcal{G}_n(\mathbf{w}) = \left(I-\Delta T\,A\right)^{-1} \left(\mathbf{w}+\Delta T\,\mathbf{g}(\mathbf{w},T_{n+1})\right), \qquad \Delta T=1/16.\tag{31}\] Thus the stiff linear part is implicit and the nonlinear part explicit, with one fixed-matrix solve per coarse step.

The first problem is the periodically forced Stuart–Landau reaction–diffusion system [31]. For \(w=p+\imath q\), \[\label{eq:numerics-stuart-landau} \begin{align} \partial_t w &= D\,\partial_{xx}w+(\lambda+\imath\omega)w -(1+\imath\beta)|w|^2w+f(x,t),\\ f(x,t) &= 0.5\left(\sin2\pi x+0.35\sin4\pi x\right)\cos2\pi t, \end{align}\tag{32}\] where \(D=0.05\), \(\lambda=0.5\), \(\omega=2\pi\), and \(\beta=1\).

The second problem is the periodically forced Brusselator [32], with \(a=1\), \(b=3\), \(D_1=10^{-2}\), \(D_2=5\times10^{-3}\), and the same forcing \(f\): \[\label{eq:numerics-brusselator} \begin{align} \partial_t u &= D_1\,\partial_{xx}u+a-(b+1)u+u^2v+f(x,t),\\ \partial_t v &= D_2\,\partial_{xx}v+bu-u^2v. \end{align}\tag{33}\] Although the forcing contains only the temporal frequency \(\ell=1\), the nonlinear terms generate higher harmonics. Both tests therefore use the fixed 7-mode band \(\mathcal{I}_3=\{0,1,2,3,N-3,N-2,N-1\}\).

The two problems exercise the two stability mechanisms evaluated below. For Stuart–Landau, the measured linearization of \(\mathcal{H}_{P^k,n}\) along the fine reference is contractive, so its evaluated bound uses the loop-gain specialization associated with 5. For the Brusselator, the corresponding linearization is noncontractive, so its evaluated bound uses the general cyclic-inverse form. The local nondegeneracy condition of 3 and the smallness condition ?? are examined below. The cyclic-inverse analysis therefore covers noncontractive settings in which the linearized cyclic operator remains invertible, beyond the contraction-based stability mechanism used in earlier quantitative nonlinear PP-PC estimates [19].

5.2.0.1 Convergence comparison

[fig:nonlinear-convergence,tab:nonlinear-summary] show that Fourier-aware PP-PC reaches \(E_\infty^k\le10^{-10}\) first, at \(k=2\) for Stuart–Landau and \(k=4\) for the Brusselator, and is the only method to do so in both tests within \(K=8\). Under the same rank-revealing procedure, its maximal projection-space dimension is smaller than that of Krylov-enhanced PP-PC in both tests. The panels show the histories through \(k=6\), while the table summarizes the complete \(K=8\) runs by their final errors, threshold iterations, maximal dimensions, and mode counts.

5.2.0.2 Projection-level decomposition of the local one-step factor

For \(P=P_{\mathcal{I}^k}\), set \(Q=I-P_{\mathcal{I}^k}\) and define \(m_Q^k:=\max_{0\le n<N}\left\lVert Q\mathbf{e}_n^k \right\rVert_2\) and \(m_P^k:=\max_{0\le n<N}\left\lVert P_{\mathcal{I}^k}\mathbf{e}_n^k \right\rVert_2\); abbreviate \(U_k^\perp:=U_k^\perp(P_{\mathcal{I}^k})\). For the Fourier-aware projection, these definitions and 2 give \[\begin{align} m_Q^k &=E_\infty^k\theta_Q^k(P_{\mathcal{I}^k}) \le \widetilde{m}_Q^k :=\frac{\operatorname{Tail}_{\mathcal{I}^k}(\mathbf{e}^k) +\operatorname{Leak}_{\mathcal{I}^k}^{P_{\mathcal{I}^k}}(\mathbf{e}^k)}{\sqrt N} =E_\infty^k\Theta_Q^k,\\ m_P^k&=E_\infty^k\Theta_P^k. \end{align}\] Thus, after multiplication by \(E_\infty^k\), the factor in 1 has the three numerator contributions \[\mu_k\,m_Q^k, \qquad \nu_k\,U_k^\perp m_P^k, \qquad \frac{\nu_k}{2}m_Q^k m_P^k.\] 5 reports this majorant, the three unweighted factors, and the resulting \(E_\infty^{k+1}\).

In both problems, \(m_Q^k\) is the largest unweighted factor. The measured weights preserve this dominance. The discrepancy-curvature factor \(m_Q^km_P^k\) is smaller and quadratic in the error, while \(U_k^\perp m_P^k\) is negligible. Over the reported steps, \(\Theta_Q^k=\widetilde{m}_Q^k/E_\infty^k\le4.93\times10^{-2}\), so the normalized tail–leak sum is small. Since \(\widetilde{m}_Q^k\) majorizes the dominant unresolved factor, these data provide a posteriori evidence that a small tail together with a small leak is associated with the observed rapid next-step decay; the complete estimate is evaluated below.

5.2.0.3 Evaluation of the nonlinear tail–leak estimate

Thus 1 replaces \(m_Q^k\) by its tail–leak majorant \(\widetilde{m}_Q^k\) in the unresolved-error and discrepancy-curvature contributions, while leaving the base-state contribution unchanged.

The corresponding diagnostic weighted sum is \[\mathcal{S}_{k,\mathrm F} :=\hat{\mu}\,\widetilde{m}_Q^k +\hat{\nu}\,U_k^\perp m_P^k +\frac{\hat{\nu}}{2}\widetilde{m}_Q^k m_P^k.\]

For Stuart–Landau, the computed finite-difference (FD) reference-Jacobian norm \[\hat{\alpha}_k :=\max_n\left\lVert J_{\mathcal{H}_{P_{\mathcal{I}^k},n}}^{\rm FD} (\mathbf{u}_n^\star) \right\rVert_2 \approx0.956<1,\] giving the loop-gain specialization \[\label{eq:numerics-estimated-nl-bound} B_{k+1,\mathrm F}^{\mathrm{nl}} :=\frac{\mathcal{S}_{k,\mathrm F}}{(1-\hat{\alpha}_k)-\tfrac12\hat{\kappa}_kR}.\tag{34}\] For the Brusselator, whose reference linearization is noncontractive, we use the general cyclic-inverse form \[\label{eq:numerics-nk-nl-bound} B_{k+1,\mathrm F}^{\mathrm{nl}} :=\frac{\overline{\widehat\Gamma}_k}{1-\tfrac12\overline{\widehat\Gamma}_k\hat{\kappa}_kR} \mathcal{S}_{k,\mathrm F}.\tag{35}\]

Here \(\overline{\widehat\Gamma}_k\) is defined analogously by the block-row formula in 29 , using the assembled finite-difference cyclic matrix. With exact Jacobians, the same construction majorizes \(\widehat\Gamma_k=\left\lVert \widehat{\mathcal{C}}_{P_{\mathcal{I}^k}}^{-1} \right\rVert_{\infty\to\infty}\); the prefactor is increasing in this norm while its denominator is positive, so the majorization preserves the exact bound.

We take \(R=E_\infty^0\). The Stuart–Landau denominator is positive, while for the Brusselator \(\overline{\widehat\Gamma}_k\in[62,66]\) and \(\tfrac12\overline{\widehat\Gamma}_k\hat{\kappa}_kR \in[0.50,0.66]\); moreover, \(E_\infty^{k+1}\le R\) holds throughout both runs. All diagnostic Jacobians use central differences with step \(10^{-6}(1+\left\lVert \mathbf{x} \right\rVert_\infty)\). Those at \(\mathbf{u}^\star\) give \(\hat{\alpha}_k\) and the cyclic matrix used to form \(\overline{\widehat\Gamma}_k\), while \(\hat{\mu}\) maximizes discrepancy-Jacobian norms at \(\mathbf{u}^\star\) and \(\mathbf{U}^0\). The estimate \(\hat{\nu}\) uses Jacobian difference quotients between these two histories and along two directions generated with seed \(0\) and scaled to radius \(E_\infty^0\) at each time point; the same directional sampling of the Jacobians of \(\mathcal{H}_{P_{\mathcal{I}^k},n}\) gives the sampled curvature \(\hat{\kappa}_k\). These finite-difference and sampled quantities are diagnostic estimates rather than certified bounds for the exact local constants.

7 shows that both evaluated estimates lie above the realized errors and follow their decay. Relative to the projection-level estimate, the tail–leak majorization adds, over the evaluated iterations, at most factors \(3.97\) and \(1.45\) for Stuart–Landau and the Brusselator, respectively. Thus 5 isolates the local projection mechanism, while 7 evaluates the corresponding Fourier tail–leak estimate of 1, including the noncontractive regime.

5.3 Summary and discussion↩︎

[tab:linear-summary,tab:nonlinear-summary] show that Fourier-aware PP-PC reaches \(E_\infty^k\le10^{-10}\) in the fewest outer iterations in all four examples. Its maximal projection-space dimension is smaller than that of Krylov-enhanced PP-PC in three examples. For the wave problem it is larger, but the tolerance is reached at \(k=8\) rather than \(k=15\) (3). The frozen-projection experiment illustrates the dependence on \(P\) in 1: at the same history depth, the mixed Fourier spaces leave smaller unresolved errors and produce smaller one-step errors under comparable computed cyclic majorants (3 5).

For the linear problems, the measured tail–leak quantities and evaluated bound are consistent with 2 (2 4). For the nonlinear problems, 5 identifies \(m_Q^k\) as the largest unweighted projection-level factor. Over the reported steps, \(\Theta_Q^k\le4.93\times10^{-2}\). 7 evaluates the nonlinear tail–leak estimate using diagnostic estimates of the local constants for the Stuart–Landau and Brusselator cases, whose computed finite-difference reference linearizations are contractive and noncontractive, respectively. Taken together, these results indicate that the mixed Fourier construction yields a more refined projection space in the reported tests: it is better aligned with error-relevant directions but not always smaller.

6 Conclusion and future work↩︎

We introduced Fourier-aware projection-based PP-PC, a projection-based periodic parareal method that combines a discrepancy-based correction scheme with a Fourier-aware projection space. We developed a convergence analysis for general nonlinear time-periodic problems, in which a temporal tail–leak bound controls both the unresolved-error term and the nonlinear remainder in a local one-step estimate, and showed that the linear case reduces to global bounds under weaker assumptions. We demonstrated on linear and nonlinear problems that Fourier-aware PP-PC converges faster than Krylov-enhanced PP-PC, with the measured errors consistent with the analysis.

Our comparisons measure outer-iteration convergence, as in [19], [23]. Since the corrected periodic solve is more expensive per iteration than the coarse solve of the original PP-PC, a work-normalized timing study and dedicated fast solvers for it are subjects of further study.

Acknowledgments↩︎

This work was supported by the National Key R&D Program of China under Grant Nos. 2020YFA0711900 and 2020YFA0711902. During the preparation of this manuscript, the authors used OpenAI’s ChatGPT to improve its readability and assist with the design of code for the numerical experiments. The authors assume responsibility for all content.

References↩︎

---

1. School of Mathematical Sciences, Fudan University, Shanghai, China (, , , ).↩︎