Adaptive spectral solver for Riesz fractional reaction–diffusion equations via penalized minimum residual iteration

A high-order solver is presented for two-dimensional Riesz fractional nonlinear reaction–diffusion equations. It employs a midpoint starter and a three-point backward differentiation formula (BDF2) to achieve second-order temporal accuracy, together with a weighted Jacobi spectral approximation that delivers nearly exponential spatial convergence for analytic solutions. After Newton linearization, each correction is obtained via a penalized Levenberg–Marquardt minimum residual method (PLM-MRM). This iteration adaptively enforces boundary conditions without requiring boundary-fitted basis functions. We establish stability and rigorous a priori error bounds. Numerical experiments over a wide range of fractional orders confirm these rates and drive the residual to machine precision within a few PLM-MRM sweeps. Compared with a conventional LM update, global errors are reduced by up to 35%, and by one to two orders of magnitude relative to Galerkin-BDF or Crank–Nicolson (CN) baselines. For a given accuracy, the scheme allows time steps up to about four times larger than a recent fourth-order CN method.