A modified domain decomposition spectral collocation method for parabolic partial differential equations

1.   Introduction

The parabolic partial differential equation (PDE)

equipped with the initial condition

and boundary condition

is a class of important mathematical models, where $\Omega = (\overline{a}, \overline{b})\times (\tilde{a}, \tilde{b})$, and $k^*_i, \; i = 1, 2, \cdots, 5$ are all smooth functions of the variables $x, \; y, \; t$. This class of models plays important roles in many engineering problems. For example, when $k_1 = k_2 = 1, \; k_3 = k_4 = k_5 = 0$, and $f^*(u, x, y, t) = g(x, y, t)$, (1.1) becomes the heat conduction equation ^[1]. Taking $k_1 = k_2 = 1, \; k_3 = k_4 = k_5 = 0$, and $f(u, x, y, t) = au(1-u^2)+g(x, y, t)$, (1.1) represents the Allen-Cahn equation ^[2]. When describing the Brownian motion of particles, denote by $y$ the velocity of the particles, and let $k_1$ be a nonzero constant, $k_2 = k_3 = 0, \; k_4 = -y, \; k_5 = \beta y-\gamma x$, and $f(u, x, y, t) = \beta u$. Then, (1.1) corresponds to the Fokker-Planck equation ^[3]. Moreover, some reaction-diffusion models and convection-diffusion models can also be seen as special cases of (1.1) ^[4,5].

Usually, it is difficult to exactly solve Eqs (1.1)–(1.3), and hence finding its numerical solutions is a more popular choice for many researchers. In recent decades, finite difference methods (FDM) and finite element methods (FEM) have been increasingly developed for solving various PDEs, including the parabolic case Eqs (1.1)–(1.3) ^{[6,7,8,9,10,11,12,13,14]}. The most remarkable feature of these two methods is that they can lead to sparse linear equations when some iterative methods (such as Netwon iterative method) are employed to deal with the nonlinear term. However, we all know that these methods are usually inferior to spectral methods in terms of approximating precision.

Spectral methods have attracted extensive attentions since they were first proposed. As an excellent competitor for numerically solving PDEs, they possess great approximation accuracy and are very easy to implement ^{[15,16,17,18]}. Especially, Lui et al. ^[19,20,21] have recently presented spectral discretization in both time and space, and this can greatly accelerate the solving of PDEs containing both space and time factors. However, the drawback of common spectral methods is that the resultant coefficient matrix is dense, which leads to difficulties when solving the corresponding linear equations. This fact is very obvious when spectral methods are used for high-dimensional PDEs or PDEs with large intervals in time and/or space. For remedying the above deficiency, some researchers have considered the combination of domain decomposition with spectral method, and presented some domain decomposition spectral methods ^[3,22,23].

In this paper, we study a novel domain decomposition spectral collocation method employing Legendre polynomials as basis functions. Comparing with the existing techniques, our proposed method enjoys an additional advantage, that is, the resultant linear system has a saddle structure, for which a specially designed preconditioning technique is suitable. In the Sections 3 and 4, we will see that implementing this preconditioner only requires solving the sub-problem $Ax = y$, with $A$ a block diagonal matrix. This means that a lot of CPU time can be saved.

The rest of this paper is as follows. In Section 2, the modified domain decomposition spectral method (shortly, MDDSM) is introduced, and the difference between MDDSM and the conventional domain decomposition spectral method (shortly, CDDSM) is discussed. Later, as an important component, in Section 3, according to ^[24,25], we apply a dimension expanding preconditioning technique to iteratively solve the linear system obtained in Section 2. In Section 4, numerical examples are given to show the efficiency of the proposed MDDSM and the preconditioner. Some conclusions are finally drawn in the last section.

2.   Description of MDDSM

We first divide the whole interval $[0, T]\times\Omega$ into some big uniform sub-domains. For the fixed positive integers $\overline{p}, \; \tilde{p}, \; \hat{p}$, let $\overline{q} = (\overline{b}-\overline{a})/\overline{p}, \; \tilde{q} = (\tilde{b}-\tilde{a})/\tilde{p}, \; \hat{q} = T/\hat{p}$, and we give a uniform partition $[T_0, T_1], \; [T_1, T_2], \; \cdots, \; [T_{\hat{p}-1}, T_{\hat{p}}]$ in time and $\Omega_{i, j} = [\overline{a}_{i-1}, \overline{a}_{i}]\times [\tilde{a}_{j-1}, \tilde{a}_{j}], \; i = 1, \; 2, \; \cdots, \; \overline{p}, \; j = 1, \; 2, \; \cdots, \; \tilde{p}$ in space, where $T_0 = 0, \; T_{\hat{p}} = T$, $T_{k}-T_{k-1} = \hat{q}, \; k = 1, \; 2, \; \cdots, \; \hat{p}$, $\overline{a}_0 = \overline{a}, \; \overline{a}_{\overline{p}} = \overline{b}$, $\tilde{a}_0 = \tilde{a}, \; \tilde{a}_{\tilde{p}} = \tilde{b}$, and $\overline{a}_i-\overline{a}_{i-1} = \overline{q}, \; i = 1, \; 2, \; \cdots, \; \overline{p}$, $\tilde{a}_j-\tilde{a}_{j-1} = \tilde{q}, \; j = 1, \; 2, \; \cdots, \; \tilde{p}$.

Next, at the $k_{th}$ level $[T_{k-1}, T_k], \; k = 1, \; 2, \; \cdots, \; \hat{p}$, and in the big sub-domains $\Omega^k_{i, j} = \Omega_{i, j}\times [T_{k-1}, T_k]$, we do the transformation $x = \frac{\overline{q}}{2}x^*+(i-\frac{1}{2})\overline{q}+\overline{a}, \; y = \frac{\tilde{q}}{2}y^*+(j-\frac{1}{2})\tilde{q}+\tilde{a}, \; t = \frac{\hat{q}}{2}t^*+(k-\frac{1}{2})\hat{q},$ and denote

Then, by some computations, Eq (1.1) is transformed into an equivalent system of the form

where $k_i = k_i(x^*, y^*, t^*), \; i = 1, \; \cdots, \; 5$. For ease of description, we uniformly replace the variables $x^*, \; y^*, \; t^*$ with $x, \; y, \; t$, and rewrite the equivalent equation of (2.1) as

with $k_i = k_i(x, y, t), \; i = 1, \; \cdots, \; 5$.

Similarly, the initial value (1.2) and boundary condition (1.3) can be equivalently transformed into

and

where $\varphi(x, y), \psi(x, y, t)$ are known functions.

Now we start to compute the approximating solution of $v(x, y, t)$ in Eqs (2.2)–(2.4) at the first level $[T_0, T_1]$. In each big sub-domain $\Omega^1_{i, j}, \; i = 1, \; 2, \; \cdots, \; \overline{p}, \; j = 1, \; 2, \; \cdots, \; \tilde{p},$ we choose the Legendre-Gauss-Lobatto type points $x^i_0, \; \cdots, \; x^i_{\overline{N}}, \; y^j_0, \; \cdots, \; y^j_{\tilde{N}}, \; t^1_0, \; \cdots, \; t^1_{\hat{N}}$ as the nonuniform grid nodes in $x, \; y, \; t$ directions, respectively; this can be seen in Page 20–21 in ^[15]. It is easy to understand that $x^i_{\overline{N}} = x^{i+1}_0, \; i = 1, \; 2, \; \cdots, \; \overline{p}-1, \; y^j_{\tilde{N}} = y^{j+1}_0, \; j = 1, \; 2, \; \cdots, \; \tilde{p}-1$.

For conveniently implementing the coming domain decomposition spectral collocation method, we denote by $P_{all}$ the set of all grid nodes $(x^i_l, y^j_m, t^1_n), \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}, \; l = 0, \; \cdots, \; \overline{N}, \; m = 0, \; \cdots, \; \tilde{N}, \; n = 0, \; \cdots, \; \hat{N}$, and we group this set into the following four subsets.

($g_1$) $P_{inner}:$ interior points: $(x^i_l, y^j_m, t^1_n), \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}, \; l = 1, \; \cdots, \; \overline{N}-1, \; m = 1, \; \cdots, \; \tilde{N}-1, \; n = 1, \; \cdots, \; \hat{N}.$

($g_2$) $P_{init}:$ initial values points: $(x^i_l, y^j_m, t^1_0), \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}, \; l = 0, \; \cdots, \; \overline{N}, \; m = 0, \cdots, \; \tilde{N}.$

($g_3$) actual boundary points $P_{act}$ (namely grid points that lie on $\partial \Omega$):

($g_4$) ghost boundary points $P_{gho}$: the interior grid points which lie on $\partial(\Omega^1_{i, j})$, namely $P_{gho} = P_{all}-(P_{inner}\cup P_{init}\cup P_{act})$. Specifically, when $l = 0$ or $l = \overline{N}$, $(x^i_l, y^j_m, t^1_n)$, they are called a $x$-type ghost boundary points (shortly, $x$-type gbp), alternatively, when $m = 0$ or $m = \tilde{N}$, $(x^i_l, y^j_m, t^1_n)$, they are called a $y$-type ghost boundary points (shortly, $y$-type gbp). (Figure 1 shows an illustrative example of the ghost boundary points $P_{gho}$, where the green line represents an $x$-type gbp, and the red shows a $y$-type gbp.)

Moreover, in order to ensure the number of unknowns is the same as that of equations, we should employ the fact that the left-half derivative is equivalent to the right-half derivative, and consequently, we impose the following ghost boundary conditions on each ghost boundary points:

(a) the x-type ghost boundary condition (for $x$-type gbp):

(b) the y-type ghost boundary condition (for $y$-type gbp):

where $v^{gx}_{i, (j-1)\widetilde{N}+m, n}, \; v^{gy}_{(i-1)(\overline{N}-1)+l, j, n}$ are all undetermined unknowns.

Let

be the approximation of the exact solution $v(x, y, t)$ of (2.2)-(2.4), where $L_l(x), \; L_m(y)$ and $L_n(t)$ are all Legendre polynomials with degrees of $l, \; m$ and $n$, respectively, and $v^{i, j, 1}_{l, m, n}$ are the unknowns to be determined. Then, in each sub space $\Omega^1_{i, j}$, taking $x = x^i_{l}$, $y = y^j_m$, $t = t^1_n$ as the collocation points, and substituting (2.7) into Eqs (2.2)–(2.6), we can obtain the corresponding algebraic equations. Concretely, MDDSM can be listed as the following algorithm.

For describing the resultant linear equations, we define all the unknowns as $\mathbf{v}_h = (\mathbf{v}^a_h, \mathbf{v}^g_h)$, with

where $\mathbf{v}^a_{i, j}$ is the solution vector involved in $\Omega^1_{i, j}$, and $\mathbf{v}^g_{i, j}$ is the ghost value vector included in Eqs (2.5) and (2.6) on $\partial\Omega^1_{i, j}$. That is,

and

Then, by deleting the error term caused by the approximating solution function $v_h(x, y, t)$ in (2.7), Algorithm 1 will result in $M = \overline{p}\tilde{p}(\overline{N}+1)(\tilde{N}+1)(\hat{N}+1)+ (\overline{p}-1)(\tilde{p}-1)(\tilde{N}\hat{N}+(\overline{N}-1)\hat{N}) +(\overline{p}-1)(\tilde{N}-1)\hat{N}+(\tilde{p}-1)(\overline{N}-1)$ algebraic equations

When $f(v, x, y, t)$ in Eq (2.2) is a linear term of $v$, it is easy to know that Eq (2.9) becomes the following linear system

where the coefficient matrix $H$ has a saddle structure of the form

Concretely, $A = diag(A_{1, 1}, \; A_{2, 2}, \; \cdots, \; A_{N_1, N_1}), \; N_1 = \overline{p}\tilde{p}$, $size(A_i) = (\overline{N}+1)(\tilde{N}+1)(\hat{N}+1) \times (\overline{N}+1)(\tilde{N}+1)(\hat{N}+1), \; i = 1, \; 2, \; \cdots, \; N_1$, where $size(A)$ expresses the numbers of rows and columns of $A$, and $B$ is a full column rank matrix that has at most one entry per row equal to $-1$, at least two entries per column equal to $-1$, and at most four entries per column equal to $-1$. Specifically, when $k_i, \; i = 1, \; \cdots, \; 5$ in (2.2) are all constants, there is $A_{1, 1} = A_{2, 2} = \cdots = A_{N_1, N_1}$. The left one in Figure 2 shows an illustration of $H$ taking $\overline{a} = \tilde{a} = 0, \overline{b} = \tilde{b} = T = 12, \overline{N} = \tilde{N} = \hat{N} = 8, \overline{q} = \tilde{q} = \hat{q} = 2$.

When $f(v, x, y, t)$ in Eq (2.2) is a nonlinear term of $v$, then, after dealing with the nonlinear algebraic equations by employing some iterative methods (for instance, the classic Newton iterative methods), (2.9) will still result in the saddle linear system (2.10).

Considering we have computed the values of $v^{i, j, k-1}_{l, m, n}, \; l = 0, \; \cdots, \; \overline{N}, \; m = 0, \; \cdots, \; \tilde{N}, \; n = 0, \; \cdots, \; \hat{N}, \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}$, we now substitute the points $(x^i_l, y^j_m, t^{k-1}_{\hat{N}}), \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}, \; l = 0, \; \cdots, \; \overline{N}, \; m = 0, \; \cdots, \; \tilde{N}$ into

to get $v_h(x^i_l, y^j_m, t^{k-1}_{\hat{N}})$. These are also the values of $v_h(x^i_l, y^j_m, t^{k}_{0})$ since $(x^i_l, y^j_m, t^{k-1}_{\hat{N}}) = (x^i_l, y^j_m, t^{k}_0), \; i = 1, \; \cdots, \; \overline{p}, \; j = 1, \; \cdots, \; \tilde{p}, \; l = 0, \; \cdots, \; \overline{N}, \; m = 0, \; \cdots, \; \tilde{N}$. Consequently, in this way, we obtain the initial values at the $k_{th}$ level $[T_{k-1}, T_k]$. By virtue of this group of initial values, we can solve

using the same idea of Algorithm 1.

For investigating the difference between MDDSM and CDDSM, here we additionally give the idea of CDDSM and the resultant algebraic equations. The actual efficiency of these two algorithms will be shown in the numerical examples in Section 4.

In the CDDSM, the ghost boundary conditions are replaced by

(c) the x-type ghost boundary condition (for $x$-type gbp):

(d) the y-type ghost boundary condition (for $y$-type gbp):

Then, CDDSM can be expressed as follows

From these two algorithms, we can clearly find the unique difference between MDDSM and CDDSM is the ghost boundary conditions in step 4. In MDDSM, the participation of ghost unknowns $v^{gx}_{i, (j-1)\widetilde{N}+m, n}, \; v^{gy}_{(i-1)(\overline{N}-1)+l, j, n}$ leads to a saddle structure of the coefficient matrix $H$, for which we will see that a dimension expanding preconditioner is suitable in the next section. However, in CDDSM, the resulting linear system is

where $\mathbf{v}_h^a$ is the unknown vector, and the coefficient matrix $\hat{H}$ is no longer a saddle matrix. The right figure in Figure 2 shows an illustration of $\hat{H}$ corresponding to $\overline{a} = \tilde{a} = 0, \; \overline{b} = \tilde{b} = T = 12, \; \overline{N} = \tilde{N} = \hat{N} = 8, \; \overline{q} = \tilde{q} = \hat{q} = 2$.

3.   The application of dimension expanding preconditioning technique on the saddle problem of MDDSM

In this section, we introduce a preconditioning technique for quickly solving the saddle system (2.10). This kind of preconditioning technique was first presented by Luo et al. ^[24,25]. When this technique is used combined with some Krylov subspace method such as the GRMES method ^[26], an advantage is that, in each iteration, one only needs to solve a sub problem $Ax = r,$ where $A$ is a block diagonal matrix.

For ease of description, we rewrite (2.10) as

and suppose the size of $A, \; B, \; C$ are respectively $m\times m$, $m\times n$, and $n\times m$. Clearly, (3.1) can be equivalently augmented as

where

with $I$ an $n\times n$ identity matrix, and $\mathbf{u}_2 = -\mathbf{u}_3$.

Denoting by $\tilde{A} = A+\frac{\tilde{\alpha}}{\tilde{\alpha}-1}BC$, and letting $\tilde{\alpha}$ be a chosen parameter, we take the preconditioner as

This preconditioner $P_{DE}$ can be implemented by

where

and $Q_3^{-1}$ can be further expressed as

Remark 1. From Eq (3.5) we see that when implementing the preconditioner $P_{DE}$ in Eq (3.3), some parallel computations can be employed. In fact, when computing $P_{DE}^{-1}r$, one of important step is to solve $A^{-1}v$, and this can be independently completed by multiple processors, considering the fact that $A$ is a block diagonal matrix.

4.   Numerical results

In this section, we give some examples to investigate the actual efficiency of the presented MDDSM and the preconditioner $P_{DE}$. For reflecting the superiority, we compare MDDSM with CDDSM from the perspectives of approximating accuracy and the costed CPU time. When dealing with the nonlinear PDEs (Example 3), we transform the resultant nonlinear algebraic equations into linear equations by the classical Newton iterative method, always taking an initial value vector $x = ones(N)$. All the obtained linear algebraic equations are solved by the popular GMRES method ^[26]. In MDDSM, the GMRES method is equipped with the preconditioner $P_{DE}$, and in CDDSM, the common ILU preconditioner is used. All these examples are implemented using MATLAB R2016(a) on a PC equipped with an Intel(R)Core(TM)i5-8265U processor (CPU@1.60GHz).

The following Table 1 is given only to investigate the resultant accuracy of the presented MDDSM in both time and space, hence, the convergence tolerance for GMRES (in the left column) is taken as $1e-15$, and the maximum number of iterations is set equal to 1000. Furthermore, in this table, we additionally investigate the efficiency of the common ILU preconditioner (setup.type = 'crout'; setup.milu = 'row'; setup.droptol = 0.1) for MDDSM with convergence tolerance $1e-6$, and the results are reported in the right column.

In the Tables 2–8, 'MDDSM+DE' means that the PDE is solved by the MDDSM method, and the obtained linear system is solved by the GMRES method equipped with DE preconditioner, 'CDDSM+ILU' means that the PDE is solved by the CDDSM method, and the obtained linear system is solved by the GMRES method equipped with ILU preconditioner (setup.type = 'crout'; setup.milu = 'row'; setup.droptol = 0.1), and 'CDDSM+direct' means that the PDE is solved by the CDDSM method, and the obtained linear system is solved by the direct method $x = H\backslash b$. This direct method is listed only to check the accuracy of CDDSM. In all the runs in Tables 2–8, the convergence tolerance for GMRES is taken as $1e-6$, and the maximum number of iterations is still set equal to 1000. 'size' shows the number of rows of the coefficient matrix $H$, 'iter' denotes the number of iterations in the GMRES method, 'time1(s)' and 'time2(s)' respectively denote the time for discretizing the PDE (2.2)–(2.4) using MDDSM/ CDDSM and the time for solving the resultant linear system (2.10) (using $P_{DE}$), (2.13) (using ILU preconditioner), or (2.13) (using the direct method). '*' indicates that the GMRES method fails. 'error' is the error between the exact solution and the approximating solution at all the interior mesh points, initial value mesh points, and actual boundary mesh points, and this is computed using the infinite norm of a vector.

Moreover, in Tables 7 and 8, 'iter(Newton)' is the number of iterations in the Newton iterative method, and 'iter(gmres)' is the sum of all the iterations of GMRES involved in the whole Newton iterative method.

For convenience, by taking $\overline{a} = \tilde{a}, \overline{b} = \tilde{b}$ in the space domain $\Omega$, we consider the following examples.

Example 1. (the heat conduction equation):

$k_1^* = k_2^* = c, \; k_3^* = k_4^* = k^*_5 = 0$, $f(u, x, y, t) = cos(x+y+t)+2c sin(x+y+t)$, $u(x, y, t) = sin(x+y+t)+2.$

Example 2. (with variable coefficients):

$k_1^* = 0.5exp(x+2y+5t), \; k_2^* = cos(x-y+2t)/exp(x+y+2t)$, $k_3^* = 0.1, \; k_4^* = 0.2, \; k^*_5 = -0.1$, $f(u, x, y, t) = 20u+g(x, y, t), \; u(x, y, t) = sin(x+y+t)+2.$

Example 3. (the Allen-Cahn equation):

$k_1^* = k_2^* = 1, \; k_3^* = k_4^* = k^*_5 = 0, \; f(u, x, y, t) = \frac{1}{\varepsilon}u(1-u^2)+g(x, y, t)$, $g(x, y, t) = cos(x+y+t)+2sin(x+y+t)-\frac{1}{\varepsilon} (sin(x+y+t)+2-(sin(x+y+t)+2)^3)$, $u(x, y, t) = sin(x+y+t)+2.$

From the numerical results in Tables 1–8, we can get the following observations:

1. From the left column in Table 1, we see that for a fixed interval $[0, T]\times \Omega$, as the degree of the Legendre polynomial (also the number of collocation points) increases, the resulting error decreases accordingly. Moreover, the results in the right column of Table 1 reflect that the common ILU preconditioner is ineffective for MDDSM.

2. In all the tested intervals in time and space, the expected accuracy $1e-6$ can be obtained by only taking $\overline{N} = \tilde{N} = \hat{N} = 7$ or $8$ in MDDSM. This means that MDDSM does well in large time-space interval $\Omega\times [0, T]$ with little cost.

3. The results in Tables 4, 5, and 8 reflect that MDDSM is stable in time level, namely, supposing we have obtained the values in $[T_0, T_1], \; [T_1, T_2], \; \cdots, \; [T_{k-2}, T_{k-1}]$, then, the initial value of the $k_{th}$ level $[T_{k-1}, T_{k}]$ can be chosen from the value of $T_{k-1}$ computed in $[T_{k-2}, T_{k-1}]$, and the error in $[T_{k-1}, T_{k}]$ will be of the same order of magnitude of $[T_{k-2}, T_{k-1}]$.

4. Tables 2–8 show that the DE preconditioner is efficient in convergence, and that it is stable, namely, it appears to be independent of the size of interval in time/space.

5. From Tables 7 and 8, we can see that the classical Newton iterative method is suitable for the nonlinear algebraic equations which are obtained by using MDDSM to solve the tested model.

6. In Tables 2 and 3, we find that ''MDDSM+DE" is obviously superior to ''CDDSM+ILU"; in fact, in the heat equation, ILU has failed in CDDSM.

7. In all these results, we find that, for a fixed $N$, the number of iterations of the GMRES method is basically independent of the time interval $[0, T]$ and space region $\Omega$, as reflected in Tables 2–8. Consequently, for a given numerical experiment, we can choose the parameter $\tilde{\alpha}$ by testing it in a small region $[0, T]\times \Omega$.

5.   Conclusions

In this paper, we presented a modified domain decomposition spectral collocation method, and explored a dimension expanding preconditioner for the obtained linear system. Numerical examples show that this MDDSM is efficient for model (1.1) with proper boundary values and initial values in a large time-space interval $\Omega\times [0, T]$, and that the classical Newton iterative method and DE preconditioner are suitable for the proposed MDDSM.

How to choose an optimal parameter $\tilde{\alpha}$ in (3.3) is still a concern of our future work. Also, we will focus on improving this kind of dimension expanding preconditioner for saddle problem (2.10).

Author contributions

Wei Hua Luo conceived the idea for this paper and derived the algebraic equations; Liang Yin conducted all the numerical experiments; and Jun Guo was responsible for editing the paper.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is supported by the Scientific Research Fund of Hunan Provincial Science and Technology Department (2022JJ30416), the Scientific Research Fund of Hunan Provincial Education Department (22A0483). The work of Jun Guo is supported by the Sichuan National Applied Mathematics co-construction project (2022ZX004), CUIT (KYTD202243) and the Scientific Research Foundation (KYTZ202184).

Conflict of interest

The authors declare there is no conflict of interest.