  The purpose of this paper is twofold. Firstly, we provide explicit and
compact formulas for computing both Caputo and (modified) Riemann-Liouville
(RL) fractional pseudospectral differentiation matrices (F-PSDMs) of any order
at general Jacobi-Gauss-Lobatto (JGL) points. We show that in the Caputo case,
it suffices to compute F-PSDM of order $\mu\in (0,1)$ to compute that of any
order $k+\mu$ with integer $k\ge 0,$ while in the modified RL case, it is only
necessary to evaluate a fractional integral matrix of order $\mu\in (0,1).$
Secondly, we introduce suitable fractional JGL Birkhoff interpolation problems
leading to new interpolation polynomial basis functions with remarkable
properties: (i) the matrix generated from the new basis yields the exact
inverse of F-PSDM at "interior" JGL points; (ii) the matrix of the highest
fractional derivative in a collocation scheme under the new basis is diagonal;
and (iii) the resulted linear system is well-conditioned in the Caputo case,
while in the modified RL case, the eigenvalues of the coefficient matrix are
highly concentrated. In both cases, the linear systems of the collocation
schemes using the new basis can solved by an iterative solver within a few
iterations. Notably, the inverse can be computed in a very stable manner, so
this offers optimal preconditioners for usual fractional collocation methods
for fractional differential equations (FDEs). It is also noteworthy that the
choice of certain special JGL points with parameters related to the order of
the equations can ease the implementation. We highlight that the use of the
Bateman's fractional integral formulas and fast transforms between Jacobi
polynomials with different parameters, are essential for our algorithm
development.
