Numerical approximation of Caputo-type advection-diffusion equations via Sylvester equations

TL;DR

This work develops a high-order numerical scheme for Caputo-type advection-diffusion equations with time nonlocal derivatives $D_t^{\alpha}$ on $\mathbb{R}$ or bounded intervals. It combines a $3-\alpha$-order time discretization with spectral space discretization and recasts the PDE as a Sylvester equation $\mathbf A\mathbf U + \mathbf U\mathbf B = \mathbf C$ to compute the solution at all time levels simultaneously; a fast FFT-based convolution and an explicit operational matrix $\mathbf D_t^{\alpha}$ accelerate the time discretization, while careful boundary handling enables accurate bounded-domain simulations. The approach is demonstrated on unbounded and bounded domains, including a literature-inspired test, achieving high accuracy in MATLAB and enabling extension to more general domains and higher-order spatial operators. Overall, the framework provides an efficient and reproducible tool for solving time-fractional PDEs in engineering and physics, with potential applicability to semi-infinite domains and broader discretizations.

Abstract

In this paper, we approximate numerically the solution of Caputo-type advection-diffusion equations of the form $D_t^α u(t,x) = a_1(x)u_{xx}(t,x) + a_2(x)u_x(t,x) + a_3u(t,x) + a_4(t,x)$, where $D_t^α u$ denotes the Caputo fractional derivative of order $α\in(0,1)$ of $u$ with respect to $t$, $t\in[0, t_f]$ and the spatial domain can be the whole real line or a closed interval. First, we propose a method of order $3 - α$ to approximate Caputo fractional derivatives, explain how to implement an FFT-based fast convolution to reduce the computational cost, and express the numerical approximation in terms of an operational matrix. Then, we transform a given Caputo-type advection-diffusion equation into a Sylvester equation of the form $\mathbf A\mathbf U + \mathbf U \mathbf B = \mathbf C$, and special care is given to the treatment of the boundary conditions, when the spatial domain is a closed interval. Finally, we perform several numerical experiments that illustrate the adequacy of our approach. The implementation has been done in Matlab, and we share and explain in detail the majority of the actual codes that we have used.

Figures (2)

• Figure 1: Errors $|D_{t,num}^{\alpha} e^{2t_j} - D_t^{\alpha} e^{2t_j}|$ (in black) and $|D_{t,num}^{\alpha,*} e^{2t_j} - D_t^{\alpha} e^{2t_j}|$ (in red) corresponding respectively to the numerical approximations $D_{t,num}^{\alpha} e^{2t_j}$ (given by \ref{['e:cptftj']}) and $D_{t,num}^{\alpha,*} e^{2t_j}$ (given by \ref{['e:cptftjstar']}) of $D_t^{\alpha} e^{2t}$, taking $t_f = 1.2$, $\alpha = 0.17$, and $N\in\{100, 200, 400, 800, 1600\}$.
• Figure 2: Left: $\log_{10}(\max_{1\le j\le N} |D_{t,num}^{\alpha} e^{2t_j} - D_t^{\alpha} e^{2t_j}|)$ versus $\log_2(N)$, for $\alpha\in\{0.05, 0.1, \ldots, 0.95\}$ and $N\in\{2^1, 2^2, \ldots, 2^{20}\} = \{2, 4, \ldots, 1048576\}$. Right: numerical approximation of the order of convergence (in red) for $\alpha\in\{0.05, 0.1, \ldots, 0.95\}$ compared to the theoretical values $3 - \alpha$ (dash-dotted black line).

Theorems & Definitions (2)

• Remark
• Remark