Krylov Integration Factor Method on Sparse Grids for High Spatial Dimension Convection–Diffusion Equations

Abstract

Krylov implicit integration factor (IIF) methods were developed in Chen and Zhang (J Comput Phys 230:4336–4352, 2011) for solving stiff reaction–diffusion equations on high dimensional unstructured meshes. The methods were further extended to solve stiff advection–diffusion–reaction equations in Jiang and Zhang (J Comput Phys 253:368–388, 2013). Recently we studied the computational power of Krylov subspace approximations on dealing with high dimensional problems. It was shown that the Krylov integration factor methods have linear computational complexity and are especially efficient for high dimensional convection–diffusion problems with anisotropic diffusions. In this paper, we combine the Krylov integration factor methods with sparse grid combination techniques and solve high spatial dimension convection–diffusion equations such as Fokker–Planck equations on sparse grids. Numerical examples are presented to show that significant computational times are saved by applying the Krylov integration factor methods on sparse grids.

Appendix: Linear Stability Analysis of the IIF2 Scheme (7) for CDR Equations

To analyze the linear stability of IIF schemes, we use the following scalar linear test equation

$$\begin{aligned} u_t=au-du+ru, \qquad \text{ with } r\in \mathcal {C}, \text{ and } a, d \in \mathcal {R}, d>0. \end{aligned}$$

(23)

In the context of solving CDR equations, a and d actually represent spatial discretizations for the convection term and the diffusion term respectively. Following the stability analysis approach in [22], we show boundaries of the stability regions in the complex plane for $r\Delta t$, a family of curves for different values of $d\Delta t$ and $a\Delta t$, and indicate the corresponding stability regions. Here we present the analysis of the IIF2 scheme (7) as an example. More details and analysis results can be found in [12].

Fig. 8

Linear stability regions of the IIF2 scheme (7) for different values of $d\Delta {t}$ under a fixed value of $a\Delta {t}$. a $a\Delta {t}=1.0$; b $a\Delta {t}=10.0$; c $a\Delta {t}=-1.0$; d $a\Delta {t}=-10.0$

Fig. 9

Linear stability regions of the IIF2 scheme (7) for different values of $a\Delta {t}$ under a fixed value of $d\Delta {t}$. a $d\Delta {t}=1.0$; b $d\Delta {t}=2.0$; c $d\Delta {t}=10.0$; d $d\Delta {t}=20.0$

Applying the IIF2 scheme (7) to the Eq. (23) with a uniform time step size $\Delta t$, then substituting $u_n=e^{in\theta }$ into the resulting equation, we obtain

$$\begin{aligned} \left( 1-\frac{\lambda }{2}\right) e^{2i\theta }=e^{-d\Delta t}\left( 1+\frac{\lambda }{2}+\frac{3}{2}a\Delta t\right) e^{i\theta } -\frac{a}{2}\Delta t e^{-2d \Delta t}, \end{aligned}$$

(24)

where $\lambda =r\Delta t$ has a real part $\lambda _r$ and imaginary part $\lambda _i$. Solve the Eq. (24) for $\lambda _r$ and $\lambda _i$ to have

$$\begin{aligned} {\left\{ \begin{array}{ll} \lambda _r=\dfrac{B_1C_2-B_2C_1}{A_1B_2-A_2B_1};\\ \lambda _i=\dfrac{A_1C_2-A_2C_1}{A_2B_1-A_1B_2}, \end{array}\right. } \end{aligned}$$

(25)

where

$$\begin{aligned} {\left\{ \begin{array}{ll}A_1=e^{-d \Delta t}\dfrac{1}{2}\cos \theta +\dfrac{1}{2}\cos 2 \theta ,\\ B_1=-e^{-d \Delta t}\dfrac{1}{2}\sin \theta -\dfrac{1}{2}\sin 2 \theta ,\\ C_1=-\dfrac{a}{2}\Delta t e^{-2d \Delta t}+e^{-d \Delta t}\left( 1+\dfrac{3}{2}a \Delta t\right) \cos \theta -\cos 2\theta ,\\ A_2=e^{-d \Delta t}\dfrac{1}{2}\sin \theta +\dfrac{1}{2}\sin 2 \theta ,\\ B_2=e^{-d \Delta t}\dfrac{1}{2}\cos \theta +\dfrac{1}{2}\cos 2 \theta ,\\ C_2=e^{-d \Delta t}\left( 1+\dfrac{3}{2}a \Delta t\right) \sin \theta -\sin 2\theta . \end{array}\right. } \end{aligned}$$

(26)

Stability regions in the complex plane of $r\Delta t$ for different values of $d\Delta t$ under a fixed value of $a \Delta t$ are presented in Fig. 8. As examples we choose four different $a \Delta t$ values: $a \Delta t=1.0$, $a \Delta t=10.0$ , $a \Delta t=-1.0$ and $a \Delta t=-10.0$. The points on boundaries of stability regions are obtained by varying $\theta$ from 0 to $2\pi$ in (25) and (26). A stability boundary curve divides the whole complex plane into the stable region and the unstable region for a pair of fixed values of $d \Delta t$ and $a \Delta t$. Based on analyzing the growth factor of the scheme (7) for some special values of $d\Delta t$, $a \Delta t$ and $\lambda$, we find that the stable regions always include the point $\lambda =(-20,0)$ for any values of $d \Delta t$ and $a \Delta t$ used in Fig. 8. Then stable and unstable regions are determined and shown in Fig. 8. From Fig. 8, we can see that the whole regions outside of the stability boundary curves are stable regions, which shows that the IIF2 scheme (7) has large stability regions. For a fixed $a\Delta t$, the stable region becomes larger with the increase of the value of $d \Delta t$. Next we show stability regions for different values of $a\Delta t$ under a fixed value of $d \Delta t$ in Fig. 9. $d \Delta t=1.0$ , $d \Delta t=2.0$, $d \Delta t=10.0$ and $d \Delta t=20.0$ are chosen as examples. Again, analysis of the growth factor of the scheme (7) for some special values of $d\Delta t$, $a \Delta t$ and $\lambda$, we find that the stable regions always include the point $\lambda =(-10,0)$ for any values of $a \Delta t$ and $d \Delta t$ used in Fig. 9. Stable regions for the cases shown in Fig. 9 are the whole regions outside of the stability boundary curves. For a fixed $d\Delta t$, the stable region becomes smaller with the increase of the value of $|a| \Delta t$ which corresponds to the convection terms. Based on the linear stability analysis, we conclude that the diffusion term tends to stabilize the scheme, while the convection term gives constraints on time step sizes. Due to the implicit property of the scheme, the stability regions are quite large and often include the whole left complex plane, with a relatively large size diffusion parameter d and a mild size convection parameter a.