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

In this paper, we approximate numerically the solution of Caputo-type advection-diffusion equations of the form $D_t^{\alpha} 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^{\alpha} u$ denotes the Caputo fractional derivative of order $\alpha\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 - \alpha$ 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.