  We propose a novel numerical algorithm utilizing model reduction for
computing solutions to stationary partial differential equations involving the
spectral fractional Laplacian. Our approach utilizes a known characterization
of the solution in terms of an integral of solutions to classical elliptic
problems. We reformulate this integral into an expression whose continuous and
discrete formulations are stable; the discrete formulations are stable
independent of all discretization parameters. We subsequently apply the reduced
basis method to accomplish model order reduction for the integrand. Our choice
of quadrature in discretization of the integral is a global Gaussian quadrature
rule that we observe is more efficient than previously proposed quadrature
rules. Finally, the model reduction approach enables one to compute solutions
to multi-query fractional Laplace problems with order of magnitude less cost
than a traditional solver.
