<figure><img src="image_1.png" /></figure>

HAL
open science

# Projective integration schemes for nonlinear degenerate parabolic systems

Tommaso Tenna

<figure><img src="image_2.png" /></figure>

## To cite this version:

Tommaso Tenna. Projective integration schemes for nonlinear degenerate parabolic systems. 2025. hal-04981107

HAL Id: hal-04981107

https://hal.science/hal-04981107v1

Preprint submitted on 6 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Open licence - etalab
PROJECTIVE INTEGRATION SCHEMES FOR NONLINEAR
DEGENERATE PARABOLIC SYSTEMS

TOMMASO TENNA

ABSTRACT. A general high-order fully explicit scheme based on projective integration methods is here presented to solve systems of degenerate parabolic equations in general dimensions. The method is based on a BGK approximation of the advection-diffusion equation, where we introduce projective integration method as time integrator to deal with the stiff relaxation term. This approach exploits the clear gap in the eigenvalues spectrum of the kinetic equation, taking into account a sequence of small time steps to damp out the stiff components, followed by an extrapolation step over a large time interval. The time step restriction on the projective step is similar to the CFL condition for advection-diffusion equations. In this paper we discuss the stability and the consistency of the method, presenting some numerical simulations in one and two spatial dimensions.

2020 MATHEMATICS SUBJECT CLASSIFICATION: 65M08, 65L04, 65M12,

# 1. INTRODUCTION

The numerical study of degenerate parabolic equations has gained particular interest in the last decades. These problems arise in several physical applications, such as flow of gas in porous media [13], radiative transport [42] or cell chemotactic movement [40]. The numerical approximation of these equations is challenging due to the presence of possible degenerate points, in which the parabolic system loses its properties, exhibiting hyperbolic behavior with finite speed of propagation and sharp fronts. In general, classical solutions might not exist even for $C^2$ initial data and weak solutions must be considered. Many different schemes have been proposed to solve degenerate parabolic equations, such as local discontinuous Galerkin methods [47], mixed finite element methods [1], central finite volume schemes [24], finite volume schemes preserving steady-states [8], high order finite difference WENO schemes [32] and WENO schemes based on deep learning techniques [23]. Among these methods, a discrete kinetic approximation of these equations have been proposed by Aregba-Driollet, Natalini and Tang [5], inspired by relaxation schemes [20, 4] and kinetic approximation of hydrodynamic equations in the hyperbolic setting [30]. The idea relies on a Bhatnagar-Gross-Krook (BGK) approximation of the macroscopic equation, where both the kinetic velocities and the source term depend singularly on the relaxation parameter. Different relaxation systems for the approximation

<footnote>(Tommaso Tenna) UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE</footnote>

<footnote>E-mail address: tommaso.tenna@univ-cotedazur.fr.</footnote>

<footnote>Key words and phrases. Advection-diffusion equations; Relaxation schemes; High-order methods; Projective integration;</footnote>
of degenerate parabolic equations have been proposed in [35, 37], with high-order extension developed in [11], where the authors coupled ENO and WENO schemes for space discretization with IMEX schemes for time advancement. This approach is inspired by kinetic schemes for the Carleman model, but the numerical integration is directly performed at the macroscopic level. In the context of relaxation schemes, other diffusive kinetic models and approximations have been proposed, see for instance [31, 19, 36]. The main advantage of relaxation schemes is the possibility to treat the scalar and the system case in the same way at the numerical level, with the feasibility of parallelizing the algorithm due to the diagonal form of the approximating problems. Indeed discrete kinetic approximations are based on linear hyperbolic terms, whose structure is very likely for numerical and theoretical purpose [4, 5, 38]. The nonlinearity inside the derivatives is replaced by a semilinearity: the differential part becomes linear and all the nonlinearity is concentrated inside the source term. On the other hand, the main disadvantage of this approach is the appearance of a stiff source term, describing the relaxation of the new variable towards the equilibrium, which requires suitable time integrators.

In this perspective, the idea is to use an efficient explicit time integrator which is able to speed up the computation due to the particular eigenvalues separation in the spectrum of the kinetic formulation. Projective integration method was introduced in [15] as a technique to accelerate a brute-force integration of stiff systems of ordinary differential equations, which require restrictive conditions on the time step. Such stiff problems are characterized by the presence of different scales, that correspond to a clear gap in the eigenvalue spectrum, in which the fast modes (corresponding to the eigenvalues with large negative part) decay quickly, whereas the slow modes are the active components, still present in the solution. The idea, extended in [27] for kinetic equations with a diffusive scaling, is based on an extrapolation in time over a large time step, performed after taking a few inner steps to damp out the transient corresponding to the fast modes. Several extensions have been proposed both for hyperbolic problems [26] and for high order approximation of kinetic equations [25]. We mention also telescopic projective integration techniques, where the projective integration method have been extended to deal with problems having multiple time scales, both for ODEs systems [16] and for kinetic equations [34, 33, 6].

In [5, 11] the methods are based on splitting techniques, restricting the order in time to 2. The order of convergence in time could be increased by solving explicitly the BGK equation, without projecting over the Maxwellian space. Many asymptotic-preserving schemes based on IMEX techniques have been proposed to integrate stiff kinetic equations and we refer to the recent survey [9] for a detailed analysis on the topic.

The idea of this paper is to propose a fully explicit and robust projective integration scheme to numerically solve the stiff BGK system with arbitrary order of accuracy in time, extending the idea of [26] to the case of possibly degenerate parabolic systems.

The plan of the paper is the following: in Section 2 we introduce discrete kinetic relaxation schemes for convection-diffusion equations, investigating different models. We construct the numerical scheme, focusing on the phase space discretization. In Section 3 we introduce projective integration method to treat the time discretization with a fully explicit approach, by exploiting the particular structure of the BGK operator. Based
on [25], we investigate the spectrum of the operator and we perform stability and consistency analysis of the method. In Section 4 we present some numerical simulations in both one and two space dimensions for different problems. All the details about the numerical schemes are reported in the Appendix.

## 2. RELAXATION SCHEMES FOR NONLINEAR DEGENERATE PARABOLIC SYSTEMS

In this section, we study discrete kinetic schemes for systems of conservation laws with possibly degenerate diffusion, referring to [5], but other kinetic models and approximations can be built, see for instance [19]. Let $u: \mathbb{R}^D \times \mathbb{R}^+ \to \mathbb{R}^K$, $D > 0$, $K > 0$, be a weak solution of the following system

$$ (1) \qquad \partial_t u + \sum_{d=1}^{D} \partial_{x_d} A_d(u) = \Delta_x B(u), \quad (x,t) \in \mathbb{R}^D \times (0, +\infty) $$

for a given initial condition

$$ (2) \qquad u(x, 0) = u_0(x), \quad x \in \mathbb{R}^D. $$

Here $t$ defines the time, $x_d$ defines the spatial variable for each dimension $d = 1, \dots, D$ and $u_0: \mathbb{R}^D \to \mathbb{R}^K$. The functions $A_d: \mathbb{R}^K \to \mathbb{R}^K$ and $B: \mathbb{R}^K \to \mathbb{R}^K$ are assumed to be Lipschitz-continuous. In addition, for all $u$ lying in a fixed domain $\Omega \in \mathbb{R}^D$ it holds true:

* For all $\xi \in \mathbb{R}^D$, $\sum_{d=1}^{D} \xi_d A'_d(u)$ has real eigenvalues and it is diagonalizable,

* The real part of the eigenvalues of $B'(u)$ is non-negative.

In [4] a family of discrete diffusive kinetic schemes has been proposed, in order to construct a numerical approximation of the original system (1). Let us consider $f^\epsilon: \mathbb{R}^D \times \mathbb{R}^+ \to \mathbb{R}^L$ as the solution of the following kinetic model in general form:

$$ (3) \qquad \left\{ \begin{array}{l} \partial_t f^\epsilon + \sum_{d=1}^{D} \Gamma_d^\epsilon \partial_{x_d} f^\epsilon = \frac{1}{\epsilon}(M(\epsilon, u^\epsilon) - f^\epsilon), \\ u^\epsilon(x,t) = \sum_{l=1}^{L} f_l^\epsilon(x,t), \end{array} \right. $$

where $L > 0$ is the number of discrete velocities and let us take the initial condition

$$ (4) \qquad f^\epsilon(x, 0) = f_0^\epsilon(x). $$

Here, $M: \mathbb{R}^+ \times \mathbb{R}^K \to \mathbb{R}^L$ are the Maxwellians and $\Gamma_d^\epsilon \in \mathbb{R}^{L \times L}$ are constant diagonal matrices, where the diagonal entries $\gamma_{ld}^\epsilon$ are given by

$$ (5) \qquad \gamma_{ld}^\epsilon := \lambda_{ld} + \frac{\theta_{ld}}{\sqrt{\epsilon}}, $$

for some fixed real constants $\lambda_{ld}, \theta_{ld}$. The relaxation parameter $\epsilon \in \mathbb{R}^+$ can be a physical parameter or an artificial one and the kinetic system (3) converges (at least formally) to the original macroscopic system as $\epsilon \to 0$, under an analogous of the *subcharacteristic condition*.

The Maxwellian $M$ essentially embeds the macroscopic variable $u^\epsilon$ in the kinetic space
and the right hand side models a relaxation of the kinetic variable $f^\epsilon$ towards the equilibrium state given by the Maxwellian $\mathcal{M}$. From now on, let us consider $\Omega \in \mathbb{R}^K$ a bounded domain, such that $u(x,t) \in \Omega$ for all $(x,t) \in \mathbb{R}^D \times \mathbb{R}^+$. When $\epsilon \to 0$, the kinetic approximation formally converges to the macroscopic model under some fundamental assumptions on the Maxwellian functions:

$$ (M1) \quad \sum_{l=1}^{L} \mathcal{M}_l(\epsilon, w) = w \quad \text{for all } \epsilon \in (0,1] \text{ and for all } w \in \Omega, \text{ for all } d=1, \dots, D $$

$$ (M2) \quad \sum_{l=1}^{L} \gamma_{ld}^{\epsilon} \mathcal{M}_l(\epsilon, w) = A_d(w) \quad \text{for all } \epsilon \in (0,1] \text{ and for all } w \in \Omega, \text{ for all } d=1, \dots, D $$

$$ (M3) \quad \sum_{l=1}^{L} \theta_{ld} \theta_{lj} \mathcal{M}_l(0, w) = \delta_{dj} B(w) \quad \text{for all } w \in \Omega, \text{ for all } d, j = 1, \dots, D $$

$$ (M4) \quad \mathcal{M}_l(\epsilon, w) \to \mathcal{M}_l(0, w) \quad \text{when } \epsilon \to 0 \text{ uniformly for } w \in \Omega. $$

If $\mathcal{M}$ satisfies the above conditions it is called a local Maxwellian Function for system (1) on $\Omega$. The first property tells us that the projection of the Maxwellian is the identity, which means that considering a given macroscopic variable embedded in the microscopic space and projecting it back, we obtain the original state. The second property and the third property are necessary to guarantee the consistency with the macroscopic fluxes and, therefore, to preserve the structure of the macroscopic model in the limit of the kinetic model (3) towards (1). The last property is fundamental to obtain a uniform convergence with respect to $\epsilon$, when $\epsilon \to 0$.

The stability condition for this model deals with the following definition.

**Definition 2.1.** A local Maxwellian function is a monotone Maxwellian function (MMF) if for all $l \in \{1, \dots, L\}$, for all positive $\epsilon \le 1$ and for all $u \in \Omega$, the real parts of the eigenvalues of $M'_l(u)$ are nonnegative.

Rigorous theoretical results about the convergence of the solution to the BGK model (3) towards a weak solution to the system (1) have been proven for the scalar case in [10] and for a class of one dimensional strongly degenerate parabolic systems in [28]. In the next subsection we give a formal derivation of the macroscopic equation, starting from the Chapman-Enskog expansion and assuming the properties (M1)-(M4).

**2.1. Chapman-Enskog expansion.** Let us define the projection matrix $P \in \mathbb{R}^{L \times K}$ and the auxiliary variables

$$ (6) \quad u^{\epsilon} := P f^{\epsilon} = \sum_{l=1}^{L} f_{l}^{\epsilon}, \quad v_{d}^{\epsilon} = P \Gamma_{d} f^{\epsilon} = \sum_{l=1}^{L} \gamma_{ld}^{\epsilon} f_{l}^{\epsilon} \quad d = 1, \dots, D. $$
Then from (3) and the compatibility assumptions (M1)-(M4), we have

$$ (7) \qquad \left\{ \begin{array}{l} \partial_t u^\varepsilon + \sum_{d=1}^D \partial_{x_d} v_d^\varepsilon = 0, \\ \partial_t v_d^\varepsilon + \sum_{j=1}^D \partial_{x_j} (P \Gamma_d \Gamma_j f^\varepsilon) = \frac{1}{\varepsilon} (A_d(u^\varepsilon) - v_d^\varepsilon), \\ \sum_{l=1}^L \theta_{ld} \theta_{lj} f_l^\varepsilon = \delta_{dj} B(u^\varepsilon) \end{array} \right. $$

Let us consider a formal Chapman-Enskog expansion of $f^\varepsilon$ in the form

$$ (8) \qquad f^\varepsilon = \mathcal{M}(u^\varepsilon) + \varepsilon g^\varepsilon + \mathcal{O}(\varepsilon^2), $$

Then, from (7), one has

$$ (9) \qquad v_d^\varepsilon = A_d(u^\varepsilon) - \varepsilon \left( \partial_t v_d^\varepsilon + \sum_{j=1}^D \partial_{x_j} (P \Gamma_d \Gamma_j f^\varepsilon) \right) + \mathcal{O}(\varepsilon^2), $$

from which we have

$$ (10) \qquad v_d^\varepsilon = A_d(u^\varepsilon) - \varepsilon \left( \partial_t v_d^\varepsilon + \sum_{j=1}^D \partial_{x_j} \sum_{l=1}^L \gamma_{dj} \gamma_{lj} f_l^\varepsilon \right) + \mathcal{O}(\varepsilon^2). $$

This implies

$$ (11) \qquad \partial_t u^\varepsilon + \sum_{d=1}^D \partial_{x_d} A_d(u^\varepsilon) = \varepsilon \sum_{d=1}^D \partial_{x_d} \left( \partial_t v_d^\varepsilon + \sum_{j=1}^D \partial_{x_j} \sum_{l=1}^L \gamma_{dj} \gamma_{dl} f_l^\varepsilon \right) + \mathcal{O}(\varepsilon^2). $$

Then, we have

$$ (12) \qquad \partial_t v_d^\varepsilon = A'_d(u^\varepsilon) \partial_t u^\varepsilon + \mathcal{O}(\varepsilon) = - \sum_{j=1}^D A'_d(u^\varepsilon) A'_j(u^\varepsilon) \partial_{x_j} u^\varepsilon + \mathcal{O}(\varepsilon^2), $$

which leads

$$ \begin{aligned} & \partial_t u^\varepsilon + \sum_{d=1}^D \partial_{x_d} A_d(u^\varepsilon) \\ & = \sum_{d=1}^D \sum_{j=1}^D \partial_{x_d} \partial_{x_j} \sum_{l=1}^L \theta_{dj} \theta_{dl} f_l^\varepsilon + 2\sqrt{\varepsilon} \sum_{d=1}^D \sum_{j=1}^D \partial_{x_d} \sum_{l=1}^L \lambda_{dj} \theta_{dl} f_l^\varepsilon \\ & + \varepsilon \left( \sum_{d=1}^D \sum_{j=1}^D \partial_{x_d} \sum_{l=1}^L \lambda_{dj} \lambda_{dl} f_l^\varepsilon + \sum_{d=1}^D \partial_{x_d} \left( \sum_{j=1}^D A'_d(u^\varepsilon) A'_j(u^\varepsilon) \partial_{x_j} u^\varepsilon \right) \right) + \mathcal{O}(\varepsilon^2). \end{aligned} $$

Using (M3), when $\varepsilon \to 0$, we finally obtain

$$ (13) \qquad \partial_t u^\varepsilon + \sum_{d=1}^D \partial_{x_d} A_d(u^\varepsilon) = \Delta_x B(u^\varepsilon), $$

requiring that the leading order terms in the $\varepsilon$-expansion are positive, in the same spirit of the so called *subcharacteristic condition* for the hyperbolic case. This formal justification
shows that the natural convergence rate in (M4) would be $\sqrt{\epsilon}$. The rigorous proof of this result has been obtained in [10], using a priori bounds on the discrete velocity BGK approximation and kinetic entropy inequalities.

2.2. Discrete kinetic models. In the pure hyperbolic setting, system (3) reduces to the one developed in [38] and [4]. Let us consider the following hyperbolic system of conservation laws

$$ (14) \qquad \partial_t u + \sum_{d=1}^{D} \partial_{x_d} A_d(u) = 0. $$

and the following approximating discrete BGK model

$$ (15) \qquad \left\{ \begin{array}{l} \partial_t f_l^\epsilon + \sum_{d=1}^D \lambda_{ld} \partial_{x_d} f_l^\epsilon = \frac{1}{\epsilon} (\tilde{M}_l(u^\epsilon) - f_l^\epsilon), \quad l = 1, \dots, J, \\ u^\epsilon = P f^\epsilon = \sum_{l=1}^J f_l^\epsilon, \end{array} \right. $$

where $J$ is the number of velocities needed to approximate the hyperbolic equation, complemented with the compatibility conditions

$$ (16) \qquad \sum_{l=1}^{J} \tilde{M}_l(u) = u, \qquad \sum_{l=1}^{J} \lambda_{ld} \tilde{M}_l(u) = A_d(u), \quad \forall u \in \Omega, \quad d = 1, \dots, D, $$

In [38] it has been shown that the solution of the discrete BGK model (15), under the hypothesis of nondecreasing Maxwellian functions, converges to the unique entropy solution of the hyperbolic system (14), see also [39].

In the same spirit, a model for (1) can be obtained by adding linear combinations of $B$ to the $\tilde{M}_l$, together with supplementary equations to take into account the diffusive scaling. More precisely, fixing $J' \ge D+1$ the number of velocities you need to deal with the parabolic term and $L = J + J'$, we can design the model as

$$ (17) \qquad \left\{ \begin{array}{l} \partial_t f_l^\epsilon + \sum_{d=1}^D \lambda_{ld} \partial_{x_d} f_l^\epsilon = \frac{1}{\epsilon} (\tilde{M}_l(u^\epsilon) - f_l^\epsilon), \quad 1 \le l \le J, \\ \partial_t f_{J+m}^\epsilon + \left( \mu + \frac{\theta \sqrt{J'}}{\sqrt{\epsilon}} \right) \sum_{d=1}^D \sigma_{md}^\epsilon \partial_{x_d} f_{J+m}^\epsilon = \frac{1}{\epsilon} \left( \frac{B(u^\epsilon)}{J'\theta^2} - f_{J+m}^\epsilon \right), \quad 1 \le m \le J', \\ u^\epsilon(x,t) = \sum_{l=1}^L f_l^\epsilon(x,t), \end{array} \right. $$

In particular, for $1 \le l \le J$:

$$ (18) \qquad \mathcal{M}_l(u) = \tilde{M}_l(u) + b_l B(u), $$

where $b_l \in \mathbb{R}$. On the other hand, for the last $J'$ equations, we have

$$ (19) \qquad \theta \ge 0, \mu \ge 0 $$

and $\{\sigma_d\}_{d=1}^D$ is an orthonormal family in $\{X \in \mathbb{R}^{J'} \text{ s.t. } \sum_{m=1}^{J'} X_m = 0\} = H$.
2.2.1. Diagonal Relaxation Model 1 (DRM1). Let us first focus on the one-dimensional case, where we choose $J = 2$ and $J' = 2$. Then the kinetic model reads as

$$ (20) \qquad \partial_t f^\varepsilon + \Gamma^\varepsilon \partial_x f^\varepsilon = \frac{1}{\varepsilon}(\mathcal{M}(u^\varepsilon) - f^\varepsilon), $$

where

$$ (21) \qquad \left\{ \begin{array}{l} \Gamma^\varepsilon = \operatorname{diag} \left( -\lambda, \lambda, -\frac{\mu}{\sqrt{2}} - \frac{\theta}{\sqrt{\varepsilon}}, \frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\varepsilon}} \right)^T, \\ \mathcal{M}(u) = \frac{1}{2} \left( u - \frac{A(u)}{\lambda} - \frac{B(u)}{\theta^2}, u + \frac{A(u)}{\lambda} - \frac{B(u)}{\theta^2}, \frac{B(u)}{\theta^2}, \frac{B(u)}{\theta^2} \right)^T \end{array} \right. $$

Note that the Maxwellian functions do not depend anymore on $\varepsilon$.

Let us suppose that for all $u$, $A'(u)$ and $B'(u)$ have a basis of common eigenvectors, denoting $\lambda_k(U)$ and $\theta_k^2(u)$ their respective eigenvalues for $1 \le k \le K$ and let us assume that the eigenvalues of $B'$ are real, then $\mathcal{M}$ is a Monotone Maxwellian Function if

$$ (22) \qquad \sup_{u \in \Omega} \sup_{1 \le k \le K} \frac{|\lambda_k(u)|}{\lambda} + \frac{\theta_k^2(u)}{\theta^2} \le 1, $$

which is a generalized subcharacteristic condition [5].

Let us consider now the two-dimensional model, in which $J' \ge D + 1 = 3$, $\theta > 0$, $\mu \ge 0$. Let $\{\sigma^{(1)}, \sigma^{(2)}\}$ be orthonormal in $H$.

$$ (23) \qquad \gamma_1^\varepsilon = \left( \begin{array}{c} \lambda_1(-1, 0, 1)^T \\ \left( \mu + \frac{\theta \sqrt{J'}}{\sqrt{\varepsilon}} \right) \sigma^{(1)} \end{array} \right), \qquad \gamma_2^\varepsilon = \left( \begin{array}{c} \lambda_2(-1, 0, 1)^T \\ \left( \mu + \frac{\theta \sqrt{J'}}{\sqrt{\varepsilon}} \right) \sigma^{(2)} \end{array} \right). $$

The Maxwellian functions are:

$$ (24) \qquad \mathcal{M}(u) = \frac{1}{3} \left( \begin{array}{c} u - 2 \frac{A_1(u)}{\lambda_1} + \frac{A_2(u)}{\lambda_2} \\ u - 2 \frac{A_1(u)}{\lambda_1} + \frac{A_2(u)}{\lambda_2} \\ u - 2 \frac{A_1(u)}{\lambda_1} + \frac{A_2(u)}{\lambda_2} \\ (0, \dots, 0)^T \end{array} \right) + \frac{B(u)}{\theta^2} \left( \begin{array}{c} -\frac{1}{3} (1, 1, 1)^T \\ \frac{1}{J'} (1, \dots, 1)^T \end{array} \right). $$

The Maxwellian $\mathcal{M}$ is a MMF in the scalar case ($K = 1$) if the following condition is satisfied:

$$ (25) \qquad \max_{u \in I} \left( \frac{2 A_1'}{\lambda_1} - A_2' \lambda_2, -\frac{A_1'}{\lambda_1} + 2 A_2' \lambda_2, -\frac{A_1'}{\lambda_1} - A_2' \lambda_2 \right) \le 1 - \frac{B'}{\theta^2}, $$

which is equivalent to the subcharacteristic condition presented in [4], for the hyperbolic case $\theta = 0$.
2.2.2. *Diagonal Relaxation Model 2 (DRM2)*. In this model, we consider the same structure of (DRM1) replacing $\lambda$ by two distinct values $\lambda_m$ and $\lambda_p$. In the one-dimensional case the kinetic model is given by (20), where now

$$ (26) \qquad \left\{ \begin{array}{l} \Gamma^{\varepsilon} = \operatorname{diag} \left( \lambda_m, \lambda_p, -\frac{\mu}{\sqrt{2}} - \frac{\theta}{\sqrt{\varepsilon}}, \frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\varepsilon}} \right)^{\mathrm{T}}, \\ \mathcal{M}_1(u) = \frac{1}{\lambda_p - \lambda_m} \left( \lambda_p \left( u - \frac{B(u)}{\theta^2} \right) - A(u) \right), \\ \mathcal{M}_2(u) = \frac{1}{\lambda_p - \lambda_m} \left( -\lambda_m \left( u - \frac{B(u)}{\theta^2} \right) + A(u) \right), \\ \mathcal{M}_3(u) = \mathcal{M}_4(u) = \frac{B(u)}{2\theta^2}. \end{array} \right. $$

Supposing w.l.o.g. that $\lambda_m < \lambda_p$, the Maxwellian $\mathcal{M}$ is a MMF if

$$ (27) \qquad \lambda_m \left(1 - \frac{B'(u)}{\theta^2}\right) \le A'(u) \le \lambda_p \left(1 - \frac{B'(u)}{\theta^2}\right), \quad \forall u \in I. $$

This condition allows for a sharper approximation of $A'$ than the previous one, suggesting that the numerical approximation could provide better results [5]. More specifically, we take

$$ (28) \qquad \lambda_m = \inf_{u \in I} \frac{A'(u)}{1 - B'(u)/\theta^2}, \qquad \lambda_p = \sup_{u \in I} \frac{A'(u)}{1 - B'(u)/\theta^2}, $$

to guarantee the maximality of $\lambda_m$ and the minimality of $\lambda_p$.

2.2.3. *A 3 velocities model (OVM)*. For the sake of completeness, we show another choice for the kinetic model, characterized by orthogonal velocities and $\varepsilon$-dependent Maxwellians

$$ (29) \qquad \partial_t f^{\varepsilon} + \Gamma^{\varepsilon} \partial_x f^{\varepsilon} = \frac{1}{\varepsilon} (\mathcal{M}(u^{\varepsilon}) - f^{\varepsilon}), $$

where

$$ (30) \qquad \left\{ \begin{array}{l} \Gamma^{\varepsilon} = \operatorname{diag}(0, \lambda, -\lambda), \\ \mathcal{M}_1(\varepsilon, u) = u - \frac{B(u)}{\theta^2}, \\ \mathcal{M}_2(\varepsilon, u) = \frac{\left(\lambda + \frac{\theta}{\sqrt{\varepsilon}}\right)}{2\left|\lambda + \frac{\theta}{\sqrt{\varepsilon}}\right|^2} A(u) + \frac{B(u)}{\theta^2}, \\ \mathcal{M}_3(\varepsilon, u) = -\frac{\left(\lambda + \frac{\theta}{\sqrt{\varepsilon}}\right)}{2\left|\lambda + \frac{\theta}{\sqrt{\varepsilon}}\right|^2} A(u) + \frac{B(u)}{\theta^2}. \end{array} \right. $$

In this case, if the equation (14) does not degenerate or for instance

$$ (31) \qquad B'(u) \ge |A'(u)| $$

holds, then the monotonicity conditions are verified for $\varepsilon$ small if $\theta^2 > \max B'$.
**Remark 2.1.** It is important to observe that by approximating the original system (1) with a kinetic model of the form (3), a modeling error is introduced and it is strictly $\epsilon$-dependent. However, in any numerical approximation of this model, a finite value of $\epsilon$ must be taken. This choice is dictated by the consistency of the numerical schemes with respect to the original problem.

2.3. **Numerical Approximation.** The main idea is to avoid splitting techniques to numerically solve (17). First, we introduce a phase space discretization using a discrete velocity model among the ones introduced above, for a general matrix

$$ (32) \qquad \Gamma^{\epsilon} = \operatorname{diag} \left( \lambda_1, \dots, \lambda_J, -\frac{\mu}{\sqrt{2}} - \frac{\theta}{\sqrt{\epsilon}}, \frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\epsilon}} \right), $$

coupled with a finite volume scheme for the advection term. Then, we propose a projective integration scheme for the time discretization.

2.3.1. *Velocity discretization.* Let us rewrite explicitly the semi-discrete problem in the one dimensional version:

$$ (33) \qquad \left\{ \begin{array}{ll} \partial_t f_l^\epsilon + \lambda_l \partial_x f_l^\epsilon = \frac{1}{\epsilon} (M_l(u^\epsilon) - f_l^\epsilon), & 1 \le l \le J, \\ \\ \partial_t f_{J+1}^\epsilon - \left(\frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\epsilon}}\right) \partial_x f_{J+1}^\epsilon = \frac{1}{\epsilon} \left(\frac{B(u^\epsilon)}{2\theta^2} - f_{J+1}^\epsilon\right), & \\ \\ \partial_t f_{J+2}^\epsilon + \left(\frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\epsilon}}\right) \partial_x f_{J+2}^\epsilon = \frac{1}{\epsilon} \left(\frac{B(u^\epsilon)}{2\theta^2} - f_{J+2}^\epsilon\right), & \\ \\ u^\epsilon(x,t) = \sum_{l=1}^{J+2} f_l^\epsilon(x,t). \end{array} \right. $$

The Maxwellian function $M_l$ for the $l$-th equation of the system is chosen as one of the Maxwellian described in Section 2.2.

In a BGK approximation (3), the minimal number $L$ of discrete velocities depends on the spatial dimension of the problem. For the one-dimensional case we consider $L = 4$ possible velocity directions, whereas for the two-dimensional case we set $L = 6$. It has been shown that increasing the number of discrete velocities $L$ does not provide better results in terms of accuracy [5] and since the computational complexity for the numerical approximation of the relaxation system clearly depends on the number of velocities, we will take those values $L = 4$ (in the 1D case) and $L = 6$ (in the 2D case) fixed from now on.

2.3.2. *Spatial discretization.* The spatial discretization of Equation (33) is performed by treating the advection term alone. Here we detail the approximation for a second-order (and fourth-order) finite volume scheme, but this can be easily extended to higher order methods. In particular, as we will observe in Section 3.2, we need a central finite volume scheme for the approximation of the last $J'$ equations.

Assuming a cartesian grid of the form $[0, L_1] \times [0, L_2]$ for the two-dimensional case, we discretize each dimension into $I_2$ volumes of size $\Delta x^2$. Since the system consists of decoupled equations along each direction, we detail the spatial discretization based on
finite volume method in the one-dimensional case, by considering a spatially uniform grid. The domain $[0, L_1]$ is divided in $I$ cells $C_i = [x_{i-1/2}, x_{i+1/2}]$, where $x_{i\pm 1/2} = x_i \pm \Delta x/2$ with constant size $\Delta x$, over which the cell average of the solution is approximated. A numerical flux function $\mathcal{F}$ is defined to approximate the continuous flux at the interface $x_{i+1/2}$ of each cell, such that the semi-discretized system reads as

$$ (34) \qquad \frac{d \bar{f}_{l,i}(t)}{dt} = -\frac{1}{\Delta x} [\mathcal{F}_{l,i+1/2}(t) - \mathcal{F}_{l,i-1/2}(t)], \quad i = 1, \dots, I, \quad l = 1, \dots, L. $$

This provides an approximation $\bar{f}_{l,i}$ of the cell average $f_{l,i}$, where the numerical flux satisfies for the linear advection with constant velocity $\Gamma_l$

$$ (35) \qquad \mathcal{F}_{l,i+1/2}(t) \approx \Gamma_l f_{l,i+1/2}(t), \quad i = 1, \dots, I, \quad l = 1, \dots, L. $$

For the first $J$ equations we consider $k$-th order upwind scheme. The numerical flux for a first order upwind scheme could be defined as

$$ (36) \qquad \mathcal{F}_{l,i+1/2}(t) = \begin{cases} \Gamma_l f_{l,i} & \text{if } \Gamma_l > 0, \\ \Gamma_l f_{l,i+1} & \text{if } \Gamma_l < 0 \end{cases} $$

Higher order extension of upwind schemes are easily recovered (see for instance [44]) and a third order upwind biased scheme reads as

$$ (37) \qquad \mathcal{F}_{l,i+1/2}(t) = \begin{cases} \Gamma_l \frac{2 f_{l,i+1} + 3 f_{l,i} - 6 f_{l,i-1} + f_{l,i-2}}{6} & \text{if } \Gamma_l > 0, \\ \Gamma_l \frac{-2 f_{l,i+2} + 6 f_{l,i+1} - 3 f_{l,i} + f_{l,i-1}}{6} & \text{if } \Gamma_l < 0, \end{cases} $$

for $l = 1, \dots, J$. The last $J'$ equations are instead approximated using a centered numerical flux, in order to guarantee the consistency of the scheme with the diffusive BGK scheme (3). The second-order centered flux scheme is given by

$$ (38) \qquad \mathcal{F}_{l,i+1/2}(t) = \Gamma_l \frac{f_{l,i+1} + f_{l,i-1}}{2}, \quad l = J + 1, \dots, L. $$

Fourth order discretization can be obtained by considering the following fluxes

$$ \mathcal{F}_{l,i+1/2}(t) = \begin{cases} \Gamma_l \frac{-f_{l,i+2} + 8 f_{l,i+1} - 8 f_{l,i-1} + f_{l,i-2}}{12} & \text{if } \Gamma_l > 0, \\ \Gamma_l \frac{f_{l,i+3} - 6 f_{l,i+2} + 18 f_{l,i+1} - 10 f_{l,i} - 3 f_{l,i-1}}{12} & \text{if } \Gamma_l < 0, \end{cases} \quad l = 1, \dots, J. $$

$$ \mathcal{F}_{l,i+1/2}(t) = \Gamma_l \frac{f_{l,i-2} - 8 f_{l,i-1} + 8 f_{l,i+1} - f_{l,i+2}}{12} \quad l = J + 1, \dots, L. $$

**Remark 2.2.** The violation of a maximum principle for a centered flux scheme could lead to unphysical oscillations, resulting in an unstable approximation of the advection equation. However, the consistency in $\varepsilon$ of the kinetic approximation (3) with the advection-diffusion equation (1) quickly stabilize the oscillations as $\varepsilon \to 0$, as already observed in [27].
<figure><img src="image_4.png" /><figcaption>FIGURE 1. Idea of Projective Integration method.</figcaption></figure>

### 3. PROJECTIVE INTEGRATION

The main idea behind the construction of projective integration methods is to combine a few small time steps of a naive (inner) time-stepping method, with a projective step (outer), in which an extrapolation over a large time interval is performed. The idea is sketched in Fig. 3.

In [27], the authors show how the choice of an outer time step much larger than the inner time step leads to a splitting of the stability region of the method in two circle-like connected components. In that case, the projective integration parameters can be taken in order to have all the fast eigenvalues in the first stability region and all the slow (dominant) eigenvalues in the dominant stability region.

The natural choice for the inner step is to use a simple, explicit method, such as a direct Euler discretization. This approach allows to obtain a fully explicit, arbitrary order in time method for stiff systems, sharing important features with asymptotic-preserving schemes. Indeed, despite the impossibility of evaluating the scheme for $\epsilon = 0$, the computational cost of this method is often independent from the stiffness of the problem. This implies that it is possible to solve stiff problems without requiring restrictive CFL condition, as done in the series of papers [6, 33, 34] for different stiff kinetic equations. The idea is to use the same approach to numerically approximate the BGK approximation (3) in the limit for $\epsilon \to 0$ with an asymptotic-preserving scheme.

Let us consider the semidiscrete version of the system (3) obtained through the phase space discretization of the form:

$$ (39) \qquad \frac{df}{dt} = D_t^\epsilon(f). $$

In our setting, $D_t^\epsilon$ is the semidiscrete operator given by

$$ (40) \qquad D_t^\epsilon(f) = -\Phi_\alpha(f) + \frac{1}{\epsilon}(M(u) - f), $$
where $\Phi_{\alpha}$ is a suitable discretization of the convective derivative, as discussed in Section 2.3.

Let us now define step sizes $\delta t$ and $\Delta t$, which represent the inner and the outer time steps respectively. Then, the solution $f(t)$ is approximated by $f^{n,k}$ at the time $t = n\Delta t + k\delta t$. As inner integrator, we choose an explicit time integration scheme $S_{\delta t}$. In particular, the inner step could be computed through the forward Euler method, obtaining

$$ (41) \qquad f^{n,k+1} = S_{\delta t}(f^{n,k}) = f^{n,k} + \delta t D_t^{\epsilon}(f^{n,k}) \quad k = 0, 1, \dots $$

The inner step has to be chosen as $\delta t = O(\epsilon)$ to ensure the stability of the method. This is problematic for $\epsilon \to 0$, since the condition becomes too restrictive and the spatial discretization in the kinetic space may lead unstabilities in the macroscopic space. When $\epsilon \to 0$, the system formally converges to a limiting equation for which every finite volume scheme needs to sastisfy a weaker stability restriction. The idea of projective integration schemes is to exploit the structure of the kinetic operator to accelerate the integration.

Starting from a computed numerical solution $f^n$ at time $t^n = n\Delta t$, one first takes $K+1$ inner steps of size $\delta t$, using the inner integrator (41) and then approximate the time derivative of $f$, computing $f^{n+1}$ via extrapolation in time. The last two updates of (41), $f^{n,K}$ and $f^{n,K+1}$ are used to project the solution forward

$$ (42) \qquad \begin{aligned} f^{n+1} &= f^{n,K+1} + (\Delta t - (K+1)\delta t) \frac{f^{n,K+1} - f^{n,K}}{\delta t} \\ &= f^{n,K+1} + M\delta t \frac{f^{n,K+1} - f^{n,K}}{\delta t}, \end{aligned} $$

in which we have defined the relative size of the extrapolation $M$ as:

$$ (43) \qquad M = \frac{\Delta t}{\delta t} - (K + 1). $$

This method (42) is called projective forward Euler (PFE), due to the fact that the extrapolation is performed using a first-order approximation of the derivative.

*High-order Projective Runge-Kutta methods.* To obtain higher-order accuracy in time we can employ any Runge-Kutta method as inner integrator, replacing each time derivative evaluation in a classical Runge-Kutta method by $K+1$ steps of an inner integrator, as discussed in [25].

However, the order of the scheme is dominated by the outer integrator. The idea is to use a particular higher order extension of projective integration, based on Runge-Kutta methods. Here, we briefly recall the construction of the projective Runge-Kutta (PRK) method. Using the same coefficients introduced in Appendix A for classical Runge-Kutta methods, for an explicit S-stage PRK scheme, we compute the values $k_s$ as

$$ (44) \qquad s = 1: \begin{cases} f^{n,k+1} = f^{n,k} + \delta t D_t(f^{n,k}) & 0 \le k \le K, \\ k_1 = \frac{f^{n,K+1} - f^{n,K}}{\delta t}, & \end{cases} $$
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Idea of the Projective Runge-Kutta Integration method. The scheme first computes $f^{n,K+1}$</figcaption></figure>

$$ (45) \qquad 2 \le s \le S : \left\{ \begin{array}{ll} f^{n+c_s,0} = f^{n,K+1} + (c_s \Delta t - (K+1)\delta t) \sum_{l=1}^{s-1} \frac{a_{s,l}}{c_s} \mathbf{k}_l, \\ f^{n+c_s,k+1} = f^{n+c_s,k} + \delta t D_t(f^{n+c_s,k}) & 0 \le k \le K, \\ \mathbf{k}_s = \frac{f^{n+c_s,K+1} - f^{n+c_s,K}}{\delta t}, & \end{array} \right. $$

where $f^{n+c_s,k} \approx f(t^n + c_s \Delta t + k \delta t)$, with $a_{s,l}$ and $c_s$ the Runge-Kutta matrix and the Runge-Kutta nodes, respectively. Finally, the numerical solution at time $t^{n+1}$ is computed as

$$ (46) \qquad f^{n+1} = f^{n,K+1} + (\Delta t - (K+1)\delta t) \sum_{s=1}^{S} b_s \mathbf{k}_s, $$

where $b_s$ are the Runge-Kutta weights. The explicit expression for the coefficients of the Runge-Kutta schemes is given in the Appendix A. In Fig. 3 we sketch the idea of Projective Runge-Kutta schemes for a second order case.

3.1. **Stability analysis.** The main advantage of projective integration approach is its great speed up for problems with great scale separation, but it is fundamental to estimate the spectrum of the inner integrator for the equation to guarantee the stability of the scheme. The estimation of the spectrum for a general BGK operator in the discrete kinetic approximation framework could be extremely hard. For this reason, we propose here an analysis of the spectrum for a simplified setting, namely taking $A(u) = u$ and $B(u) = u$ in (1) and the diagonal relaxation model (DRM1).
Let us consider the scalar BGK approximation in Equation (20) with a linear Maxwellian in the form of (21)

$$ (47) \qquad \mathcal{M}_{\lambda,\theta}(u) = \left( u + \frac{u}{\lambda} - \frac{u}{\theta^2}, u - \frac{u}{\lambda} - \frac{u}{\theta^2}, \frac{u}{\theta^2}, \frac{u}{\theta^2} \right). $$

We can rewrite the semidiscrete kinetic equation (40) in the Fourier (spatial) domain as

$$ (48) \qquad \partial_t \hat{f}(\xi_k) = \mathcal{K} \hat{f}(\xi_k), \qquad \mathcal{K} = \frac{1}{\varepsilon} (-\varepsilon D + MP - I) $$

where $\hat{f} \in \mathbb{R}^4$, $I \in \mathbb{R}^{4 \times 4}$ is the identity matrix, $P \in \mathbb{R}^{4 \times 4}$ is the orthogonal projection matrix and $M \in \mathbb{R}^{4 \times 4}$ is given by

$$ (49) \qquad M = \frac{1}{2} \begin{pmatrix} 1 + \frac{1}{\lambda} - \frac{1}{\theta^2} & 0 & 0 & 0 \\ 0 & 1 - \frac{1}{\lambda} - \frac{1}{\theta^2} & 0 & 0 \\ 0 & 0 & \frac{1}{\theta^2} & 0 \\ 0 & 0 & 0 & \frac{1}{\theta^2} \end{pmatrix}. $$

The advection matrix $D \in \mathbb{R}^{4 \times 4}$ can be written as

$$ (50) \qquad D = \begin{pmatrix} \alpha + i\beta & 0 & 0 & 0 \\ 0 & \alpha - i\beta & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{\varepsilon}}(\xi + i\gamma) & 0 \\ 0 & 0 & 0 & \frac{1}{\sqrt{\varepsilon}}(\xi - i\gamma) \end{pmatrix}. $$

Here $D$ represents the Fourier matrix of the spatial discretisation chosen for the advection term, and $\alpha, \beta$ are respectively the real and the imaginary part of the upwind approximation

$$ (51) \qquad \left\{ \begin{aligned} \alpha &= -\frac{|\lambda|}{6 \Delta x} (3 - 4 \cos(\zeta) + \cos(2\zeta)), \\ \beta &= -\frac{\lambda}{6 \Delta x} (8 \sin(\zeta) - \sin(2\zeta)), \end{aligned} \right. $$

and $\xi$ and $\gamma$ are respectively the real and the imaginary part of the centered approximation

$$ (52) \qquad \left\{ \begin{aligned} \xi &= 0, \\ \gamma &= \left( \frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\varepsilon}} \right) \frac{8 \sin(\zeta) - \sin(2\zeta)}{6 \Delta x}. \end{aligned} \right. $$

Let us define an auxiliary matrix $\mathcal{A} = -\varepsilon D + MP$, such that $\mathcal{K} = \varepsilon^{-1}(\mathcal{A} - I)$. Since the matrix has the following structure

$$ (53) \qquad \mathcal{A} = D + u e^T, $$
where e is the canonical vector in $\mathbb{R}^4$, we can use Sherman-Morrison formula to compute the determinant $\chi_A$ [17], which is given by

$$ (54) \qquad \chi_A(\zeta) = \prod_{i=1}^{4} (-\varepsilon D_i - \zeta) \left( 1 + \frac{1}{4} \sum_{j=1}^{4} \frac{M_j}{-\varepsilon D_j - \zeta} \right), $$

where $D_i$ denotes the i-th diagonal entries of $D$ and $M_j$ denotes the j-th diagonal entries of $M$.

Following [25], it is possible to localize the eigenvalues, by using Rouché's Theorem. The main differences here are given by the non-constant $\varepsilon$-dependency of the coefficients $D_j$ and the non-uniform structure of $M_j$.

**Theorem 3.1.** The spectrum of the matrix $\mathcal{A}$ satisfies

$$ (55) \qquad \mathrm{Sp}(\mathcal{A}) \subset \mathcal{D}(0, C\sqrt{\varepsilon}(|\xi| + |\gamma|)) \cup \{\zeta(\varepsilon)\}, $$

where $C$ is a constant depending on the entries of $D$, $\xi$ and $\gamma$ are defined in (50) and $\zeta(\varepsilon)$ is the dominant eigenvalue.

For the sake of completeness we retrace the proof of Theorem 3.1 in the Appendix B, extending the results in [25]. As a consequence, the spectrum of $\mathcal{K}$ satisfies

$$ (56) \qquad \mathrm{Sp}(\mathcal{K}) \subset \left( \mathcal{D}\left(-\frac{1}{\varepsilon}, C\frac{1}{\sqrt{\varepsilon}}(|\xi| + |\gamma|)\right) \right) \cup \{\zeta(\varepsilon)\}, $$

where $\zeta$ is the dominant eigenvalue, whose explicit expression, in the case of an upwind scheme (51) for the hyperbolic part and a centered scheme (52) for the parabolic part, is given by

$$ (57) \qquad \left\{ \begin{array}{l} \mathrm{Re}(\zeta(\varepsilon)) = \left(1 - \frac{1}{\theta^2}\right) \alpha + \mathcal{O}(\varepsilon), \\ \mathrm{Im}(\zeta(\varepsilon)) = \frac{\beta}{\lambda} + \mathcal{O}(\varepsilon). \end{array} \right. $$

We first observe that the spectrum of $\mathcal{K} = \varepsilon^{-1}(\mathcal{A} - I)$ consists of two well separated clusters for $\varepsilon \to 0$ and its localization on the real axis is given by two contributions: the dominant eigenvalue, whose real part is $\mathcal{O}(1)$, and the remaining eigenvalues, contained in a region around $-1/\varepsilon$. This is coherent with the derivation of the spectrum eigenvalues for a linearized BGK operator, as shown in [12].

Moreover, we assess our theoretical result by introducing a numerical approximation of the spectrum, as done in [6]. The approach is based on a discretization of the problem in phase space, obtaining a semi-discrete system (39). The set of eigenvalues of the Jacobian of $D_t^\varepsilon$ is precisely the spectrum of our problem. Introducing a finite difference approximation and implementing, for instance, the forward Euler time discretization of (39), the semidiscrete operator could be rewritten as

$$ (58) \qquad G(f^{n,k}) = \frac{F(f^{n,k}) - f^{n,k}}{\delta t}, $$

where $F(f^{n,k}) = f^{n,k+1}$. The function $G(f^{n,k})$ represents the discretization of the temporal derivative, which, according to (39), gives us a discrete of $D_t^\varepsilon$. The Jacobian of this operator and its eigenvalues can be now computed numerically, to obtain a numerical
approximation of the spectrum of our problem.

In Fig. 3 we show the spectrum of the operator for $\varepsilon = 10^{-7}$ with a central difference scheme for the parabolic terms. Similarly to the spectrum of a BGK equation in the diffusive scaling, here the imaginary part is governed by $\sqrt{\varepsilon}$, in accordance to (56).

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Spectrum of the operator for $\varepsilon = 10^{-7}$.</figcaption></figure>

In [15] it has been shown that the stability region of the PFE consists of the union of two separated disks $\mathcal{D}_1^{\text{PFE}} \cup \mathcal{D}_2^{\text{PFE}}$, where

$$ (59) \qquad \mathcal{D}_1^{\text{PFE}} = \mathcal{D} \left( 1 - \frac{\delta t}{\Delta t}, \frac{\delta t}{\Delta t} \right), \qquad \mathcal{D}_2^{\text{PFE}} \left( 0, \left( \frac{\delta t}{\Delta t} \right)^{\frac{1}{K}} \right), $$

and $\mathcal{D}(c, r)$ denotes the ball with center $(c, 0)$ and radius $r$ in the complex $\lambda$-plane. The eigenvalues lying in $\mathcal{D}_1^{\text{PFE}}$ corresponds to the slowly decaying modes, whereas the eigenvalues $\mathcal{D}_2^{\text{PFE}}$ corresponds to the fast modes, quickly damped by the inner integrator. Since the dominant set of eigenvalues is given by $\mathcal{D}_1^{\text{PFE}}$, the approach of projective integration schemes is to suitably choose the method parameters $\delta t$, $\Delta t$ and $K$ to allow for accurate integration of the modes in $\mathcal{D}_1^{\text{PFE}}$, while ensuring stability for the modes in $\mathcal{D}_2^{\text{PFE}}$.

In [25], it has been shown that the choice of stable parameters for the PFE guarantees also the stability of high order projective integration methods. Thus, according to 3.1, we set the parameters for projective integration method as follows:

$$ (60) \qquad \delta t = \varepsilon, \quad K \ge 2. $$

The $\Delta t$ is chosen according to the parabolic CFL condition.

We refer to the Appendix C for a detailed description of the stability analysis of PRK schemes.

### 3.2. Consistency analysis.
The consistency analysis for our numerical method is performed again in the scalar one-dimensional setting, following the derivation in [27]. For
this purpose, we consider the expression (8). The scheme is given by

$$ (61) \qquad f^{n+1} = S_{\delta t} f^n = f^n + \delta t \left( -\Phi^{\epsilon}(f^n) + \frac{M(u^n) - f^n}{\epsilon} \right), $$

where $\Phi^{\epsilon}$ is the spatial discretization operator.

If we consider $\delta t = \epsilon$ (according to the stability condition) and a centered flux for the last $J'$ equations, the inner scheme is then reduced to the following expression

$$ f^{n+1} = S_{\epsilon} f^n = M(u^n) - \epsilon \Phi^{\epsilon}(f^n). $$

The corresponding density $u^{N+1}$ at time $t^{N+1}$ satisfies

$$ (62) \qquad u^{N+1} = u^{N,K+1} + (\Delta t - (K+1)\epsilon) \frac{u^{N,K+1} - u^{N,K}}{\epsilon}, $$

where $K + 1$ is the number of inner time steps taken to guarantee the stability of the method.

The truncation error could be defined for the macroscopic quantities as

$$ (63) \qquad E^{N+1} = \frac{\tilde{u}^{N+1} - u^{N+1}}{\Delta t}, $$

where $\tilde{u}^N$ is the exact solution for $u$ at time $t^N$.

By considering the projective step, we can rewrite the truncation error as

$$ (64) \qquad E^{N+1} = \frac{\tilde{u}^{N+1} - u^{N,K+1}}{\Delta t} - \left( \frac{\Delta t - (K+1)\epsilon}{\Delta t} \right) \frac{u^{N,K+1} - u^{N,K}}{\epsilon}. $$

This means that the local truncation error of the projective scheme depends on its inner integrator too. Thus, we introduce the local truncation error of the inner integrator related to the kinetic variable $f$ as

$$ (65) \qquad e_f^{N,k+1} = \frac{\tilde{f}^{N,k+1} - f^{N,k+1}}{\delta t} = \frac{\tilde{f}^{N,k+1} - f^{N,k+1}}{\epsilon}. $$

We can rewrite both quantities in terms of their solution at time $t^{N,k}$. For $f^{N,k+1}$ we simply use the expression of the inner integrator, which, in this case, is chosen as the Forward Euler scheme

$$ (66) \qquad f^{N,k+1} = S_{\delta t} f^{N,k} = M(u^{N,k}) - \epsilon \Phi^{\epsilon}(f^{N,k}). $$

The exact solution $\tilde{f}^{N,k+1}$ is expanded through a Taylor series around $t^{N,k}$, using the kinetic formulation of the model

$$ (67) \qquad \begin{aligned} \tilde{f}^{N,k+1} &= \tilde{f}^{N,k} + \epsilon \partial_t \tilde{f}^{N,k} + \mathcal{O}(\epsilon^2) \\ &= \tilde{f}^{N,k} - \epsilon \left( \sum_{d=1}^{D} \Gamma_d \partial_{x_d} \tilde{f}^{N,k} \right) + \left( M(\tilde{u}^{N,k}) - \tilde{f}^{N,k} \right) + \mathcal{O}(\epsilon^2) \\ &= S_{\delta t}(\tilde{f}^{N,k}) + \epsilon \left( \Phi(\tilde{f}^{N,k}) - \sum_{d=1}^{D} \Gamma_d \partial_{x_d} \tilde{f}^{N,k} \right) + \mathcal{O}(\epsilon^2). \end{aligned} $$

The third equality is obtained from adding and subtracting $\Phi(\tilde{f}^{N,k})$ to the expression and using the forward Euler operator $S_{\delta t}$.
This leads to

$$ (68) \qquad e_f^{N,k+1} = S_{\delta t}(e_f^{N,k}) + \left( \Phi(\tilde{f}^{N,k}) - \sum_{d=1}^{D} \Gamma_d \partial_{x_d} \tilde{f}^{N,k} \right) + \mathcal{O}(\epsilon). $$

By recursion we could rewrite

$$ (69) \qquad e_f^{N,k+1} = \sum_{k=0}^{K} S_{\delta t}^{k} \left( \Phi(\tilde{f}^{N,K-k}) - \sum_{d=1}^{D} \Gamma_d \partial_{x_d} \tilde{f}^{N,K-k} \right) + \mathcal{O}((K+1)\epsilon). $$

Let us split the expression as

$$ (70) \qquad e_f^{N,k+1} = \sum_{k=0}^{K} S_{\delta t}^{k} \left( \mathcal{A}(\tilde{f}^{N,K-k}) - \sum_{d=1}^{J} \lambda_d \partial_{x_d} \tilde{f}^{N,K-k} \right) \\ \qquad \qquad + \left( \frac{\mu}{\sqrt{2}} + \frac{\theta}{\sqrt{\epsilon}} \right) \sum_{k=0}^{K} S_{\delta t}^{k} \left( \mathcal{B}(\tilde{f}^{N,K-k}) - \sum_{d=J+1}^{J'} \partial_{x_d} \tilde{f}_d^{N,K-k} \right) + \mathcal{O}((k+1)\epsilon). $$

The difference inside the first brackets in equation (70) precisely corresponds to the spatial discretization error of the hyperbolic contribution. More specifically, choosing a second order upwind scheme, we can obtain a second order spatial discretization. Unfortunately, if we compute directly the spatial truncation error between the second brackets for the parabolic term using an upwind scheme, a term $\Delta x^2/\sqrt{\epsilon}$ appears. Following the idea of [27], an approximation based on centered difference schemes allows us to recover consistency also for the parabolic term. Estimating explicitly $S_{\epsilon}\Delta_B f$, where $\Delta_B : f^{\epsilon} \in C_c^{\infty}(\Omega \times (0, T); \mathbb{R}^L) \rightarrow \mathcal{B}(f^{\epsilon}) - \sum_{d=J+1}^{J'} \partial_{x_d} \tilde{f}_d^{N,K-k}$. Then, we rewrite

$$ (71) \qquad S_{\epsilon}(\Delta_B f^{\epsilon}) = \mathcal{M} \left( \sum_i \Delta_B f_i^{\epsilon} \right) - \epsilon \mathcal{B} \Delta_B f^{\epsilon}. $$

Consequently, we have

$$ (72) \qquad S_{\epsilon}^{2}(\Delta_{B}f^{\epsilon}) = \mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}^{\epsilon}\right) - \epsilon\mathcal{M}\left(\sum_{i} \mathcal{B}\Delta_{B}f_{i}^{\epsilon}\right) - \epsilon\mathcal{B}\mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}^{\epsilon}\right) + \epsilon^{2}\mathcal{B}^{2}\Delta_{B}f^{\epsilon}, $$

which yields, for $k \ge 3$

$$ (73) \qquad S_{\epsilon}^{k}(\Delta_{B}f^{\epsilon}) = \mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}^{\epsilon}\right) - \epsilon\mathcal{M}\left(\sum_{i} \mathcal{B}\Delta_{B}f_{i}^{\epsilon}\right) - \epsilon\mathcal{B}\mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}^{\epsilon}\right) \\ \qquad \qquad + \epsilon^{2}\mathcal{M}\left(\sum_{i} \mathcal{B}^{2}\Delta_{B}f_{i}^{\epsilon}\right) + \epsilon^{2}\mathcal{B}^{2}\mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}^{\epsilon}\right) \\ \qquad \qquad + \Phi\mathcal{M}\left(\sum_{i} \mathcal{B}\Delta_{B}f_{i}\right) + \epsilon^{2}(k-3)\mathcal{M}\left(\mathcal{B}\mathcal{M}\left(\sum_{i} \Delta_{B}f_{i}\right)\right) + \mathcal{O}(\epsilon^{3}). $$

Let us now compute the average value of $S_{\epsilon}^{k}(\Delta_B f^{\epsilon})$, using the fact that $\Delta_B$ is odd.
For $k = 0$,

$$ (74) \qquad \sum_{i} \Delta_B f_i^{\epsilon} = \epsilon \sum_{i} \Delta_B g_i^{\epsilon} = \epsilon \mathcal{O}(\Delta x^2). $$
For $k = 1$,

$$ (75) \qquad \sum_i S_{\varepsilon}(\Delta_B f_i^{\varepsilon}) = \varepsilon \sum_i \Delta_B g_i^{\varepsilon} - \varepsilon \sum_i \mathcal{B} \Delta_B g_i^{\varepsilon} = \varepsilon \mathcal{O}(\Delta x^2) + \varepsilon^2 \mathcal{O}(\Delta x^2). $$

More generally, for $k \ge 2$, we obtain

$$ (76) \qquad \sum_i S_{\varepsilon}^k (\Delta_B f_i^{\varepsilon}) = \varepsilon \mathcal{O}(\Delta x^2) + \varepsilon^2 \mathcal{O}(\Delta x^2) + \mathcal{O}(\varepsilon^3). $$

The last estimate together with Young's inequality leads to

$$ (77) \qquad \sum_i e_i^{N,k+1} = \mathcal{O}\left(\sqrt{\varepsilon(k+1)\Delta x^2}\right) + \mathcal{O}\left((k+1)^{5/2}\right), \quad \forall k \ge 1. $$

This implies that

$$ E^{N+1} = \frac{\tilde{u}^{N+1} - u^{N,K+1} + \varepsilon \sum_i e_i^{N,K+1}}{\Delta t} - \left(1 - (K+1)\frac{\varepsilon}{\Delta t}\right) \frac{\tilde{u}^{N,K+1} - \tilde{u}^{N,K}}{\varepsilon} \\ - \left(1 - (K+1)\frac{\varepsilon}{\Delta t}\right) \sum_i (e_i^{N,K+1} - e_i^{N,K}) $$

$$ (78) \qquad \begin{aligned} & = \left(1 - (K+1)\frac{\varepsilon}{\Delta t}\right) \partial_t u(t^{N,K+1}) + \mathcal{O}(\Delta t) + \varepsilon^{1/2} \mathcal{O}\left(\frac{\Delta x^2 + \varepsilon^2}{\Delta t}\right) \\ & \quad - \left(1 - (K+1)\frac{\varepsilon}{\Delta t}\right) \partial_t u(t^{N,K+1}) + \mathcal{O}(\varepsilon) + \mathcal{O}\left(\frac{\varepsilon^2}{\Delta t}\right) \\ & \quad + \mathcal{O}(\varepsilon^{1/2} \Delta x^2) + \mathcal{O}(\varepsilon^{5/2}) + \mathcal{O}\left(\frac{\varepsilon^{3/2} \Delta x^2}{\Delta t}\right) + \mathcal{O}\left(\frac{\varepsilon^{7/2}}{\Delta t}\right) \end{aligned} $$

In conclusion, we summarize the above results on consistency of the method:

**Theorem 3.2.** Under the CFL condition $\Delta t = \mathcal{O}(\Delta x^2)$, the PFE scheme with second order upwind flux for $f_i$ with $i = 1, \dots, J$ and centered flux for $f_j$ with $j = J+1, \dots, L$ is consistent with (1) at order 1 in $\varepsilon$ and, as $\varepsilon \to 0$, the limiting scheme is consistent at order 1 in time and 2 in space, namely

$$ (79) \qquad E^{N+1} = \mathcal{O}(\varepsilon) + \mathcal{O}(\Delta t) + \mathcal{O}\left(\frac{\varepsilon^2}{\Delta t}\right). $$

# 4. NUMERICAL SIMULATIONS

## 4.1. Order of convergence: linear diffusion equation.
The test for the order of convergence is first performed for the DRM1 kinetic model by considering the linear diffusion problem [45]

$$ (80) \qquad \partial_t u = \xi \partial_{xx} u, $$
where $\xi = 10^{-2}$ and the initial condition is $u_0$, on $(0, T) \times \Omega = (0, T) \times [0, 1]$. In particular, we consider as initial condition

$$ (81) \qquad u_0(x) = 1 + 0.01 \exp \left( - \frac{(x - 0.5)^2}{\delta^2} \right) $$

where $\delta = 0.1$ and with periodic boundary conditions. We compare the numerical solution with the exact solution

$$ (82) \qquad u(x,t) = 1 + 0.01 \sqrt{\frac{\delta^2}{\delta^2 + 4\xi t}} \exp \left( - \frac{(x - 0.5)^2}{\delta^2 + 4\xi t} \right). $$

We first test the spatial order of convergence, where the time step is taken equal to $\Delta t = 10^{-7}$ and $\delta t = \varepsilon = 10^{-10}$. With this choice, we are able to study the decrease of the error with respect to the spatial step size $\Delta x$ (since the truncation error dueto the temporal discretization becomes negligible). In Fig 4 we compare the order of convergence in $L^1$, $L^2$ and $L^\infty$ norm of the 3rd order scheme and 4th order scheme, plotting the computed error against the expected error in log-log scale. The truncation error behaves as $O(\Delta x^3)$ and $O(\Delta x^4)$, respectively, for large $\Delta x$, until a plateau is reached. Indeed, for small values of $\Delta x$, the contribution of the inner integrator $O(\varepsilon)$ becomes dominant. This is clearly seen in the 4th order approximation, where the error does not decrease more than $10^{-10}$, which corresponds to $\varepsilon$ in this test case.

The temporal order of projective integration method is proven using a different ap-

<figure><img src="image_7.png" /><figcaption>FIGURE 4. Order of convergence for the 3rd order scheme (on the left) and the 4th order scheme (on the right) for different norms, using the linear diffusion equation (80).</figcaption></figure>

proach. Indeed, since we have to satisfy the parabolic CFL condition $\Delta t = O(\Delta x^2)$, the contribution of the spatial discretization error is dominant. Our idea is based on Richardson extrapolation [43, 22], where two solution computed with different time step $\Delta t/k$, for $k \in \mathbb{N}$ are compared.

Assuming a time discretization of order $p$, the solution $u_{\Delta t/k}^{n+1}$ is deviating from the exact
solution $u_*^{n+1}$ of the semidiscrete system (40) according to

$$ (83) \qquad u_{\Delta t/k}^{n+1} = u_{*}^{n+1} + k C_n \left(\frac{\Delta t}{k}\right)^{p+1} + \mathcal{O}(\Delta t^{p+2}), $$

where $C_n$ is the error constant. Therefore, we can rewrite the solution $u_*^{n+1}$ as

$$ (84) \qquad u_{*}^{n+1} = \frac{u_{\Delta t}^{n+1} - 2^p u_{\Delta t/2}^{n+1}}{1 - 2^p}, $$

from which we obtain an estimation of $u_{\Delta t}^{n+1} - u_{*}^{n+1}$ in terms of $u_{\Delta t}^{n+1} - u_{\Delta t/2}^{n+1}$, namely

$$ (85) \qquad u_{\Delta t}^{n+1} - u_{*}^{n+1} = \frac{2^p (u_{\Delta t}^{n+1} - u_{\Delta t/2}^{n+1})}{2^p - 1}. $$

Since to compute the order of convergence we are interested in the ratio between two "consecutive" errors, i.e. errors computed with $\Delta t/k$ and $\Delta t/(2k)$, we can give the following definition of the relative error for a general norm at time $T$

$$ (86) \qquad e_{\Delta t} = ||u_{\Delta t}(T) - u_{\Delta t/2}(T)|| $$

One observes that the truncation error behaves as expected and the results are depicted in Fig. 5, where plot the slope of the expected error n log-log scale. In this case, as for the spatial discretization error, the truncation error reaches a plateau, since the contribution of the inner integrator $\mathcal{O}(\varepsilon)$ becomes dominant.

<figure><img src="image_8.png" /><figcaption>FIGURE 5. Order of convergence for the PRK3 (on the left) and the PRK4 (on the right) for different norms, using the linear diffusion equation (80).</figcaption></figure>

4.2. **Order of convergence: advection-diffusion equation.** Following again [45], we propose a test for the advection-diffusion equation with constant velocity $c$

$$ (87) \qquad \partial_t u + c \partial_x u = \partial_{xx} u, $$

with the same initial condition as (81). The choice for the discretization parameters is analogous to the previous test. Setting $c = 10$, we reproduce the same test case, where now the solution is advected with a positive velocity. We first test again the spatial
order of convergence, in the same setting of the previous case. In Fig. 6 we compare again the order of convergence in $L^1$, $L^2$ and $L^\infty$ norm of the 3rd order and 4th order schemes, plotting the computed error against the expected error in log-log scale.

<figure><img src="image_9.png" /><figcaption>FIGURE 6. Order of convergence for the 3rd order scheme (on the left) and the 4th order scheme (on the right) for different norms, using the advection-diffusion equation (80).</figcaption></figure>

The temporal order of convergence is tested again using the definition in (86). In Fig. 7 we report the slope of the expected error in log-log scale. Differently from the previous test, the plateau is not reached here, since the spatial or the temporal error is dominant compared to the contribution of the inner integrator $O(\varepsilon)$.

<figure><img src="image_10.png" /><figcaption>FIGURE 7. Order of convergence for the PRK3 (on the left) and the PRK4 (on the right) for different norms, using the advection-diffusion equation (80).</figcaption></figure>

### 4.3. Testing the speedup: the viscous LWR equation.
Let us consider the viscous LWR equation

$$ (88) \qquad \partial_t u + \partial_x \left( \frac{u - u^2}{2} \right) = \xi \partial_{xx} u, $$
where $\xi$ assumed different values. We consider the following initial condition, as in [14]

$$ (89) \qquad u_0(x) = 0.6 + 0.25 \sin(2 \pi x), $$

in the domain $[0, 1]$ with periodic boundary conditions. The discretization is performed using $\Delta x = 5 \cdot 10^{-3}$ and $\Delta t$ according to the parabolic CFL condition. The projective integration parameters are $\delta t = \varepsilon$ and $K = 2$, where $\varepsilon = 10^{-7}$ from now on. In Fig. 8 we report the solution for $\xi = 10^{-2}$ and $\xi = 10^{-3}$. Let us now test the speedup of the numerical approximation for this test case. From a theoretical point of view, if we assume that the overhead due to extrapolation is negligible and time stepping with the innermost integrator is computationally most demanding, the speedup can be computed as

$$ (90) \qquad S_{\text{PI}} = \frac{\Delta t}{\delta t(K + 1)}, $$

that is the ratio of the total number of naive forward Euler time steps within one outermost time step over the number of real innermost steps in the method. In Table 1 we report the CPU time of the projective integration method, compared to the classical explicit integration. When $\varepsilon \to 0$, the projective integration approach becomes much more efficient than the direct integration. Moreover, in Table 2 we compare a PFE method with a first order IMEX scheme for stiff BGK equation proposed in [41]. The numerical scheme based on projective integration method exhibits a computational cost of the same order of the IMEX approach (although slightly slower) and it proves to be *asymptotic preserving* with respect to $\varepsilon$. The great advantage of the projective integration approach is that high order extension though a PRK scheme (44)-(45) is easily constructed.

<table><caption>TABLE 1. CPU time and improvement factor</caption><thead><tr><th rowspan="2">ε</th><th colspan="2">CPU time (s)</th><th colspan="2">Improvement factor</th></tr><tr><th>Direct</th><th>Projective Integration</th><th>Real</th><th>Theoretical</th></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>1.51</td><td>0.89</td><td>1.69</td><td>1.70</td></tr><tr><td>10<sup>-7</sup></td><td>152</td><td>0.94</td><td>162</td><td>171</td></tr></tbody></table>

<table><caption>TABLE 2. Comparison between the performance of PFE and IMEX schemes.</caption><thead><tr><th rowspan="2">ε</th><th colspan="2">CPU time (s)</th></tr><tr><th>Projective Integration</th><th>IMEX Scheme</th></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>0.89</td><td>0.73</td></tr><tr><td>10<sup>-7</sup></td><td>0.94</td><td>0.62</td></tr></tbody></table>

### 4.4. The viscous Burgers' equation.
Let us consider the viscous Burgers' equation

$$ (91) \qquad \partial_t u + \partial_x \left( \frac{1}{2} u^2 \right) = \xi \partial_{xx} u, $$
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Numerical approximation of the viscous LWR equation for two different values of $\xi$, respectively $\xi = 10^{-2}$ (on the left) and $\xi = 10^{-3}$ (on the right).</figcaption></figure>

where $x_i = 10^{-3}$. We first test our method for the “steady shock” whose exact solution is given by [7, 46]

$$ (92) \qquad u(x) = -\frac{2\xi}{\delta} \tanh((x - 0.5)/\delta), $$

where $\delta = 0.01$. In Fig. 9 we show the results obtained as numerical approximation with schemes of different order. The discretization is performed with step size $\Delta x = 3 \cdot 10^{-3}$ and $\Delta t$ is chosen according to the parabolic CFL condition. The projective integration parameters are chosen as $\delta t = \varepsilon$ and $K = 2$ projective levels. Compared to first-order approximation, fourth order PRK4 scheme is able to reduce numerical dissipation errors, preserving the steady state. The simulation has a perfect agreement with [46].

### 4.5. The viscous Burgers' equation with strongly degenerate diffusion.
Let us consider the Burgers equation

$$ (93) \qquad \partial_t u + \partial_x (u^2) = \xi \partial_x (\nu(u) \partial_x u), $$

with strongly degenerate diffusion

$$ (94) \qquad \nu(u) = \begin{cases} 0, & |u| \le 0.25, \\ 1, & |u| > 0.25, \end{cases} $$

and the following initial condition

$$ (95) \qquad u_0(x) = \begin{cases} 1, & \text{if } -\frac{1}{\sqrt{2}} - 0.4 \le x \le -\frac{1}{\sqrt{2}} + 0.4, \\ -1, & \text{if } \frac{1}{\sqrt{2}} - 0.4 \le x \le \frac{1}{\sqrt{2}} + 0, \\ 0, & \text{otherwise.} \end{cases} $$

We set $\xi = 0.1$ and the discretization parameters are chosen $\Delta x = 0.015$, $\delta t = \varepsilon$, $\Delta t = C \Delta x^2$ and $K = 2$. In Fig. 10 we show the numerical results at final time $t = 0.7$,
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Comparison of 1st order and 4th order approximation of the “steady shock” solution to the viscous Burgers’ equation (92).</figcaption></figure>

compared with the inviscid case ($\xi = 0$). This test shows the ability of the scheme to deal with strong degenerate diffusions.

4.6. **The two-dimensional viscous Burgers' equation.** We propose here an extension of (93) to higher dimensions. Let us focus on the 2D viscous Burgers' equation, explicitly written as

$$ (96) \qquad \partial_t u + \partial_x (u^2) + \partial_y (u^2) = \xi \partial_x (\nu(u) \partial_x u) + \xi \partial_y (\nu(u) \partial_y u), $$

again with strongly degenerate diffusion (94) and initial condition

$$ (97) \qquad u_0(x, y) = \begin{cases} -1 & \text{if } (x - 0.5)^2 + (y - 0.5)^2 \le 0.16, \\ 1 & \text{if } (x + 0.5)^2 + (y + 0.5)^2 \le 0.16, \\ 0 & \text{otherwise,} \end{cases} $$

The projective integration parameters are $K = 2$, $\delta t = \varepsilon$ and $\Delta t = \varepsilon$. The spatial discretization is obtained using $\Delta x = \Delta y = 0.03$ in the domain $[-1.5, 1.5] \times [-1.5, 1.5]$. The simulation in Fig. 11 shows the good approximation with PRK4 scheme. We underline again that a great advantage of discrete kinetic BGK approximation is its feasibility of easily extending the strategy to the multidimensional case without increasing the complexity of the structure of system (3).

4.7. **A three-phase flow model.** Let us consider the following three phase flow model introduced in the paper of Karlsen et al. [21] and approximated also in [5] to validate the
<figure><img src="image_13.png" /><figcaption>FIGURE 10. Numerical approximation of the viscous Burgers' equation compared with the inviscid solution at final time $t = 0.7$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 11. Numerical approximation of the two-dimensional viscous Burgers' equation: initial condition (left) and final solution at time $t = 0.5$ (right).</figcaption></figure>

relaxation scheme. It consists of one dimensional system of advection-diffusion equations modeling a multi-phase flow in a porous medium, which reads as follows

$$ (98) \qquad \begin{cases} \partial_t u + \partial_x f(u) = \xi \partial_{xx} B(u), \\ \partial_t v + \partial_x g(u, v) = \xi \partial_{xx} B(v), \end{cases} $$
with initial condition given by

$$ (99) \qquad (u_0(x), v_0(x)) = \begin{cases} (0.4, 0.6), & \text{if } x < 1, \\ (0, 0), & \text{otherwise.} \end{cases} $$

The variables $u$ and $v$ represent the phase saturations (gas and water, respectively), which evolves according to different dynamics. Indeed, the decoupling between the gas phase and the other phase is a reasonable assumption [21] and the flux function $(f(u), g(u, v))$ is given by

$$ (100) \qquad \left\{ \begin{aligned} f(u) &= \frac{u^2}{u^2 + (1-u)^2/10}, \\ g(u, v) &= \frac{(1-u)^2 + u^2/10}{10u^2 + (1-u)^2} f(v), \end{aligned} \right. $$

and the diffusion function is chosen such that

$$ (101) \qquad B'(w) = 4w(1-w). $$

The choice of this particular diffusion function is due to the possibility of reproducing the shape of reservoir models. Moreover, we set the scaling factor $\xi = 0.1$ and the relaxation parameter $\varepsilon = 10^{-7}$. In Fig. 12 we show the numerical results for $\Delta x = 0.005$ in [0, 2.5]. The projective integration parameters are chosen as $\delta t = \varepsilon$ and $\Delta t = C\Delta x^2$, with $K = 2$ projective levels. The results are in perfect agreement with the one shown in [5]. We observe that the choice of a 4th order scheme leads to a less diffusive approximation of the solution, as expected.

We also perform a comparison between DRM1 and DRM2 models, as shown in Fig. 13 with a zoom on the shock. For $\Delta x = 0.05$ the approximation given by the DRM2 model is less dissipative, since the choice of $\lambda_1 \neq -\lambda_2$ improves the kinetic approximation of the original system. Nevertheless, using $\Delta x = 0.0125$, we observe that the two approximations give equivalent results.

### 4.8. A three-phase flow model with gravitational effect.
We consider again the advection-diffusion Buckley-Leverett equation of the same form

$$ (102) \qquad \partial_t u + \partial_x f(u) = \xi \partial_{xx} B(u), $$

where now the flux function is given by

$$ (103) \qquad f(u) = \frac{u^2}{u^2 + a(1-u)^2} (1 - g(1-u)^2), $$

where $a$ is the ratio of the viscosities of the two fluids and $g$ is the gravitational effect. In this test the parabolic term has the form

$$ (104) \qquad B(u) = \begin{cases} 0, & \text{if } u < 0, \\ \xi \left( 2u^2 - \frac{4}{3}u^3 \right) & \text{if } 0 \le u \le 1, \\ \frac{2}{3}\xi & \text{if } u > 1. \end{cases} $$
<figure><img src="image_15.png" /><figcaption>FIGURE 12. Numerical approximation of the Buckley-Leverett equation with diffusion for the two phase saturations: gas (in blue) and water (in orange).</figcaption></figure>

Following [23], we set the following initial condition

$$ (105) \qquad u_0(x) = \begin{cases} 0, & \text{if } 0 \le x \le 1 - \frac{1}{\sqrt{2}}, \\ 1, & \text{if } \frac{1}{\sqrt{2}} < x \le 1. \end{cases} $$

We propose some numerical simulations for different values of g. We set the scaling factor $\xi = 0.01$ and the relaxation parameter $\epsilon = 10^{-7}$. The spatial numerical grid is chosen $\Delta x = 0.01$ in $[0, 1]$ and the projective integration parameters are $\delta t = \epsilon$, $\Delta t = C \Delta x^2$ and
<figure><img src="image_16.png" /><figcaption>FIGURE 13. Comparison between the DRM1 (circles) and the DRM2 (crosses) model close to the shock for the Buckley Leverett equation (98).</figcaption></figure>

K = 2. In Fig. 14 and Fig. 15 we observe the results for g = 0 and g = 5 respectively at time t = 0.05 and t = 0.1.

<figure><img src="image_17.png" /><figcaption>FIGURE 14. Numerical approximation of the Buckley-Leverett with diffusion (104) with g = 0 at two different times t = 0.05 and t = 0.1.</figcaption></figure>

4.9. **A two-dimensional three-phase model.** As last example, we consider the two-dimensional Buckley-Leverett equation with diffusion, of the form

$$ (106) \qquad \partial_t u + \partial_x f_1(u) + \partial_y f_2(u) = \xi (\partial_{xx}^2 u + \partial_{yy}^2 u), $$

where the scaling factor $\xi = 0.01$ and the flux functions are

$$ (107) \qquad f_1(u) = \frac{u^2}{u^2 + (1-u)^2}, \qquad f_2(u) = f_1(u) (1 - 5(1-u)^2). $$
<figure><img src="image_18.png" /><figcaption>FIGURE 15. Numerical approximation of the Buckley-Leverett with diffusion (104) with $g = 5$ at two different times $t = 0.05$ and $t = 0.1$.</figcaption></figure>

The initial condition is given by

$$ (108) \qquad u(x, y, 0) = \begin{cases} 1, & \text{if } x^2 + y^2 < 0.5, \\ 0, & \text{otherwise.} \end{cases} $$

The domain is $[-1.5, 1.5] \times [-1.5, 1.5]$ and the spatial discretization is obtained taking $\Delta x = \Delta y = 0.015$. The parameters for the projective integration method are the same as before, namely $\delta t = \varepsilon$, $\Delta t = C\Delta x^2$ and $K = 2$, where $C = 0.4$. Since this simulation has a discontinuous initial condition in two dimensions which may lead to oscillatory behavior, we employ a CWENO3 approximation for the hyperbolic variables, following the scheme [29]. The results at time $t = 0.5$ agree with the ones of [24]. We underline again that the advantage of using relaxation schemes relies on their easily extensibility to the multi-dimensional case, without increasing the complexity of the numerical scheme.

## 5. CONCLUSION

In this paper, we have proposed an explicit scheme for the solution of discrete kinetic approximations of advection-diffusion equations. Indeed, the stiffness due to the presence of a relaxation operator could be treated with a projective integration approach, in order to avoid implicit time integration. The idea of using projective integration methods for discrete kinetic approximation could be easily extended to other types of equations. Indeed, relaxation schemes have a deep connection with lattice Boltzmann Methods (LBM), see for instance [18, 46, 2, 3]. In a very recent paper, the theory of discrete kinetic BGK approximation of hyperbolic problems has been used to prove monotonicity properties for the LBM with overrelaxation [2]. This approach could be extended to the case of possibly degenerate parabolic equations and the computational performance of LBM could be compared (or coupled) with the projective integration
<figure><img src="image_19.png" /><figcaption>FIGURE 16. Numerical approximation of the Buckley-Leverett with diffusion (106) at final time $t = 0.5$.</figcaption></figure>

approach.

**Acknowledgement.** The author would like to thank Thomas Rey and Giuseppe Visconti for comments on this manuscript and for their careful reading. He would also like to thank Roberto Natalini and Giovanni Samaey for fruitful discussions on this topic. The author is funded by the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No.101072546) and he is member of GNCS-INdAM research group.

**Data availability.** Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

## DECLARATIONS

**Conflict of interest.** The author declares that he has no financial or non-financial conflicts of interest that are relevant to the content of this article.

## REFERENCES

[1] ARBOGAST, T., WHEELER, M. F., AND ZHANG, N.-Y. A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media. *SIAM J. Numer. Anal.* **33**, 4 (1996), 1669–1687.

[2] AREGBA-DRIOLLET, D. Convergence of lattice Boltzmann methods with overrelaxation for a non-linear conservation law. *ESAIM, Math. Model. Numer. Anal.* **58**, 5 (2024), 1935–1958.

[3] AREGBA-DRIOLLET, D., AND BELLOTTI, T. Monotonicity and convergence of two-relaxation-times lattice Boltzmann schemes for a non-linear conservation law. Preprint, arXiv:2501.07934 [math.NA] (2025), 2025.

[4] AREGBA-DRIOLLET, D., AND NATALINI, R. Discrete kinetic schemes for multidimensional systems of conservation laws. *SIAM J. Numer. Anal.* **37**, 6 (2000), 1973–2004.

[5] AREGBA-DRIOLLET, D., NATALINI, R., AND TANG, S. Explicit diffusive kinetic schemes for nonlinear degenerate parabolic systems. *Math. Comput.* **73**, 245 (2004), 63–94.

[6] BAILO, R., AND REY, T. Projective and telescopic projective integration for non-linear kinetic mixtures. *J. Comput. Phys.* **458** (2022), 24.
[7] BENTON, E. R., AND PLATZMAN, G. W. A table of solutions of the one-dimensional Burgers equation. *Q. Appl. Math.* **30** (1972), 195–212.

[8] BESSEMOULIN-CHATARD, M., AND FILBET, F. A finite volume scheme for nonlinear degenerate parabolic equations. *SIAM J. Sci. Comput.* **34**, 5 (2012), b559–b583.

[9] BOSCARINO, S., PARESCHI, L., AND RUSSO, G. *Implicit-explicit methods for evolutionary partial differential equations (to appear)*, vol. 24 of *Math. Model. Comput. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM)*, 2024.

[10] BOUCHUT, F., GUARGUAGLINI, F. R., AND NATALINI, R. Diffusive BGK approximations for non-linear multidimensional parabolic equations. *Indiana Univ. Math. J.* **49**, 2 (2000), 723–749.

[11] CAVALLI, F., NALDI, G., PUPPO, G., AND SEMPLICE, M. High-order relaxation schemes for non-linear degenerate diffusion problems. *SIAM J. Numer. Anal.* **45**, 5 (2007), 2098–2119.

[12] CERCIGNANI, C. *The Boltzmann equation and its applications*, vol. 67 of *Appl. Math. Sci.* New York etc.: Springer-Verlag, 1988.

[13] DIAZ, J. I., AND KERSNER, R. On a nonlinear degenerate parabolic equation in infiltration or evaporation through a porous medium. *J. Differ. Equations* **69** (1987), 368–403.

[14] FERRETTI, R. A fully semi-Lagrangian technique for viscous and dispersive conservation laws. *J. Comput. Phys.* **526** (2025), 113784.

[15] GEAR, C. W., AND KEVREKIDIS, I. G. Projective methods for stiff differential equations: Problems with gaps in their eigenvalue spectrum. *SIAM J. Sci. Comput.* **24**, 4 (2003), 1091–1106.

[16] GEAR, C. W., AND KEVREKIDIS, I. G. Telescopic projective methods for parabolic differential equations. *J. Comput. Phys.* **187**, 1 (2003), 95–109.

[17] GOLUB, G., AND VAN LOAN, C. F. *Matrix computations.*, 3rd ed. ed., vol. 43 of *Texts Read. Math. New Delhi: Hindustan Book Agency*, 2007.

[18] GRAILLE, B. Approximation of mono-dimensional hyperbolic systems: a lattice Boltzmann scheme as a relaxation method. *J. Comput. Phys.* **266** (2014), 74–88.

[19] JIN, S., PARESCHI, L., AND TOSCANI, G. Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations. *SIAM J. Numer. Anal.* **35**, 6 (1998), 2405–2439.

[20] JIN, S., AND XIN, Z. The relaxation schemes for systems of conservation laws in arbitrary space dimensions. *Commun. Pure Appl. Math.* **48**, 3 (1995), 235–276.

[21] KARLSEN, K. H., LIE, K.-A., NATVIG, J. R., NORDHAUG, H. F., AND DAHLE, H. K. Operator splitting methods for systems of convection-diffusion equations: nonlinear error mechanisms and correction strategies. *J. Comput. Phys.* **173**, 2 (2001), 636–663.

[22] KOELLERMEIER, J., AND SAMAEY, G. Projective integration methods in the Runge-Kutta framework and the extension to adaptivity in time. *J. Comput. Appl. Math.* **454** (2025), 26.

[23] KOSSACZKÁ, T., EHRHARDT, M., AND GÜNTHER, M. A neural network enhanced weighted essentially non-oscillatory method for nonlinear degenerate parabolic equations. *Physics of Fluids* **34**, 2 (2022), 026604.

[24] KURGANOV, A., AND TADMOR, E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. *J. Comput. Phys.* **160**, 1 (2000), 241–282.

[25] LAFITTE, P., LEJON, A., AND SAMAEY, G. A high-order asymptotic-preserving scheme for kinetic equations using projective integration. *SIAM J. Numer. Anal.* **54**, 1 (2016), 1–33.

[26] LAFITTE, P., MELIS, W., AND SAMAEY, G. A high-order relaxation method with projective integration for solving nonlinear systems of hyperbolic conservation laws. *J. Comput. Phys.* **340** (2017), 1–25.

[27] LAFITTE, P., AND SAMAEY, G. Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit. *SIAM J. Sci. Comput.* **34**, 2 (2012), a579–a602.

[28] LATTANZIO, C., AND NATALINI, R. Convergence of diffusive BGK approximations for nonlinear strongly parabolic systems. *Proc. R. Soc. Edinb., Sect. A, Math.* **132**, 2 (2002), 341–358.

[29] LEVY, D., PUPPO, G., AND RUSSO, G. A third order central WENO scheme for 2d conservation laws. *Appl. Numer. Math.* **33**, 1-4 (2000), 415–421.

[30] LIONS, P. L., PERTHAME, B., AND TADMOR, E. A kinetic formulation of multidimensional scalar conservation laws and related equations. *J. Am. Math. Soc.* **7**, 1 (1994), 169–191.
[31] LIONS, P. L., AND TOSCANI, G. Diffusive limit for finite velocity Boltzmann kinetic models. *Rev. Mat. Iberoam.* **13**, 3 (1997), 473–513.

[32] LIU, Y., SHU, C.-W., AND ZHANG, M. High order finite difference WENO schemes for nonlinear degenerate parabolic equations. *SIAM J. Sci. Comput.* **33**, 2 (2011), 939–965.

[33] MELIS, W., REY, T., AND SAMAEY, G. Projective and telescopic projective integration for the nonlinear BGK and Boltzmann equations. *SMAI J. Comput. Math.* **5** (2019), 53–88.

[34] MELIS, W., AND SAMAEY, G. Telescopic projective integration for linear kinetic equations with multiple relaxation times. *J. Sci. Comput.* **76**, 2 (2018), 697–726.

[35] NALDI, G., AND PARESCHI, L. Numerical schemes for kinetic equations in diffusive regimes. *Appl. Math. Lett.* **11**, 2 (1998), 29–35.

[36] NALDI, G., AND PARESCHI, L. Numerical schemes for hyperbolic systems of conservation laws with stiff diffusive relaxation. *SIAM J. Numer. Anal.* **37**, 4 (2000), 1246–1270.

[37] NALDI, G., PARESCHI, L., AND TOSCANI, G. Relaxation schemes for partial differential equations and applications to degenerate diffusion problems. *Surv. Math. Ind.* **10**, 4 (2002), 315–343.

[38] NATALINI, R. Convergence to equilibrium for the relaxation approximations of conservation laws. *Commun. Pure Appl. Math.* **49**, 8 (1996), 795–823.

[39] NATALINI, R. A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws. *J. Differ. Equations* **148**, 2 (1998), 292–317.

[40] NATALINI, R., RIBOT, M., AND TWAROGOWSKA, M. A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis. *J. Sci. Comput.* **63**, 3 (2015), 654–677.

[41] PARESCHI, L., AND RUSSO, G. Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. *J. Sci. Comput.* **25**, 1-2 (2005), 129–155.

[42] POMRANING, G. C., AND FGLESONG, G. M. Transport-diffusion interfaces in radiative transfer. *J. Comput. Phys.* **32** (1979), 420–436.

[43] RICHARDSON, L. F. IX. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. *Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character* **210**, 459-470 (1911), 307-357.

[44] SHU, C.-W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In *Advanced numerical approximation of nonlinear hyperbolic equations. Lectures given at the 2nd session of the Centro Internazionale Matematico Estivo (C. I. M. E.) held in Cetraro, Italy, June 23-28, 1997*. Berlin: Springer, 1998, pp. 325–432.

[45] WISSOCQ, G., AND ABGRALL, R. A new local and explicit kinetic method for linear and non-linear convection-diffusion problems with finite kinetic speeds. I: One-dimensional case. *J. Comput. Phys.* **518** (2024), 29.

[46] WISSOCQ, G., LIU, Y., AND ABGRALL, R. A positive- and bound-preserving vectorial lattice Boltzmann method in two dimensions. Preprint, arXiv:2411.15001 [math.NA] (2024), 2024.

[47] ZHANG, Q., AND WU, Z.-L. Numerical simulation for porous medium equation by local discontinuous Galerkin finite element method. *J. Sci. Comput.* **38**, 2 (2009), 127–148.

## APPENDIX A. RUNGE-KUTTA METHODS

In Section 3 we have presented a higher order extension of projective integration schemes based on Runge-Kutta (RK) methods. Let us recall the main properties of RK methods. We consider an explicit S-stage RK method for a general equation of the form

$$ (109) \qquad y' = f(t, y), \quad y(t_0) = y_0. $$

Let us consider a time step discretization given by $\Delta t$, then we have

$$ (110) \qquad \begin{cases} f^{n+c_s} = f^n + \Delta t \sum_{l=1}^{s-1} a_{s,l} \mathbf{k}_l, \\ \mathbf{k}_s = f(t^n + c_s \Delta t, f^{n+c_s}) \end{cases} \quad 1 \le s \le S. $$
The matrix $\mathbf{A} = (a_{s,l})_{s,l=1}^S$ is the so-called Runge-Kutta matrix, the vector $\mathbf{b} = (b_s)_{s=1}^S$ and $\mathbf{c} = (c_s)_{s=1}^S$ are the Runge-Kutta weights and the Runge-Kutta nodes, respectively. The weights $\mathbf{b}_s$ and $\mathbf{c}_s$ are chosen to ensure consistency of the method, requiring $0 \le b_s \le 1$ and $0 \le c_s \le 1$ and

$$ (111) \qquad \sum_{s=1}^{S} b_s = 1, \sum_{s=1}^{S-1} a_{s,l} = c_s, \quad 1 \le s \le S. $$

Here we represent the Bucther's tableaux used in this paper, where we first recall the general form:

$$ \begin{array}{c|c} c & A \\ \hline & b^T \end{array} $$

Now, we give detailed expressions of the Runge-Kutta tableaux for the third and the fourth order methods.

* Third-order SSP Runge-Kutta

$$ \begin{array}{c|ccc} 0 & 1 & 1 \\ 1 & 1/2 & 1/4 \\ \hline & 1/2 & 1/4 \\ & & 1/6 \end{array} \begin{array}{ccc} 1/4 & 1/6 & 2/3 \end{array} $$

* Classical fourth-order Runge-Kutta method

$$ \begin{array}{c|cccc} 0 & 1/2 & 1/2 & 1/2 \\ 1/2 & 1/2 & 0 & 1/2 \\ \hline & 1/2 & 0 & 0 \\ & 1 & 0 & 1 \end{array} \begin{array}{cccc} 1/6 & 1/3 & 1/3 & 1/6 \end{array} $$

## APPENDIX B. STABILITY OF THE INNER INTEGRATORS

In this section we prove the result on the spectrum of the inner integrator, partially retracing the proof in [25]. The main difference are the non-constant $\varepsilon$-dependent coefficient in front of the entries $D_j$ and the values of the Maxwellians $M_j$.

**Theorem 3.1.** The spectrum of the matrix $\mathcal{A}$ defined in Section 3.1 satisfies

$$ (112) \qquad \mathrm{Sp}(\mathcal{A}) \subset \mathcal{D}(0, C\sqrt{\varepsilon}(|\xi| + |\gamma|)) \cup \{\zeta(\varepsilon)\}, $$

where $C$ is a constant depending on the entries of $D$, $\xi$ and $\gamma$ are defined in (50) and $\zeta(\varepsilon)$ is the dominant eigenvalue.

*Proof.* In the first step we use Rouché's Theorem to localize the eigenvalues, whereas in the second step we give the explicit asymptotic expansion for the dominant eigenvalue. We have obtained the following expression for the characteristic polynomial of $\mathcal{A}$,

$$ (113) \qquad \chi_{\mathcal{A}}(\zeta) = \prod_{i=1}^{4} (-\varepsilon D_i - \zeta) \left(1 - \frac{1}{4} \sum_{j=1}^{4} \frac{M_j}{\zeta - \varepsilon D_j}\right). $$
The expression for the characteristic polynomial of matrix $\mathcal{A}_0 := D + P$ is given instead by

$$ (114) \qquad \chi_{\mathcal{A}_0}(\zeta) = \prod_{i=1}^{4} (-\varepsilon D_i - \zeta) \left( 1 - \frac{1}{4} \sum_{j=1}^{4} \frac{1}{\zeta - \varepsilon D_j} \right). $$

The idea is to study the behavior of the rational function

$$ (115) \qquad \mathcal{F}: \zeta \to \frac{\chi_{\mathcal{A}(\zeta)} - \chi_{\mathcal{A}_0}(\zeta)}{\chi_{\mathcal{A}_0}(\zeta)} = - \frac{\frac{1}{4} \sum_{j=1}^{4} (M_j - 1) \frac{1}{\zeta - \varepsilon D_j}}{1 - \frac{1}{4} \sum_{j=1}^{4} \frac{1}{\zeta - \varepsilon D_j}} $$

and to observe the regions containing eigenvalues of $\mathcal{A}_0$ for which $|\mathcal{F}(\zeta)| < 1$.

Here, we just point out some aspects of the computation which differ from the original one.

* Let us first consider the dominant eigenvalue of $\mathcal{A}_0$, restricting the analysis to the region parametrized by $\Omega_1 := \{1 + \varepsilon r e^{i\theta}, \theta \in [0, 2\pi)\}$, performing a Taylor expansion of $\frac{1}{1 + \varepsilon r e^{i\theta} - \varepsilon D_j}$ in terms of $\varepsilon$:

$$ \frac{1}{1 + \varepsilon r e^{i\theta} - \varepsilon D_j} = 1 + (\varepsilon D_j - \varepsilon r e^{i\theta}) + \mathcal{O}(\varepsilon^2). $$

We observe that

$$ (116) \qquad \sum_{j=1}^{4} (M_j - 1) = 0. $$

The rational function $\mathcal{F}$ could be expanded as

$$ (117) \qquad \mathcal{F}(1 + \varepsilon r e^{i\theta}) = - \frac{\frac{1}{4} \sum_{j=1}^{4} (M_j - 1) D_j + \mathcal{O}(\varepsilon)}{r e^{i\theta} - \frac{1}{4} \sum_{j=1}^{4} D_j + \mathcal{O}(\varepsilon)}. $$

Thus, choosing $r > \frac{1}{4} \left| \sum_j D_j \right| + \frac{1}{2} \left| \sum_j D_j (M_j - 1) \right|$, we ensure that $|\mathcal{F}(\zeta)| < 1/2 + \mathcal{O}(\varepsilon)$. From Rouché's Theorem, in $\Omega_1$ (around $\zeta = 1$), $\chi_A$ and $\chi_{\mathcal{A}_0}$ have the same number of zeroes, which means that $\mathcal{A}$ have only one eigenvalue in $\Omega_1$.

* We focus now on the remaining eigenvalues, by considering the region around $\zeta = 0$, which contains the “fast” eigenvalues of $\mathcal{A}_0$. Then, we restrict the analysis to the region $\Omega_2 := \{\varepsilon r e^{i\theta}, \theta \in [0, 2\pi)\}$ and we perform a Taylor expansion of
$\frac{1}{r e^{i\theta} - D_j}$ in terms of $1/r$:

$$ \frac{1}{r e^{i\theta} - D_j} = \frac{e^{-i\theta}}{r} + \frac{D_j e^{-i2\theta}}{r^2} + \mathcal{O}\left(\frac{1}{r^3}\right). $$

Using again (116), the rational function $\mathcal{F}$ could be now expanded as

$$ \begin{aligned} \mathcal{F}(\varepsilon r e^{i\theta}) = & -\frac{\frac{1}{4} \sum_{j=1}^{4} (M_j - 1) D_j \frac{e^{-i2\theta}}{r^2} + \mathcal{O}\left(\frac{1}{r^3}\right)}{4\varepsilon - \frac{e^{-i\theta}}{r} + \mathcal{O}\left(\frac{1}{r^2}\right)} \\ & -\frac{\frac{1}{r} \frac{e^{-i\theta}}{4} \frac{1}{4} \sum_{j=1}^{4} (M_j - 1) D_j + \mathcal{O}\left(\frac{1}{r}\right)}{4\varepsilon r e^{i\theta} - 1 + \mathcal{O}\left(\frac{1}{r}\right)}. \end{aligned} $$

We finally perform a Taylor expansion in terms of $\varepsilon$, to obtain

$$ (118) \qquad \mathcal{F}(\varepsilon r e^{i\theta}) = \frac{e^{-i\theta}}{r} \frac{1}{4} \sum_{j=1}^{4} (M_j - 1) D_j + \mathcal{O}(\varepsilon). $$

Again, choosing $r > \frac{1}{2} \left| \sum_j D_j \right| + 1$, we ensure the inequality $|\mathcal{F}(\zeta)| < 1/2$ and we conclude using Rouché's Theorem that in $\Omega_2$ (around $\zeta = 0$), $\chi_A$ and $\chi_{A_0}$ have the same number of zeroes, which means that $A$ have $J - 1$ eigenvalues in $\Omega_2$.

In the second step we give an asymptotic expression for the largest eigenvalue in terms of $\varepsilon$. We start from the expression of $\chi_A$ given in (113). If $\zeta(\varepsilon)$ is a root of $\chi_A$, it holds

$$ (119) \qquad \tilde{\chi}_A(\zeta) := 1 - \frac{1}{4} \sum_{j=1}^{4} \frac{M_j}{\zeta - \varepsilon D_j} = 0. $$

Let us consider an expression of $\zeta$ in terms of its real coordinates, namely

$$ (120) \qquad \zeta(\varepsilon) = x(\varepsilon) + i y(\varepsilon) $$

and let us consider the asymptotic expansion of $x(\varepsilon)$ and $y(\varepsilon)$ in terms of $\varepsilon$:

$$ x(\varepsilon) = x_0 + \varepsilon x_1 + \varepsilon^2 x_2 + \mathcal{O}(\varepsilon^3), \qquad y(\varepsilon) = y_0 + \varepsilon y_1 + \varepsilon^2 y_2 + \mathcal{O}(\varepsilon^3). $$

We plug these expansions in (119) and we perform a Taylor series expansion in terms of $\varepsilon$. Since $\zeta(\varepsilon)$ is a root for $\tilde{\chi}_A$, all the coefficients of equal power in $\varepsilon$ must vanish. For the zeroth order term we have

$$ (121) \qquad 1 - \frac{1}{4} \sum_{j=1}^{4} \frac{M_j}{x_0 + i y_0} = 0, $$

which, since $\sum_j (M_j - 1) = 0$ implies

$$ \frac{x_0}{x_0^2 + y_0^2} = 1 \qquad - \frac{y_0}{x_0^2 + y_0^2} = 0, $$
namely $x_0 = 1$ and $y_0 = 0$. Let us now determine the terms of order $\epsilon$, using

$$ \frac{1}{4} \sum_{j=1}^{4} M_j [(x_1 + i y_1) + D_j] = 0, $$

which implies

$$ x_1 = \left(1 - \frac{1}{\theta^2}\right) \alpha + \frac{1}{\theta^2 \sqrt{\epsilon}} \xi, \quad y_1 = \frac{1}{\lambda} \beta. $$

Finally, let us consider the second order terms in $\epsilon^2$

$$ \frac{1}{4} \sum_{j=1}^{4} M_j \left[ (x_1 + i y_1)^2 + 2 D_j^2 + 4(x_1 + i y_1) D_j - 2(x_2 + i y_2) \right] = 0, $$

for which we obtain

$$ x_2 = 2(x_1^2 - y_1^2) + 2 \left(1 - \frac{1}{\theta^2}\right) (\alpha^2 - \beta^2) + \frac{2}{\epsilon\theta^2} (\xi^2 - \gamma^2), $$

$$ y_2 = 4x_1 y_1 + 2 \frac{1}{\lambda} \alpha \beta, $$

□

## APPENDIX C. STABILITY OF HIGHER ORDER PROJECTIVE INTEGRATION

Let us recall here the stability properties of the projective Runge-Kutta methods. In order to compute the linear stability regions, the Dahlquist equation is introduced

$$ (122) \qquad \dot{y} = \lambda y, \quad \mathrm{Re}(\lambda) < 0, $$

whose general solution does not blow up in time. Introducing a general inner integrator for this equation

$$ (123) \qquad y^{k+1} = \tau(\lambda\delta t)y^k, $$

we call $\tau(\lambda\delta t)$ the *amplification factor* of the inner integrator. Applying explicit time integration to equation (122), we observe that

$$ (124) \qquad y^{k+1} = \tau^{k+1} y_0, $$

which means that the inner integrator is linearly stable if $|\tau| \le 1$. In the case of forward Euler integrator, we have $\tau(\lambda\delta t) = 1 + \lambda\delta t$. The goal of this section is to determine the values of $\tau$ for which the projective integration method is also stable. Let us first focus on projective Forward-Euler method applied to (122), which can be written as

$$ (125) \qquad y^{n+1} = \sigma(\tau, \Delta t, \delta t, K) y^n, $$

where the amplification factor $\sigma$ is given by

$$ (126) \qquad \sigma^{\mathrm{PFE}}(\tau, \Delta t, \delta t, K) = \left[ \left( \frac{\Delta t - (K+1)\delta t}{\delta t} + 1 \right) \tau - \frac{\Delta t - (K+1)\delta t}{\delta t} \right] \tau^k. $$

The PFE method is stable if $|\sigma^{\mathrm{PFE}}(\tau, \Delta t, \delta t, K)| \le 1$. The main goal of introducing a projective integration approach is to use a time step $\Delta t = O(\Delta x^\gamma)$. In the limit $\epsilon \to 0$ for fixed $\Delta x$, we look at the limiting stability regions as $\Delta t/\delta t \to \infty$. In [15] it has
been shown that the values $\tau$ satisfying (126) lie in the union of two separated disks $D_1^{\text{PFE}} \cup D_2^{\text{PFE}}$, where

$$ (127) \qquad D_1^{\text{PFE}} = D \left( 1 - \frac{\delta t}{\Delta t}, \frac{\delta t}{\Delta t} \right), \qquad D_2^{\text{PFE}} \left( 0, \left( \frac{\delta t}{\Delta t} \right)^{\frac{1}{K}} \right), $$

and $D(c, r)$ denotes the ball with center $(c, 0)$ and radius $r$ in the complex $\lambda$-plane. For higher order projective Runge-Kutta (PRK) schemes, it is possible to give an explicit expression of the amplification factor, applying the integrator to Equation (122). For the sake of simplicity, we introduce

$$ (128) \qquad M = \frac{\Delta t - (K + 1)\delta t}{\delta t}, \quad \text{and} \quad M_s = \frac{c_s \Delta t - (K + 1)\delta t}{\delta t}, $$

where $c_s$ are the RK nodes. We also denote by $\mathbf{A} = (a_{s,l})_{s=1}^S$ and $\mathbf{b} = (b_s)_{s=1}^S$ the RK matrix and the RK weights for a general $S$-stage Runge-Kutta method. Then, for an explicit PRK method, we have

$$ (129) \qquad \left\{ \begin{array}{l} \mathbf{k}_1 = \kappa_1(\tau)y^n = \frac{\tau^{K+1} - \tau^K}{\delta t} y^n, \\[1em] \mathbf{k}_s = \kappa_s(\tau)y^n = \frac{\tau^{K+1} - \tau^K}{\delta t} \left( \tau^{K+1} + (M_s \delta t) \sum_{l=1}^{s-1} \frac{a_{s,l}}{c_s} \kappa_l \right) y^n, \quad 2 \le s \le S, \\[1em] y^{n+1} = \sigma(\tau)y^n = \left( \tau^{K+1} + (M \delta t) \sum_{s=1}^S b_s \kappa_s \right) y^n. \end{array} \right. $$

The following result on the stability of higher order projective Runge-Kutta methods holds.

**Theorem C.1** ([25], Theorem 3.2). Assume the inner integrator is stable, i.e. $|\tau| \le 1$ and $\delta t, \Delta t$ and $K$ are chosen such that the projective forward Euler method is stable. Then, a projective Runge-Kutta method is also stable if it satisfies the following assumption

$$ (130) \qquad \sum_{s=1}^{S} b_s = 1, \quad \sum_{s=1}^{S-1} a_{s,l} = c_s, \quad 1 \le s \le S, $$

with $0 \le b_s, c_s \le 1$ and the convexity condition

$$ (131) \qquad 0 \le a_{s,l} \le c_s, \quad \forall 1 \le l \le s, \quad \forall 1 \le s \le S. $$

In Fig. 17-18 we show the stability regions of projective Runge-Kutta schemes for different values of $\delta t/\Delta t$ and $K$. The main feature, as proven in Theorem C.1, is that the qualitative dependence of higher order PRK methods on the parameters is equivalent to those of PFE.
<figure><img src="image_20.png" /><figcaption>FIGURE 17. Stability region for higher order projective integrators in the complex $h\lambda$-plane for $\Delta t/\delta t = 10$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 18. Stability region for higher order projective integrators in the complex $h\lambda$-plane for $\Delta t/\delta t = 25$.</figcaption></figure>