The exposition follows in style and idea the paper "M. Weiser. Faster SDC convergence on non-equidistant grids with DIRK sweeps. BIT Numerical Mathematics 55 (4): 1219-1241, 2015", but differs in details. In particular, we focus on QSDC and [0,1] here.

From ODEs to iteration matrices

We consider \( \dot u = f(u) \), or its Picard reformulation

\[ u(t) = u(0) + \int_0^{t} f(u) \, d\tau. \]

For a given approximation \( u^k \in \mathbb{P}_{n+1} \), we consider the defect equation on a collocation grid \( 0=t_0 < \dots < t_n \le 1 \),

\[ \begin{aligned} \delta u^k(t_i) &= u(t_i)-u^k(t_i) \\ &= \delta u^k(t_{i-1}) + \int_{t_{i-1}}^{t_i} f(u) \, d\tau + u^k(t_{i-1})-u^k(t_i)\\ &\approx \delta u^k(t_{i-1}) + \int_{t_{i-1}}^{t_i} f(u^k) + f'(u^k)\delta u^k \, d\tau + u^k(t_{i-1})-u^k(t_i) \end{aligned}\]

with \( \delta u^k(t_0) = 0 \). Integrating the right hand side \( f(u^k)\) with a high order quadrature formula and the implicit term \( f'(u^k)\delta u^k \) with a lower order formula relying on values up to \( t_i \) only, we obtain

\[ \delta u^k_i - \delta u^k_{i-1} \approx \sum_{j=0}^n S_{ij} f(u^k_j) + u^k_{i-1}-u^k_i + \sum_{j=0}^i \hat S_{ij} f'(u^k_j) \delta u^k_j, \quad i=1,\dots,n, \quad \delta^k_0 = 0. \]

With matrices

\[ D = \begin{bmatrix} 1 \\ -1 & 1 \\ & \ddots & \ddots \\ && -1 & 1 \end{bmatrix} , \quad S = \begin{bmatrix} 0 & \dots & 0 \\ S_{10} & \dots & S_{1n} \\ \vdots & \ddots & \vdots \\ S_{n0} & \dots & S_{nn} \end{bmatrix}, \quad \hat S = \begin{bmatrix} 0 & & & 0 \\ \hat S_{10} & \hat S_{11} \\ \vdots & \vdots & \ddots \\ \hat S_{n0} & \hat S_{n1} & \dots & \hat S_{nn} \end{bmatrix} \in \mathbb{R}^{n+1 \times n+1},\]

this can be written in vectorized form as

\[ D \delta u - \hat S F \delta u^k = S f(u^k) - D u^k + s, \quad F = \text{diag}(f'(u^k)), \quad s = \begin{bmatrix} u^k_0 \\ 0 \end{bmatrix}. \]

With the update \( u^{k+1} = u^k + \delta u^k \), we get

\[ u^{k+1} = u^k + (D-\hat S F)^{-1}(Sf(u^k) - Du^k + s) \]

On Dahlquist's test equation \( \dot u = \lambda u \), this reads

\[ u^{k+1} = \left( I + (D-\lambda \hat S)^{-1}(\lambda S - D) \right) u^k + (D-\lambda \hat S)^{-1} s. \]

Consequently, the iteration matrix of the stationary iteration is

\[ G = I - (D-\lambda \hat S)^{-1}(D-\lambda S). \]

Choosing the approximate integration matrix

In QSDC, \( \hat S \) (and also \( D \) on one side of the equation, but we won't consider this here) can be chosen freely as long as the triangular structure is preserved. The standard choice is the right-looking rectangular rule corresponding to the implicit Euler sweep, i.e.

\[ \hat S = \begin{bmatrix} 0 & 0 & \dots & 0 \\ 0 & t_1-t_0 & 0 & \dots \\ 0 & 0 & t_2-t_1 & 0 \\ \vdots & \vdots & &\ddots \end{bmatrix}. \]

Design spaces

We can choose \( \hat S \) in restricted design spaces imposing certain favorable properties, and then optimize numerically. For example, we may select

Uniform Euler sweep

Due to the triangular structure of \( \hat S \), we need to solve systems of the form \( M - \tau \hat S_{ii} f'(u_i^k) \). Often, \(f'(u_i^k)=A\) can be assumed to be constant. If a factorization of \( M - \tau \hat S_{ii} A \) is to be used, it would be great if \( \hat S_{ii} \equiv s \) holds for some fixed \( s \), such that only one factorization needs to be constructed and can be used over all Euler steps of the sweep. We thus look for

\[ \hat S = \begin{bmatrix} 0 & 0 & \dots & 0 \\ 0 & s & 0 & \dots \\ 0 & 0 & s & 0 \\ \vdots & \vdots & &\ddots \end{bmatrix}. \]

Optimization criteria

As a means to select good approximate integration matrices, we make use of Dahlquist's test equation \( \dot u = \lambda u \), here for \( \lambda \in ]-\infty, 0[ \), as we are focusing on parabolic equations. For the iteration matrix \( G(\lambda) = I - (D-\lambda \hat S)^{-1}(D-\lambda S) \), we consider three quantities:

• its norm \( a=\|G\|_2 \) as an upper bound for the contraction
• its mean contraction \( b=\|G^m\|_2^{1/m} \) over \( m \) iterations, where \( m \) is chosen as the number of SDC iterations to be performed, usually the collocation order
• its spectral radius \( c=\max_{s\in \sigma(G)} |s| \) as a lower bound for the contraction

For each of these quantities, we average over a relevant portion of the negative real line, i.e.

\[ \bar a = \frac{1}{|\Lambda|} \int_{\lambda\in \Lambda} e^\lambda a \, d\lambda, \]

with \( \Lambda = [-10^{3}, -10^{-3}]\).