A hybrid algorithm based on parareal and Schwarz waveform relaxation

1.   Introduction

For time dependent partial differential equations (PDEs), Schwarz waveform relaxation (SWR) is a powerful technique for realizing fast computation. It belongs to the widely used domain decomposition (DD) methods ^[1,2] but with a completely new implementation strategy: the classical strategy employs the alternating or parallel Schwarz method to the elliptic PDEs which result from semi-discretization of the original PDEs in time ^[3,4,5]; but, in the SWR framework one directly decomposes the space-time domain and solve each subproblem simultaneously or alternately. The main features (or say merits) of the SWR algorithms are able to treat different subproblems numerically differently with an adapted procedure for each subdomain. One can refer to ^{[6,7,8,9,10,11,12,13]} for more details about the SWR algorithm.

The aforementioned SWR algorithm realizes parallelism across space. Next, we introduce another efficient parallel algorithm, namely 'parareal', which realizes parallelism across time steps. This method was invented by Lions, Maday and Turinici in ^[14] as a numerical method to solve evolution problems in parallel. The name was chosen to indicate that the algorithm is well suited for parallel real time computations of evolution problems whose solution cannot be obtained in real time using one processor only. The method approximates successfully the solution later in time before having fully accurate approximations from earlier times. Successful applications of this method include power system simulations ^[15,16], differential algebraic equations ^[17], complex chemical reactions ^[18,19], optimal control problems ^{[20,21,22,23]} and Volterra integro-differential equations ^[24,25,26]. Convergence analysis of the parareal algorithm can be found in ^{[27,28,29,30]}.

It is natural to combine a parareal with other algorithms to formulate a hybrid algorithm with higher parallelism. For example, in ^{[31,32,33,34,35,36,37,38,39]} the classical waveform relaxation (WR) algorithm ^[40] is embedded into the parareal framework to reduce the computational cost. The authors proposed the "parareal-WR" algorithm and defined it by using the $k_0$-iteration WR algorithm as the fine propagator for the parareal algorithm. The algorithm is analyzed at a continuous level and superlinear and linear convergence were proved on short and long time intervals, respectively. In this paper we combine the SWR technique and the parareal algorithm to construct a hybrid algorithm, namely, Para-SWR, for time dependent PDEs, which can realize parallelism both in space and time directions. Combining a parareal with SWR has already been addressed in recent years; see, e.g., ^[41,42,43]. In those previous works, the SWR algorithm was used as the so-called fine propagator for the parareal algorithm and this leads to inner-outer iterations. This strategy is different from ours, since in our Para-SWR algorithm we use a single iteration of the parareal as a time integrator for the SWR algorithm. This actually leads to a special discrete SWR algorithm which can be used in parallel in terms of both space and time.

Here, we focus on analyzing the Para-SWR algorithm at the discrete level and deriving suitable conditions for convergence. We also investigate how the convergence rate of the Para-SWR algorithm behaves with respect to the number of subdomains and the number of parallelled time intervals. It is inspiring to report that the proposed Para-SWR algorithm posses the weak scaling property, which implies that increasing the number of subdomains and the number of parallelled time intervals has no (or slight) influence on the convergence rate. An algorithm scales weakly if it can solve a larger problem in reasonable time by increasing the number of processors.

The remainder of this paper is organized as follows. In Section 2, we describe the Para-SWR algorithm investigated in this paper. In Section 3, we analyze the Para-SWR algorithm and present suitable convergence conditions. Weak scaling of the algorithm is also investigated. Section 4 provides some numerical results to validate our theoretical predictions. We finish this paper in Section 5 with some conclusion remarks.

2.   The Para-SWR hybrid algorithm

We consider the following 1-D linear reaction diffusion equation as the model problem

where $a$ and $\nu (\neq0)$ are constants. Without loss of generality, we assume $a > 0$ since otherwise we can make a change of variables $v = ue^{-\tau t}$ with $\tau+a > 0$ which leads to (2.1) with a positive reaction coefficient. We first introduce the SWR algorithm at a continuous level to (2.1). We let $\Delta x$ be a suitable small parameter and $L = [N\times(m-1)+2]\times\Delta x$ with some integers $N(\geq2)$ and $m(\geq2)$. We then decompose the space domain $[0, L]$ into $N$ subdomains

In this situation, the overlap size between every two adjacent subdomains is $2\Delta x$. Then, the $N$-subdomain SWR algorithm for (2.1) can be written as

where $p$ is a free parameter. Note that, the transmission condition

can be equivalently rewritten as

Hence, by letting $p\rightarrow +\infty$ we get

This is the so-called Dirichlet transmission condition and the corresponding algorithm is called classical SWR algorithm (see ^[9]). We next use the central finite difference method with mesh size $\Delta x$ to discretize the SWR algorithm (2.3) and this leads to the following semi-discrete SWR algorithm:

where $i$ varies from $(j-1)(m-1)+1$ to $j(m-1)+1$ and $j = 2, 3, \dots, N-1$. Similarly, for the left subdomain $\Omega_1$ and the right subdomain $\Omega_{N}$

Note that, $v^k_j(i, t) = u^k_j(i\Delta x, t)+\mathcal{O}(\Delta x^2)$ if $i\Delta x\in(L_j, R_j)$ and $v^k_j(i, t) = u^k_j(i\Delta x, t)+\mathcal{O}(\Delta x)$ if $i\Delta x = L_j$ or $i\Delta x = R_j$. However, upon convergence it holds that $v^{\infty}_j(i, t) = u^{\infty}_j(i\Delta x, t)+\mathcal{O}(\Delta x^2) = u_j(i\Delta x, t)+\mathcal{O}(\Delta x^2)$ for $i\Delta x\in[L_j, R_j]$, that is the discretization of the artificial boundary conditions does not effect the accuracy of the converged solutions ^[44]. Define

Then the semi-discrete SWR algorithm (2.4) can be equivalently rewritten as

Clearly, this can be written in a more concise form

where

We next introduce the parareal algorithm. For a general system of ordinary differential equations

the parareal algorithm can be described as follows. First, the whole time interval $[0, T]$ is divided into $N_t$ time slices $[T_n, T_{n+1}]$, $n = 0, 1, \dots, N_t-1$. We suppose that all time slices are of uniform size, i.e., $T_{n+1}-T_n = \Delta T = \frac{T}{N_t}$. Second, each large time slice $[T_n, T_{n+1}]$ is divided into $J(\geq2)$ small time slices. Then, two numerical propagators $\mathcal{G}$ and $\mathcal{F}$ are assigned to the coarse and fine time grids, respectively. We designate by the symbol $\ominus$ the time-sequential implementation, and by the symbol $\oplus$ the time parallel implementation. With the aforementioned introduction, the concrete parareal algorithm (for Problem (2.8)) proposed by Lions, Maday and Turinici ^[14] is given in Algorithm 2.1.

Clearly, the argument $\tilde{v}_{n+1}$ can be written as $\tilde{v}_{n+1} = \mathcal{F}^J\left(T_n, v^k_n, \Delta t\right)$ and therefore Algorithm 2.1 can be written compactly as

where $\Delta t = \frac{\Delta T}{J}$ and $\mathcal{F}^J\left(T_n, v^{k}_n, \Delta t\right)$ denotes a value obtained by running $J$ steps of the fine propagator $\mathcal{F}$ with the initial value $v^k_n$ and the fine step size $\Delta t$.

Our Para-SWR algorithm consists of applying the one-iteration parareal algorithm to System (2.6), which is derived from the semi-discrete SWR algorithm. We use the backward Euler method as both the coarse and the fine propagators. On each time interval $[T_n, T_{n+1}]$, we need to provide a seed value to start up the parallel computation and it is natural to use the $\mathcal{V}^{k-1}_n$ obtained in the previous iteration as the seed. Then, the parallel computation in each subinterval $[T_n, T_{n+1}]$ is

where ${\mathbb{V}^{k}_{n, 0} = \mathcal{V}^{k-1}_{n}}$ and $\Delta T = J\Delta t$. From (2.10) we get

With the values $\mathbb{V}^{k}_{n+1}$ and $\mathbb{V}^{k}_{n, 1}$, the serial correction is

where the argument $\mathcal{G}\left(T_n, \mathcal{V}^k_{n}, \Delta T\right)$ denotes the value generated by the backward Euler method with the starting value $\mathcal{V}^k_{n}$, that is

It then follows by substituting (2.11) and (2.13) into (2.12) that

Remark 1. The uniform formula (2.6) of the semi-discrete SWR algorithm leads to (2.14) which will greatly facilitate the theoretical analysis done in the next section. In practical computation, it is convenient to directly apply the parareal algorithm on each subdomain $\Omega_j\times[0, T]$, $j = 1, 2, \dots, N$.

3.   Analysis of the Para-SWR algorithm

In this section, we focus on investigating the theoretical property of the Para-SWR algorithm. The analysis is divided into two parts, i.e., the convergence of the algorithm and the weak scaling of the algorithm. In the sequel, we assume $\Delta x p\geq1$ and this implies $q\geq0$.

3.1.   Convergence analysis

To analyze the convergence of the Para-SWR iteration (2.14), for simplicity and without loss of generality we set $\tilde{f}(t) = 0$. Let

and by using the relation $\mathbb{V}^{k}_{n, 0} = \mathcal{V}^{k-1}_{n}$ we know from the second equality of (2.11) that

It then follows by substituting (3.2) into (2.14) that

Define

Under the assumptions $q\geq0$ and $a > 0$, it is easy to know that $\lambda_c < \lambda_f < 1$. From the second equality of (3.2) we get

From (3.3) we have

where $\mathbb{E}_1$ denotes the last row of the matrix $\mathbb{E}$. Define

Then from (3.5) and (3.6) it is easy to get

Define

Let the vectors $\{z^k\}$ be generated by the following recurrence relation

with the initial vector $z^0 = \begin{pmatrix}\mathcal{V}^0\\ \epsilon^0\end{pmatrix}$. Since the matrix $\mathbb{M}$ is nonnegative, it is clear that $\begin{pmatrix}\mathcal{V}^k\\ \epsilon^k\end{pmatrix}\leq z^k$ in the component sense for all $k\geq1$. Thus, it is sufficient to investigate the convergence of the iteration (3.10). It is well known that the iteration converges to zero if and only if the spectral radius of $\mathbb{M}$ is less than one, i.e., $\rho(\mathbb{M}) < 1$. To analyze the spectral radius of $\mathbb{M}$, we need the following definitions and lemmas concerning nonnegative matrices.

Definiton 1. [45,page 95] For $n\times n$ real matrices $P$, $P_1$ and $P_2$, the relation $P = P_1-P_2$ is called a regular splitting of $P$ if $P_1$ is nonsingular with $P^{-1}_1\geq0$ and $P_2\geq0$.

Definiton 2. [45,pages 53-54] Let $P$ be a $n\times n$ real matrix and $G(P_B)$ be a directed graph associated with the matrix $P$, which has $n$ vertices and an edge from the vertex $i$ to vertex $j$ precisely when $P_{ij} > 0$. Then the matrix $P$ is irreducible if and only if its associated graph $G(P)$ is strongly connected, that is for any two vertices $i$ and $j$ of $G(P)$, $G(P)$ contains a path from $i$ to $j$.

Lemma 1. [45,page 90] Let $P = (p_{ij})$ be a $n\times n$ real matrix with $p_{ij}\leq0$ for all $i\neq j$; then, the following statements are equivalent:

1) $P$ is nonsingular and $P^{-1}\geq0$;

2) the diagonal entries of $P$ are positive real numbers and letting $P_D$ be the diagonal matrix whose diagonal entries $d_{ii}$ are defined as

then the matrix $P_B: = I-P_D^{-1}P$ is nonnegative and $\rho(P_B) < 1$.

Lemma 2. [45,page 51] Let $P\in\mathbb{R}^{n\times n}$ and $P = P_1-P_2$ be a regular splitting of $P$. Then $P$ is nonsingular with $P^{-1}\geq0$ if and only if $\rho(P_1^{-1}P_2) < 1$.

Lemma 3. [45,page 50] Let $P\in\mathbb{R}^{n\times n}$ be an irreducible non-negative matrix with the spectral radius $\rho(P)$. Then, the number $\rho(P)$ is a positive real number and an eigenvalue of the matrix $P$, called the Perron-Frobenius eigenvalue.

Theorem 1. Let $q\geq0$, $a > 0$, $\mu\lambda_f < 1$ and $\beta < 1$. Then the Para-SWR algorithm based on the backward Euler method is convergent and the convergence rate can be bounded by $\rho(\mathbb{M})$, provided the following condition is satisfied

where $\lambda_f$ and $\mu$ are defined by (3.4) and $\beta$ and $\eta$ are defined by (3.9).

Proof. It is sufficient to prove that Condition (3.11) implies $\rho(\mathbb{M}) < 1$. In fact, this is an "if and only if" condition. Let

It is easy to check that $Q = (I-\lambda_fD)^{-1}$. Let $X_1 = \eta\lambda_f\mathbb{E}_1^T$ and $X_2 = \lambda_f(I-\lambda_fD)\mathbb{D}$ and we have

Define

Then we have

From Definition 1, we know that $P = P_1-P_2$ is a regular splitting of the matrix $P$. Since $\mu\lambda_f < 1$ and $\beta < 1$, it is clear that the matrix $P_B$ is nonnegative. Hence, we have

Our ultimate goal is to prove that $\rho(\mathbb{M}) < 1$, and that due to (3.15) it is sufficient to prove that $\rho(P_B) < 1$. Our proof toward $\rho(P_B) < 1$ is divided into two parts.

Part 1: irreducibility of the matrix $P_B$. The directed graph $G(P_B)$ associated with the matrix $P_B$ defined by (3.14) is shown in Figure 1. For the arbitrarily chosen vertices $i$ and $j$, the path from $i$ to $j$ can be chosen as

● if $i > j$: the path is $(i\rightarrow i-1\rightarrow i-2\rightarrow \cdots\rightarrow j)$;

● if $i < j$: the path is $(i\rightarrow i-1\rightarrow i-2\rightarrow \cdots\rightarrow 2\rightarrow1\rightarrow j)$.

Hence, $G(P_B)$ is strongly connected. From Definition 2, we know that the matrix $P_B$ is irreducible.

Part 2: the characteristic polynomial of the matrix $P_B$. Let $R(\lambda) = \Delta (P_B-\lambda I)$ be the characteristic polynomial of the matrix $P_B$, where $\lambda$ is a scale number and for any square matrix $A$, $\Delta (A)$ denotes the determinant of $A$. To derive an explicit expression of $R(\lambda)$, we recall a well known result which states that

where $R_{1, j}$ is a $J\times J$ matrix obtained by removing the first row and the $j$-th column from the matrix $P_B-\lambda I$. Clearly, $\Delta (R_{1, 1}) = (-\lambda)^J$ and $\Delta (R_{1, 2}) = \frac{\lambda_f}{1-\mu\lambda_f}(-\lambda)^{J-1}$ since both $R_{1, 1}$ and $R_{1, 2}$ are lower triangular matrices. But the matrices $R_{1, j}$ with $3\leq j\leq J+1$ are not lower triangular. For any $j\in\{3, 4, \dots, J, J+1\}$, we have

where the entries "$\sigma, 0$" corresponds to the $(j-1)$-th row of the matrix $R_{1, j}$. It follows by applying the expansion like (3.16) to $R_{1, j}$ that $\Delta (R_{1, j}) = \sigma^{j-1}(-\lambda)^{J-j+1}$. Hence, substituting these expressions into (3.16) gives

From (3.17), we know that the nonzero root of $\Delta(P_B-\lambda I) = 0$ satisfies

Now, by combining Parts 1 and 2 and by using Lemma 3 we get

since $\sum_{j = 1}^{J}\frac{1}{[(1-\mu\lambda_f)\lambda]^{j}}$ is a decreasing function of $\lambda$ for $\lambda > 0$.

Remark 2. (speed-up analysis). On each time grid, let $T_e$ and $t_e$ be the CPU time spent for applying the backward Euler method directly to the original PDE and the PDE on a subdomain, respectively. Since the number of subdomains is $N$ and the overlap size is small, it holds that $t_e\approx\frac{T_e}{N}$. Now, assume we have $N_t\times N$ processors at our disposal, where $N_t$ processors are assigned to the parareal method and $N$ processors are assigned to the SWR method. The overall cost of the Para-SWR algorithm for per iteration is given by the cost of applying the fine and the coarse propagators simultaneously on each processor|$(J+1)t_e$, the cost of applying the coarse propagator $\mathcal{G}$ sequentially|$N_tt_e$ and some necessary communication cost|$T_c$; in total

At the current iteration, say the $k$-th iteration, the communication cost $T_c$ spent on each coarse time grid $T_n$ mainly comes from transmission of solutions on fine time grids between adjacent subdomains. For example, the processor $P_{n, i}$|assigned to the $n$-th coarse time grid $T_n$ and the $i$-th subdomain, needs to get fine solutions $V_{T_n, 1, 2, \dots, J}(i-1)$ and $V_{T_n, 1, 2, \dots, J}(i+1)$ from its neighboring processor $P_{n, i-1}$ and $P_{n, i+1}$, respectively. Also, it needs to broadcast the solutions $V_{T_n, 1, 2, \dots, J}(i)$ to its two neighboring processors. Clearly, this process can be done simultaneously in all parallelled time intervals. In the longitudinal direction, this can be done sequentially by two steps: first the $N$ processors simultaneously transmit the fine solutions along the upward direction and then along the downward direction (see Figure 2 for illustration). Let $\delta$ be the unit cost spent to communicate a float number between two processors. Then, the total communication time $T_c$ can be written as

where $m$ is the number of space mesh points on a spatial subdomain. Substituting $T_c$ into (3.20) gives

The cost of applying the fine propagator sequentially and directly to the large circuit is

We therefore obtain a speed-up of using the Para-SWR algorithm as

where $K_{\max}$ is the number of iterations performed to meet the stopping criterion. The condition $T_e\gg 2Jm\delta$ occurs naturally if the space domain of the original PDE is large and the discretization mesh $\Delta x$ is small.

3.2.   Weak scaling

The analysis given above is performed on sufficiently long time domains and therefore $\rho(\mathbb{M})$ is independent of $N_t$, the number of paralleled time intervals. We now focus on proving that $\rho(\mathbb{M})$ is robust with respect to $N$, the number of subdomains. We assume $a > 0$ and consider the case $q = 0$. Let $\Delta t$, $J$, $\Delta x$ and $L$ (the length of the space domain) be fixed constants, and then we claim that the arguments $\lambda_f$ and $\lambda_c$ decrease as $N$ increases. Indeed, since $q = 0$ we have $A_+ = A_- = A$ and the eigenvalues are given by

where $m = \frac{L/\Delta x-2}{N}+1$. From (3.25), we have

Clearly, both $\lambda_f$ and $\lambda_c$ increase as $m$ increases, i.e., $\lambda_f$ and $\lambda_c$ decrease as $N$ increases, since $L = [N\times(m-1)+2]\times\Delta x$ is constant. For the matrix $\mathcal{B}$ defined by (2.7) it is easy to know that the quantity $\|\mathcal{B}\|_2$ is independent of the number $m$ and equal to $\frac{\nu^2}{\Delta x^2}$ since $q = 0$. Hence, we get $\mu = \Delta t\frac{\nu^2}{\Delta x^2}$ for all possible $N$.

Now, for the iteration matrix $\mathbb{M}$ defined by (3.9) the aforementioned analysis implies that, except for $\beta$, all of the other elements decrease as $N$ increases. From (3.25), we also have

where $\Gamma(J) = \underset{z\geq0}{\max}\left|\frac{1}{1+Jz}-\frac{1}{(1+z)^J}\right|$. The function $\Gamma(J)$ is shown in Figure 3, and it holds that $\Gamma(J)\leq0.205$ for $J\in\{2, 3, \dots, 500\}$. We therefore can bound $\beta$ as $\beta\leq \frac{\Gamma(J)}{1-\lambda_c}$ for all possible $N$.

Corollary 1 (weak scaling). Let $q = 0$ and $a > 0$. Then, for a fixed $\Delta t$, $J$, $\Delta x$ and $L$, the Para-SWR algorithm converges linearly on sufficiently long time intervals and the convergence rate can be bounded by $\rho(\mathbb{M}^*)$, which is independent of the number of subdomains, provided Condition (3.11) is satisfied with

where the matrix $\mathbb{M}^*$ is derived from $\mathbb{M}$ by letting $\lambda_f = \lambda_f^*$, $\lambda_c = \lambda_c^*$, $\eta = \eta^*$, $\beta = \beta^*$ and $\mu = \mu^*$, and $m_{\max} = 1+\left[\frac{L}{2\Delta x}\right]$ and $\left[\frac{L}{2\Delta x}\right]$ denote the integer part of $\frac{L}{2\Delta x}$.

Proof. Since $N\geq2$ and the length of the space domain $L = [N\times(m-1)+2]\Delta x$ is fixed, it is easy to know that the maximal integer $m$ in (3.26) is $m_{\max} = 1+\left[\frac{L}{2\Delta x}\right]$. Hence, for any possible integer $N$ it holds that $\mathbb{M}\leq\mathbb{M}^*$. Since both $\mathbb{M}$ and $\mathbb{M}^*$ are nonnegative matrices, we know from [45,Theorem 2.21,p. 51] that $\rho(\mathbb{M})\leq\rho(\mathbb{M}^*)$.

Note that, the choice $q = 0$ corresponds to $p = \infty$ and in this case the transmission condition used in (2.3) reduces to the Dirichlet transmission condition. For the general Robin transmission condition, i.e., $q\in(0, +\infty)$, it is difficult to prove a similar conclusion like Corollary 1.

We next investigate how the convergence rate of the Para-SWR algorithm behaves with respect to the quantities $\frac{\nu^2}{\Delta x^2}$ and $a$. It is difficult to make a theoretical analysis, and we have to rely on numerical investigation. Since $\rho(\mathbb{M}^*)$ is a suitable upper bound of $\rho(\mathbb{M})$ and independent of the number of subdomains, we plotted it in Figure 4 on the top of the quantity $\rho(\mathbb{M}^*)$ as a function of $\frac{\nu^2}{\Delta x^2}$ and $a$, where $J = 50$ and $\Delta T = 0.1$ are considered. We see that $\rho(\mathbb{M}^*)$ decreases as $a$ increases or $\frac{\nu^2}{\Delta x^2}$ decreases. This provides us with two possible acceleration strategies for the Para-SWR algorithm. The first one is the change of variables that we have described at the beginning of Section 2. For the model problem (2.1), by letting $u(x, t) = v(x, t)e^{\tau t}$ we get

and this is essentially the same PDE as (2.1), but with a new reaction parameter $a+\tau$. Due to the plot shown in Figure 4 on the top, it seems that the Para-SWR algorithm for (3.29) converges very fast if we choose a very large $\tau$. However, this is a false prediction. Briefly speaking, for a large $\tau$, we find numerically that even though $\underset{j, \; n}{\max}\left\|V^k(x_j, T_n)-V(x_j, T_n)\right\|_2$ is already very small, $\underset{j, \; n}{\max}\left\|U^k(x_j, T_n)-U(x_j, T_n)\right\|_2$ is still a huge number, where $U^k(x_j, T_n)$ is calculated by $U^k(x_j, T_n) = e^{\tau T_n}V^k(x_j, T_n)$; we denote $V^k$ as the numerical solutions generated by the Para-SWR algorithm after $k$ iterations for the transformed PDE (3.29). To illustrate this, by using the example given in (4.1) with $a = 1, \nu = 0.01$, $T = 5$, $\Delta T = 0.1, J = 50$ and $N = 8$, we plotted the diminution history of $\underset{j, \; n}{\max}\left\|V^k(x_j, T_n)-V(x_j, T_n)\right\|_2$ and $\underset{j, \; n}{\max}\left\|U^k(x_j, T_n)-U(x_j, T_n)\right\|_2$ in Figure 4 on the bottom for two choices of $\tau$.

The second possibility is the change of space coordinates. Note that, if we let $x = \omega\tilde{x}$ the PDE problem (2.1) is transformed to

where $\tilde{L} = \frac{L}{\omega}$. The new diffusion coefficient $\omega\nu$ can be sufficiently small by choosing a small positive $\omega$, which implies that the Para-SWR algorithm converges significantly fast. However, in the meantime the computational domain $\left(0, \frac{L}{\omega}\right)$ will become large and to maintain a required spatial accuracy, a uniform step size $\Delta x$ should be used for all $\omega$. Hence, for a fixed $N$|the number of subdomains, the computational cost at each iteration will also significantly increase. To reduce the computational cost at each iteration, we can increase the number of subdomains to $\tilde{N} = N\left[\frac{1}{\omega}\right]$. Since the convergence rate of the Para-SWR algorithm is robust with respect to the number of subdomains, we know that the total computational cost for the Para-SWR algorithm will be evidently reduced by using a small $\omega$. In practical computation, we can choose $\omega = \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \dots$, since in this way it is convenient to increase the subdomains multiplicatively, as $\tilde{N} = 2N, 4N, 8N, \dots$.

4.   Numerical results

In this section, we show several numerical experiments to test the efficiency of the Para-SWR algorithm proposed in this paper. We focus on observing how the convergence rate of this algorithm behaves with respect to $N$|the number of subdomains, $N_t$|the number of paralleled time intervals and $J$|the multiplication ratio between ${\Delta T}$ and ${\Delta t}$. We consider the following PDE problem

where $\nu = 0.01$ and $a = 6$. The space discretization parameter $\Delta x$ is fixed as $\Delta x = \frac{1}{128}$ and the initial solution for the Para-SWR algorithm is chosen randomly in the interval $[0, T]$. For all experiments the iteration stops when the maximal error between the current iteration and the converged solution arrives at a prescribed tolerance.

Example 4.1 (Dirichlet transmission condition). We first consider the Dirichlet transmission condition, i.e., $q = 0$. Let $J = 50$ and $\Delta T = 0.1$, and then it is easy to get $\rho(\mathbb{M}^*) = 0.8465$. Hence, Corollary 1 implies that the convergence rate of the Para-SWR algorithm is robust with respect to $N$, the number of subdomains. The results shown in Figure 5 on the top indicate that this theoretical predication is correct. The other two panels of Figure 5 are concerned with the dependence of the convergence rate of the algorithm on the multiplication ratio $J$ (bottom left) and the number of parallelled time intervals (bottom right). We see that the convergence rate is also robust with respect to $J$ and $N_t$. The results shown in Figure 5 imply that the Para-SWR algorithm really scales weakly, and that this strong robustness indicates that the strategy of change of space coordinates mentioned in Subsection 3.2 is viable.

Example 4.2 (Robin transmission condition). We next consider the Robin transmission condition. For the SWR algorithm, it is well known that the Robin transmission condition is more powerful than the Dirichlet transmission condition; see, e.g., ^[8,9]. For $T = 9$, $J = 50$ and $\Delta T = 0.1$, we show in Figure 6, on the top, the argument $\rho(\mathbb{M})$ as a function of the parameter $q$ $\left( = \frac{1}{\Delta xp-1}\right)$. We see that $\rho(\mathbb{M})$ attains its minimum at $q = 0$, i.e., $p = \infty$, which corresponds to the Dirichlet transmission condition. Does this predict the practical computation well? In Figure 6 on the bottom, we show the errors obtained by running the Para-SWR algorithms with Robin transmission conditions after $k$ iterations and various choices of the free parameters $q$. We see that, indeed, $q = 0$ is the best choice; therefore, these numerical results confirm our theoretical predication shown on the top.

The plot shown in Figure 6 on the top also predicts that increasing $q$ has a negative effect on the convergence rate of the Para-SWR algorithm, and that, for a very large $q$ the algorithm will not converge. This is verified in Figures 7 and 8. In Figure 7, we show on the top row the third, the fourth and the fifth iterations of the 8-subdomain Para-SWR algorithm with $q = 0$, and below we show the same iterations for the algorithm with $q = 4$. This clearly shows that the choice $q = 0$ is more advisable than $q > 0$, which confirms the prediction shown in Figure 6 on the top.

In Figure 8, we show the measured convergence rate of the Para-SWR algorithm in the case of two subdomains and 32 subdomains. The results shown in this figure imply that a large $q$ really slows down the convergence rate; note how the iteration stagnates for $q = 100$.

5.   Conclusions

In this paper, we have presented a hybrid algorithm based on the well understood parareal algorithm and the SWR algorithm. The algorithm, namely, Para-SWR, has been defined by first applying the SWR algorithm to time-dependent PDEs in the case of general $N$ subdomains and then applying the one-iteration parareal algorithm to each subproblem. Convergence of the algorithm was analyzed, and it has been shown theoretically and numerically that the convergence rate is robust with respect to the number of subdomains and the number of parallelled time intervals. This implies that, compared to the parareal algorithm and the SWR algorithm, the proposed hybrid algorithm really brings more parallelism.

There are still some important problems that need to be answered in further work. First of all, the spectral radius $\rho(\mathbb{M})$ is not a good bound of the practical convergence rate. It is conservative and even invalid in some cases. Briefly speaking, for the case $\frac{\nu^2}{\Delta x^2}/a\gg1$ we found that $\rho(\mathbb{M}) > 1$, while the Para-SWR algorithm still converges. For example, by letting $a = 0$ in (4.1) (i.e., the heat equation), $\nu = 0.01$, $T = 10, \Delta x = \frac{1}{128}$, $\Delta T = 0.1$ and $J = 50$, we have $\rho(\mathbb{M}) > 1$ for all $q\geq0$ and $\frac{\nu^2}{\Delta x^2} > 0$, but the numerical results indicate that the Para-SWR algorithm is actually convergent.

The other problem is how to generalize current work to other numerical methods, such as the Runge-Kutta (RK) methods. Such a generalization is not straightforward: if we choose an $s$-stage RK method with $\tilde{s}$ off-steps $\left(\text{the intermediary nodes located in the fine time slice} \left(T_{n+\frac{j\Delta t}{J}}, T_{n+\frac{(j+1)\Delta t}{J}}\right)\right)$, the iteration matrix $\mathbb{M}$ would be a $(1+\tilde{s}+J)\times(1+\tilde{s}+J)$ matrix and its structure would be completely different from the one given by (3.9). Therefore, the current approach for proving $\rho(\mathbb{M}) < 1$ is not applicable. However, one can check that the proof can be straightforwardly generalized to the linear-$\theta$ method, which can be regarded as a family of simple RK methods with $\tilde{s} = 0$.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (12101089) and the Natural Science Foundation of Sichuan Province (2022NSFSC1844).

Conflict of interest

The authors declare that there is no conflicts of interest.