Many-Stage Optimal Stabilized Runge-Kutta Methods for Hyperbolic Partial Differential Equations

TL;DR

This work introduces a nonlinear optimization framework to generate high-degree ($>100$) stability polynomials for stabilized explicit Runge-Kutta methods, aimed at semidiscretizations of hyperbolic PDEs with spectra that include imaginary components. By parameterizing the stability polynomial through complex-conjugate pseudo-extrema and leveraging a convex-hull interpolant of the spectrum, the authors overcome ill-conditioning issues inherent in monomial coefficient formulations and construct many-stage, low-storage schemes that maintain second-order accuracy on a range of linear and nonlinear problems. They extend the approach to non-convex spectra via alpha shapes or convex hulls, and demonstrate polynomial degrees up to 128 for several canonical PDEs, including Burgers, shallow water, MHD, and 2D Euler. Internal stability considerations drive a heuristic optimization for the stage coefficients, enabling realizable implementations despite the lack of SSP guarantees. The resulting methods show substantial stability advantages, enabling larger timesteps in practice, with demonstrated convergence on both linear and nonlinear test problems, though care is needed to manage oscillations and spectral effects in nonlinear regimes.

Abstract

A novel optimization procedure for the generation of stability polynomials of stabilized explicit Runge-Kutta methods is devised. Intended for semidiscretizations of hyperbolic partial differential equations, the herein developed approach allows the optimization of stability polynomials with more than hundred stages. A potential application of these high degree stability polynomials are problems with locally varying characteristic speeds as found in non-uniformly refined meshes and different wave speeds. To demonstrate the applicability of the stability polynomials we construct 2N storage many-stage Runge-Kutta methods that match their designed second order of accuracy when applied to a range of linear and nonlinear hyperbolic PDEs with smooth solutions. The methods are constructed to reduce the amplification of round off errors which becomes a significant concern for these many-stage methods.

Figures (11)

• Figure 1: Optimal second order accurate stability polynomials $P_{S, 2}(z)$ given by \ref{['eq:StabPolyCircleSecondOrderAccurate']} for even \ref{['fig:CircleEven']} and odd degree \ref{['fig:CircleOdd']} polynomial. The spectrum corresponds to the canonical first order Finite Volume Upwind/Godunov discretization godunov:hal-01620642 of the advection equation $u_t + u_x = 0$ on the periodic $[-1, 1]$ domain discretized with $500$ cells. In the even degree case \ref{['fig:CircleEven']}, the lower degree polynomial is odd and thus there is one pure real pseudo extremum at the left end of the spectrum. For odd polynomial degrees \ref{['fig:CircleOdd']}, the lower degree polynomial is of even degree and we have only complex-conjugated pseudo-extrema. Note also that the segment crossing $0$ is twice the length of the others, which follows from the fact that $P_{S, p}(0) \equiv 1 \: \forall \: S, p \geq 1$, cf. \ref{['eq:DefinitionLowerDegree']}.
• Figure 2: Collection of strictly convex spectra of nonlinear hyperbolic pdes and optimized pseudo-extrema $\widetilde{r}_j$. The spectra are scaled with the optimal timestep $\Delta t_{16, 2}$ in each case.
• Figure 3: Scaled spectrum $\sigma_{12, 2}$ of the 1D linear advection equation $u_t +u_x =0$ discretized through the dgsem on $[0, 1]$ using $16$ cells/elements with DG polynomial degree $3$ and Rusanov/Local Lax-Friedrichs flux. The pseudo-extrema of the $12$ and $24$ degree second order accurate polynomial are computed. To highlight that the relative positions of every second pseudo extremum $\tilde{r}_{2j-1}^{2S}$ agree with $\tilde{r}_j^S$, the former are scaled by $0.5$ due to the linear scaling of the timestep according to \ref{['eq:TimestepLinearScaling']}.
• Figure 4: Collection of strictly convex spectra of nonlinear hyperbolic pdes and optimized pseudo-extrema. The spectra are scaled with the optimal timestep $\Delta t_{16, 2}$ in each case.
• Figure 5: Approximation of the optimal stability boundary by the contour of an alpha shape. The spectrum $\sigma$ is obtained from a dgsem discretization of the 2D linear advection equation with velocities $a_x = 0.5, a_y = -0.1$ on the periodic domain $\Omega = [-1,1]^2$ discretized by a $6\times 6$ mesh with local polynomials of degree $3$ and Rusanov/Local Lax-Friedrichs flux.

Theorems & Definitions (2)

• Theorem 1
• Theorem 2