High dimensional Riesz space distributed-order advection-dispersion equations with ADI scheme in compression format

1.   Introduction

In recent decades, fractional differential equations (FDEs) have been widely considered in modelling anomalous diffusion. A class of FDEs, the distributed-order FDEs (DO-FDEs), have been applied in modelling ultraslow diffusion ^[1,2], mixture of delay sources ^[3], and dielectric induction and diffusion ^[4].

To solve DO-FDEs with numerical method, one may discretize the integral with certain quadrature rule, followed by applying an approximation for ordinary fractional derivative ^[5,6,7]. In particular, in ^[8], a form of DO-FDE called Riesz space distributed-order advection-diffusion equations (RSDO-ADE) is studied, with second order schemes proposed for one-dimensional and two-dimensional cases.

On the other hand, a data compression format called the Tensor-Train format (TT-format) ^[9] has been employed in various high dimensional problems. In particular, a TT-format iterative method called the TT-GMRES method is widely used in solving FDE related problems, especially for those which discretized linear systems possess low-rank structure ^[10,11,12]. Benefiting from the Toeplitz-like structure of the linear systems, many studies have been made to explore the underlying properties and to the development of preconditioners ^[13,14]. Another TT-format approach for solving FDE problems is the application of the alternating direction implicit (ADI) method ^[15,16], which reduces the linear systems into one-dimensional systems. Although this approach probably limits the classes of problems to be solved, the convergence analysis is relatively simple to perform. These results provide practical ground for the implementation of RSDO-ADE with higher dimensions, and suggest that the implementation is a trial for the latter approach.

In this work, the RSDO-ADE from ^[8] is considered in high dimensional form. As the concerning fractional orders are distributed over $(0, 1)\cup(1, 2)$ instead of subintervals of $(1, 2)$, and the derivatives are of Riesz type instead of weighted two-sided type, we consider the finite volume approximation in ^[16] not more suitable here than the approximation in ^[8]. Following the proposed discretization, we adopt the midpoint quadrature rule for the integrals, a second order approximation for the Riesz space fractional derivatives, and the Crank-Nicolson method, thus obtaining a second order scheme. Further, similar to ^[15] and ^[16], we apply the ADI method for the reduction of the dimensionality, and TT-format is adopted so that the resulting scheme can be solved in compressed form, provided that the given data possess low TT-ranks.

As compression format introduces perturbations, error analysis is performed to estimate and maintain the overall convergence order. For efficiency, consider a $d$-dimensional case, discretized as linear systems of size $N$ with $N$ time steps, and suppose the numerical solutions possess TT-ranks around $r$, then the proposed method requires storage of ${\mathcal{O}}(2N^{2}+dr^{2}N)$ and operation cost of ${\mathcal{O}}(dr^{2}N^{3}+dr^{3}N^{2})$, provided that the Gaussian Elimination (GE) method is adopted for solving the linear systems. The analysis is testified by some numerical examples with $d\leq20$.

The content is briefly described as follows. The $d$-dimensional RSDO-ADE problem is first presented in Section 2, then the implementation of the TT-format method is described in Section 3, and some numerical results are presented in Section 4.

2.   The $d$-dimensional RSDO-ADE

2.1.   The problem

For ${\bf{x}} = [x^{(1)}, \ldots, x^{(d)}]^{\intercal}$ and $\Omega = \prod\limits_{k = 1}^{d}[x_{a}^{(k)}, x_{b}^{(k)}]$, we consider the $d$-dimensional RSDO-ADE of the following form.

where

Here, $\partial^{\alpha}u/\partial|x^{(k)}|^{\alpha}$ and $\partial^{\beta}u/\partial|x^{(k)}|^{\beta}$ denote the $\alpha$-order and $\beta$-order Riesz space fractional derivatives of $u(x, t)$ respectively, with the forms:

where

where $\Gamma(\cdot)$ denotes the gamma function.

We remark that in this problem form, the boundary condition is modified as described in ^[17].

2.2.   Discretization

The discretization of Problem 2.1 is described below.

For the integrals, by taking a positive integer $S$, we define $\sigma = 1/S$, grid points $\xi_{j} = j\sigma$ for $0\leq j\leq S$, $\alpha_{j} = (\xi_{j}+\xi_{j-1})/2 = (j-0.5)\sigma$ for $1\leq j\leq S$, $\beta_{j} = ((1+\xi_{j})+(1+\xi_{j-1}))/2 = 1+(j-0.5)\sigma$ for $1\leq j\leq S$. Then we define the midpoint quadrature operators

where $v$ can be a function or a matrix. Like the integrals, the operators are linear.

With the operators, we have

For the derivatives, take positive integers $N, M$, and grid points $x_{i}^{(k)} = x_{a}^{(k)}+ih_{k}$ and $t_{m} = m\tau$, where $h_{k} = (x_{b}^{(k)}-x_{a}^{(k)})/(N+1)$, $\tau = T/M$, $0\leq i_{k}\leq N+1$, $1\leq k\leq d$, and $0\leq m\leq M$. Then, denote $u_{i_{1}, \ldots, i_{d}}^{(m)} = u(x_{i_{1}}^{(1)}, \ldots, x_{i_{d}}^{(d)}, t_{m}), \; f_{i_{1}, \ldots, i_{d}}^{(m)} = f(x_{i_{1}}^{(1)}, \ldots, x_{i_{d}}^{(d)}, t_{m})$.

As presented in ^[8,18], the second-order approximations for the Riesz space fractional derivatives have the forms

where

With the notations, for $k = 1, \ldots, d$, define the operators

With the approximations and operators, the Crank-Nicolson method is applied to obtain the second-order implicit finite difference scheme

By complementing appropriate cross terms of ${\mathcal{O}}(\tau^{2})$, we obtain the ADI scheme

2.3.   Matrix form

Denote

$I$ as the $N\times N$ identity matrix, and $\hat{A}$ as the Toeplitz matrix

where $\tilde{g}_{k} = I_{S}^{P}\left(\dfrac{1}{h^{\alpha}}g_{k}^{(\alpha)}\right)+I_{S}^{Q}\left(\dfrac{1}{h^{\beta}}g_{k}^{(\beta)}\right) = \tilde{g}_{-k}$, so $\hat{A}$ is symmetric.

Also, denote

Together with the extended zero boundary condition, the scheme (2.2) can be rewritten as the matrix form

while the ADI scheme (2.3) can be rewritten as

or equivalently,

where $\mathring{A} = \underbrace{A\otimes\cdots\otimes A}_{d}$ and $\mathring{A'} = \underbrace{A'\otimes\cdots\otimes A'}_{d}$, with $A = I+(\tau/2)\hat{A}$ and $A' = I-(\tau/2)\hat{A}$.

In ^[8], it is proved that $\rho(A^{-1}A') = \lVert A^{-1}A'\rVert_{2} < 1$. With this property, we have

such that

3.   Solving process

3.1.   TT-format method

Following ^[15] and ^[16], the linear system (2.4) can be solved in compressed form if ${\bf{u}}^{(0)}$ and ${\bf{f}}^{(m-1/2)}$ are in TT-format, or $u_{0}({\bf{x}})$ and $f({\bf{x}}, t_{m-1/2})$ possess functional TT-format (FTT-format).

Suppose $u_{0}({\bf{x}})$ and $f({\bf{x}}, t_{m-1/2})$ have FTT-formats

where $G_{k}$ and $H_{k}$ are of sizes $r_{k-1}\times r_{k}$ and $r'_{k-1}\times r'_{k}$ respectively for $1\leq k\leq d$, then ${\bf{u}}^{(0)}$ and ${\bf{f}}^{(m-1/2)}$ can be reshaped to tensors ${\mathscr{U}}^{(0)}$ and ${\mathscr{F}}^{(m-1/2)}$ with TT-formats

where $U_{k}^{(0)}$ and $F_{k}^{(m-1/2)}$ are of sizes $r_{k-1}\times N\times r_{k}$ and $r'_{k-1}\times N\times r'_{k}$ respectively for $1\leq k\leq d$.

After this, the numerical solutions can be obtained in TT-format through the following algorithm.

In the algorithm, the notation $U\times_{2}A$ denotes the multiplication of the matrix $A$ with $U$ along the second dimension, or with the mode-$2$ fibers of $U$; while "Round" denotes a recompression process called rounding for TT-format, expecting to reduce TT-ranks by introducing relative error.

3.2.   Computational efficiency

Here, we discuss the efficiency of Algorithm 3.1.

Suppose the GE method is adopted in line 1 to compute ${\mathcal{S}}$, and the numerical solutions possess TT-ranks $\approx r$.

For storage, as the major components are the matrices and TT-cores of the numerical solutions, the GE method requires storage of ${\mathcal{O}}(N^{2}+dr^{2}N)$.

For operation cost, there are three major parts as shown in Algorithm 3.1. The first part is line 1, the matrix inversion; the second part is line 6 to line 10, matrix-vector multiplications in single time step; and the third part is line 11, rounding process in single time step. Thus the GE method requires totally ${\mathcal{O}}(N^{3}+M(dr^{2}N^{2}+dr^{3}N))$ operations.

We can see that the TT-ranks $r$ is a crucial factor affecting the efficiency. Theoretically, $r$ may range from small numbers to powers of $N$, so that compressed format method may be less efficient than full storage method. With rough estimation, if the problem data possess TT-ranks $\approx r'$, then the numerical solutions can be computed with TT-ranks of at most ${\mathcal{O}}(Mr')$, which is attained when rounding process does not compress at all. Numerical examples indicate a tentative form of growth that $r\approx {\mathcal{O}}(\log^{k}N)$ for some $k > 0$, which implies a possible dominance of the rounding process over the whole solving process. Further studies about the possible forms of growth of TT-ranks are to be made.

Remark 1. By virtue of the Toeplitz structure of $\hat{A}$, hence $I+(\tau/2)\hat{A}$, the operation cost may be reduced to ${\mathcal{O}}(N\log N+M(dr^{2}N\log N+dr^{3}N))$ in the optimal case. This may be accomplished by adopting the circulant-and-skew-circulant representation method ^[19], equipped with preconditioned conjugate gradient method and the Fast Fourier Transform ^[20].

3.3.   Error analysis

For the overall convergence of the proposed method, we have the following proposition.

Proposition 1. Suppose $P(\alpha), \; Q(\beta)$, and the exact solution $u({\bf{x}}, t)$ of Problem 2.1 satisfy the following conditions.

i. $\dfrac{\partial^{2}P(\alpha)}{\partial\alpha^{2}}, \; \dfrac{\partial^{2}}{\partial\alpha^{2}}\dfrac{\partial^{\alpha}u({\bf{x}}, t)}{\partial|x^{(k)}|^{\alpha}}\in {\mathcal{C}}(\Omega\times[0, T]\times[0, 1])$ with respect to $({\bf{x}}, t, \alpha)$, and $\dfrac{\partial^{2}Q(\beta)}{\partial\beta^{2}}, \; \dfrac{\partial^{2}}{\partial\beta^{2}}\dfrac{\partial^{\beta}u({\bf{x}}, t)}{\partial|x^{(k)}|^{\beta}}\in {\mathcal{C}}(\Omega\times[0, T]\times[1, 2])$ with respect to $({\bf{x}}, t, \beta)$.

ii. For some $\rho > 0$, for $1\leq k\leq d$, the mixed partial derivatives

are in ${\mathcal{C}}(\Omega\times[0, T])\cap {\mathcal{L}}^{1}(\Omega\times[0, T])$ for all $\gamma_{k}\in[0, 5+\rho], \; \gamma_{*}\in[0, 3+\rho]$, and vanish at infinity for all $\gamma_{k}\in[0, 4+\rho], \; \gamma_{*}\in[0, 2+\rho]$.

Then for rounding relative error $\epsilon$, we have

where $E_{m}$ is the discrete 2-norm error between the exact solution in Eq (2.4) and the perturbed numerical solution at $t_{m}$ due to rounding, $K > 0$ is a constant independent of $m, \tau, h_{k}, \sigma$, and $\hat{h} = \max\limits_{1\leq k\leq d}h_{k}$.

Proof. Denote ${\bf{u}}^{(m)}$ as the exact solution in Problem 2.1, $\hat{{\bf{u}}}^{(m)}$ as the scheme solution in Eq (2.4), $\tilde{{\bf{u}}}^{(m)}$ as an intermediate solution, and $\acute{{\bf{u}}}^{(m)}$ as the perturbed numerical solution. Then they have the following relations.

$\bullet$ By similar arguments in ^[8] and ^[16], from the given conditions, we have $(h_{1}\cdots h_{d})^{1/2}\|{\bf{u}}^{(m)}\|_{2}\leq B$, and

for some constants $B, K_{0} > 0$ independent of $m, \tau, h_{k}, \sigma$.

● $\hat{{\bf{u}}}^{(0)} = \acute{{\bf{u}}}^{(0)}$.

● $\mathring{A}\hat{{\bf{u}}}^{(m)} = \mathring{A'}\hat{{\bf{u}}}^{(m-1)}+\tau{\bf{f}}^{(m-\frac{1}{2})}$.

● $\tilde{{\bf{u}}}^{(m)} = \mathring{A}^{-1}\mathring{A'}\acute{{\bf{u}}}^{(m-1)}+\tau\mathring{A}^{-1}{\bf{f}}^{(m-\frac{1}{2})}$.

● $\|\tilde{{\bf{u}}}^{(m)}-\acute{{\bf{u}}}^{(m)}\|_{F}\leq\epsilon\|\tilde{{\bf{u}}}^{(m)}\|_{F}$.

By similar arguments in ^[15], combining the relations with the norm property in inequality (Eq 2.5), we have

where $K = \max\{2K_{0}T, B\} > 0$ is independent of $m, \tau, h_{k}, \sigma$.

This allows us to take $\epsilon = 1/M^{3}$ such that the overall convergence is of second order.

Remark 2. In ^[15], the rounding perturbation analysis for a non-Crank-Nicolson scheme is demonstrated; while in ^[16], the ADI perturbation for a Crank-Nicolson scheme is analyzed, with the corresponding rounding perturbation analysis omitted. The above proof is supplemented for the rounding perturbation analysis for a Crank-Nicolson scheme.

4.   Results and discussion

In this section, some numerical examples are presented to test the methods described above. For simplicity, we set $[x_{a}^{(k)}, x_{b}^{(k)}] = [0, 1]$ for $1\leq k\leq d$, $[0, T] = [0, 1]$, and $S = M = N$. We will show the numerical errors (Error), corresponding convergence rates (Rate), CPU times in the whole process (CPU(s)), CPU times in rounding process (rCPU(s)), and $($mean, mode, maximum$)$ of the involving TT-ranks ($M_{r}$).

The numerical examples are tested in MATLAB R2014b with the configuration: Intel(R) Core(TM) i5-8300H CPU 2.30 GHz and 8 GB RAM.

Example 1. This example is tested with a TT-ranks $1$ exact solution, with

where

The performance is summarized in Table 1, with "Error'' column showing $E_{M}$ at $t_{M} = 1$.

Example 2. This example is tested with an exact solution of TT-ranks $3$, with

where

for $2\leq k\leq d-1$, $y^{[k, \ell]} = 10(x^{(k)})^{\ell}(1-x^{(k)})^{\ell}$ for $1\leq k\leq d, \; 3\leq\ell\leq 5$, and

which is defined entrywisely.

As the closed form of $G_{k}'(x^{(k)})$ is numerically unstable, each entry of $G_{k}'(x^{(k)})$ is first reduced to the following form, then the remaining integrals are approximated by the midpoint quadrature rule.

with $\Gamma(9-\gamma)/\Gamma(\ell'-\gamma) = \sum_{r' = 0}^{9-\ell'}c_{\ell', r'}\gamma^{r'}, \; \gamma = \alpha, \beta$.

The performance is summarized in Table 2, with "Error'' column showing $E_{M}$ at $t_{M} = 1$.

Example 3. This example is tested with uncertain exact solution, and the same $u_{0}({\bf{x}})$ and $f({\bf{x}}, t)$ of TT-ranks $3$ as in Example 3 in ^[15], with

The performance is summarized in Table 3, with "Error'' column showing $E_{M}$ at $t_{M} = 1$, assuming the numerical solution of the finest grid as the exact solution.

In the tables, second order convergence is observed, agreeing with the perturbation analysis. In Tables 1 and 2, the $M_{r}$ illustrate that the problems with low TT-ranks exact solutions may have numerical solutions of higher TT-ranks, yet the ranks still grow steadily of about ${\mathcal{O}}(\log N)$, and the rounding process occupy about one-fourth of the whole solving process. While in Table 3, there seems to be a nonlinear growth of ${\mathcal{O}}(\log N)$ in the TT-ranks, and the rounding process occupy about half of the whole solving process. These agree with the discussion about the effect of the TT-ranks in operation cost. Besides, the deviations of TT-ranks of the numerical solutions from the exact solutions also happen in the related works ^[15] and ^[16] where the ADI method is applied, implicating that the ADI treatment is a potential factor for this phenomenon. The growth of TT-ranks is to be furthered studied.

5.   Conclusions

In this paper, high dimensional RSDO-ADEs are discretized with a second order ADI scheme, and TT-format method is introduced to solve the scheme in compressed form. The perturbation error analysis is performed, and it is claimed that under certain conditions, with rounding relative error $\epsilon = 1/N^{3}$, Algorithm 3.1 can maintain the second order convergence of the scheme. Numerical experiments with low TT-ranks data are conducted, the results agree with the claim. Meanwhile, the TT-ranks appear to be indicating a growth of power of $\log N$ or other nonlinear function of $\log N$, hence implying a possible dominance of the rounding process over the whole solving process. Further studies about the TT-ranks of the numerical solutions are to be made.

Acknowledgments

This work was supported by University of Macau [MYRG2020-00208-FST, MYRG2018- 00025-FST], the Science and Technology Development Fund, Macau SAR under Funding Scheme for Postdoctoral Researchers of Higher Education Institutions 2021 (File No. 0032/APD/2021).

Conflict of interest

All authors declare no conflicts of interest in this paper.