On the non-global linear stability and spurious fixed points of MPRK schemes with negative RK parameters

Recently, a stability theory has been developed to study the linear stability of modified Patankar--Runge--Kutta (MPRK) schemes. This stability theory provides sufficient conditions for a fixed point of an MPRK scheme to be stable as well as for the convergence of an MPRK scheme towards the steady state of the corresponding initial value problem, whereas the main assumption is that the initial value is sufficiently close to the steady state. Initially, numerical experiments in several publications indicated that these linear stability properties are not only local, but even global, as is the case for general linear methods. Recently, however, it was discovered that the linear stability of the MPDeC(8) scheme is indeed only local in nature. Our conjecture is that this is a result of negative Runge--Kutta (RK) parameters of MPDeC(8) and that linear stability is indeed global, if the RK parameters are nonnegative. To support this conjecture, we examine the family of MPRK22($\alpha$) methods with negative RK parameters and show that even among these methods there are methods for which the stability properties are only local. However, this local linear stability is not observed for MPRK22($\alpha$) schemes with nonnegative Runge-Kutta parameters.


Introduction
Modified Patankar-Runge-Kutta (MPRK) schemes are numerical time integration schemes which are conservative and unconditionally positivity-preserving, when applied to a positive and conservative production-destruction system (PDS) where y = (y 1 , . . ., y N ) T and p ij (t, y), d ij (t, y) ≥ 0 for all y > 0, t ≥ 0. A PDS is called positive, when positive initial values imply positive solutions for all times, and is called conservative, when 1 T y ′ = 0 or equivalently 1 T y is constant for all times.
Here, 1 = (1, . . ., 1) T ∈ R N denotes the vector with all elements equal to one.Since MPRK schemes do not belong to the class of general linear methods, they can be unconditionally positivity-preserving and of higher order at the same time.For general linear methods unconditional positivity is restricted to first order schemes, see [1][2][3].So far, second and third order MPRK schemes for the integration of autonomous PDS have been constructed in [4][5][6], based on the classical Runge-Kutta (RK) formulation.
The same was done in [7,8], based on the SSP formulation of RK schemes.In [9] arbitrary high order MPRK schemes on the basis of deferred correction schemes have been introduced.MPRK schemes for time dependent PDS have been investigated in [10].For other approaches to obtain unconditional positivity we refer to [2,[11][12][13][14][15][16].
We also want to mention the novel framework [17] on order conditions for Runge-Kutta-like schemes to which MPRK methods belong.Besides a general theory for deriving order conditions, for the first time sufficient and necessary conditions for fourth order MPRK schemes were presented and reduced therein.Additionally, the authors proved that the order of an MPRK method does not depend on the signs of the RK parameters.
In addition to the order of a numerical scheme, its stability is of course crucial for its usefulness in practical applications.In the following, we will be concerned with linear stability of MPRK schemes, i. e. the stability behavior of MPRK schemes when applied to a positive and conservative linear system.Based on A = (a ij ) ∈ R N ×N , a positive and conservative linear system has the form where A − diag(A) ≥ 0 is necessary and sufficient for positivity, see [18], and 1 T A = 0 must hold to ensure conservativity.Here, the notation diag(A) ∈ R N ×N is used to denote the diagonal matrix with diagonal elements equal to those of A. Furthermore, the conservation property 1 T A = 0 together with positivity implies diag(A) ≤ 0. If we define B = A − diag(A), system (2) can be rewritten in production-destruction form as (3) Since B = (b ij ) ≥ 0 and y ≥ 0 by assumption the system's production terms are p ij (y) = b ij y j = a ij y j ≥ 0 for i ̸ = j and the corresponding destruction terms d ij are contained within − diag(A)y ≥ 0. For N = 2 all positive and conservative linear systems (3) can be represented as The central notion for linear stability of general linear methods is A-stability.A general linear method is said to be A-stable, whenever its numerical solution of y ′ = λy for an arbitrary time step size ∆t > 0 tends to zero for all λ ∈ C − = {z ∈ C | Re(z) < 0}.The choice λ ∈ C − comes from the fact, that A-stability ensures that the numerical solution of a linear system y ′ = Ay tends to 0, whenever all eigenvalues of A belong to C − .For an A-stable linear method applied to y ′ = Ay with a spectrum σ(A) ⊆ C − , the unique steady state solution y * = 0 is an asymptotically stable fixed point for arbitrary time step sizes.However, the crucial difference between A-stability and asymptotic stability is that the former is a global property, i.e., independent of the initial value, while the latter is a local property, since asymptotic stability requires that the initial value be sufficiently close to the fixed point.
Due to its importance for general linear methods, we would also like to investigate the A-stability of MPRK schemes.Unfortunately, several obstacles stand in the way.First, MPRK schemes cannot be applied to the scalar linear test equation, particularly not with λ ∈ C − , as it is unclear how the complex term λy can be split into production and destruction terms.To get around this, we can apply MPRK schemes directly to positive and conservative linear systems.But, the conservation property is in contradiction with the asymptotic stability of y * = 0, i. e. there is no conservative linear system whose solution tend to 0 for t → ∞ for initial values y 0 > 0, since 1 T y must be constant for all times.Moreover, a conservative linear system can possess several independent linear invariants, i. e. there exist K linear independent vectors n i with n T i A = 0 and hence, n T i y remains constant for all times.To avoid the issue with asymptotic stability, we can weaken the requirement and demand only stability instead, while at the same time, we can require that the numerical approximations tend to the unique steady state solution y * of y ′ = Ay, y(0) = y 0 , which always satisfies n T i y * = n T i y 0 .In summary, we are looking for MPRK schemes for which a stable steady state y * of a positive and conservative linear system y ′ = Ay, becomes a stable fixed point of the MPRK scheme for all time step sizes ∆t, and in addition, the iterates y n tend to y * for all initial values y 0 that satisfy n T i y 0 = n T i y * .The fact that an MPRK scheme y n+1 = g(y n ) is not a general linear method makes it complicated to find MPRK schemes with the desired properties, since the application of an MPRK scheme to a linear system results in a nonlinear iteration of the form y n+1 = R(∆tA, y n )y n .In addition, the conservation property implies the existence of infinitely many nonhyperbolic fixed points y * ̸ = 0 of the map g, i. e. the Jacobian Dg(y * ) has eigenvalues with absolute value equal to one.Hence, a linear stability theory for MPRK schemes must be a stability theory for non-hyperbolic fixed points of nonlinear iterations.One such approach to study stability is based on the center manifold theory of dynamical systems and was introduced in [19,20].Assuming that the initial value y 0 is sufficiently close to the fixed point y * , the theory provides sufficient conditions for the stability of y * , as well as for the convergence of the iterates to the steady state of the corresponding initial value problem.Thus, this stability theory would almost suffice to find the desired schemes if the theory did not make local statements only, with respect to the initial value.However, to the authors' knowledge this is the only approach to study stability of MPRK schemes so far.
The stability theory of [19,20] was used to investigate the linear stability of MPRK22 schemes and it was proven therein that MPRK22(α) schemes with α ≥ 0.5 are linearly stable and that their iterates tend to the correct fixed point y * , whenever the initial value y 0 is sufficiently close to y * .The stability theory was used in [21] to find SSP-MPRK schemes with the same properties.A stability analysis of the third order MPRK schemes from [6] and the MPDeC schemes of [9] was carried out in [22].We also want to note, that the stability theory of [19,20] is not only applicable to MPRK schemes and was also used in [23] to analyze the linear stability of BBKS [10,12,13] and GeCo schemes [14].Furthermore, it was used in [24] to investigate the stability of an MPRK scheme in the context of a nonlinear PDS.
Even though the stability results of [19,20] are only valid in a sufficiently small neighborhood of the fixed point, the numerical results in [19][20][21] suggested that the local stability might actually be a global one, just like it is the case for A-stable general linear methods.
However, in [22,25] it was discovered that the MPDeC(8) method with equidistant nodes is indeed only locally stable.To check the conjecture that this is due to the negative RK parameters of the MPDeC(8) scheme, the family of MPRK22(α) schemes with negative RK parameters, i. e. α < 0.5, was investigated in [26].The numerical experiments presented therein show that MPRK22(−0.5) is another MPRK scheme, which is only locally stable.
The aim of this paper is to summarize and extend the results of [26] and to justify the conjecture that this local linear stability behavior, which is unknown from general linear methods, only occurs for MPRK schemes if the Butcher tableau of the underlying RK scheme contains negative values.

Linear stability of MPRK22(α) schemes
This section summarizes some of the results of [26].
The idea of MPRK schemes is to modify explicit RK schemes by introducing additional weighting factors that ensure unconditional positivity and conservation.Destruction terms are multiplied by weights with respect to the equation they appear in.Production terms are multiplied by the same weights as their corresponding destruction counter parts.If all parameters of the underlying RK schemes are nonnegative, there is no difference between production or destruction terms on the continuous and discrete level.But if production or destruction terms are multiplied by a negative RK parameter they switch their roles from the continuous to the discrete level, which has to be dealt with appropriately within the implementation of the scheme.
The MPRK22(α) schemes were introduced in [5] and are based on general second order explicit RK schemes, i. e.
In [5] only nonnegative RK parameters were considered, which is the case, if α ≥ 1 2 .For 0 < α < 1 2 the parameter b 1 becomes negative and for α < 0 the parameters a 21 and b 1 are negative.Consequently, we need to distinguish these three different cases, for varying values of α.
In a form which is suitable for positive as well as negative RK parameters, the MPRK22(α) schemes for the solution of (1) can be reformulated as p ij (y (1) ) for i = 1, . . ., N , with i ) 1/a21 , i = 1, . . ., N and the index function The purpose of the index function is to decide, whether a term is a production or destruction term and to choose the weighting factors accordingly.We also note that the index function (6) was introduced in [9] for the definition of MPDeC schemes.It was proven in [5] that MPRK22(α) schemes with α ≥ 1 2 are unconditionally positive and conservative second order schemes.The same is true for 0 < α < 1  2 and α < 0, which was proven in [26] by a straight-forward modification of the proof given in [5].For a general framework concerning the order conditions of MPRK schemes we refer to [17].
To examine the linear stability of MPRK22(α) schemes in terms of the stability theory [19,20], we consider their application to positive and conservative linear PDS (3).In terms of a shorter notation we define denotes the diagonal matrix with the elements of the vector v = (v 1 , . . ., v n ) T ∈ R N on the diagonal.As mentioned above, we need to distinguish between the following three cases.

Case α ≥ 0.5
This case was already considered in [19,20] and we present the results for the sake of completeness.
As α = a 21 > 0 implies that all RK parameters are nonnegative, the application of an MPRK22(α) scheme ( 5) to (3) yields In comparison to the underlying RK scheme a diagonal matrix with Patankar-weights was introduced on the left of diag(A) in the destruction parts and on the right of B in the production parts.This is done in the other cases as well.
The dependence of y n+1 on y n can be expressed by an implicit function g, i. e. y n+1 = g(y n ).Each steady state y * of ( 3) is also a fixed point of g and has the eigenvalues R(∆tλ), where λ is an eigenvalue of A and the stability function is As a result, |R(z)| < 1 for all z ∈ C − , which implies stability of fixed points and convergence towards the steady state of the underlying initial value problem.Hence, the MPRK22(α) schemes with α ≥ 0.5 fulfill all desired properties, apart from the fact that stability could only be proven in a sufficiently small neighborhood of the fixed point.Furthermore, we want to emphasize that in this case, apart from ∆t, stability depends only on λ like in the continuous case.
2.2 Case 0 < α < 0.5 A similar analysis as for α ≥ 0.5 can be conducted for the situation that 0 < α < 0.5.This was done in [26] and the results will be summarized here.In this case we have b 1 < 0 and a 21 , b 2 > 0. Consequently, the scheme (5) applied to (3) reads The stability behavior is more complicated than for α ≥ 0.5 as the Jacobian Dg(y * ) in general also depends on A T and the steady state y * itself, see [26] for details.However, if we restrict ourselves to the specific system (4), then A technical computation shows that for every α there exists a z * < 0 for which |R(z * )| > 1.Hence, the MPRK22(α) schemes are only conditionally stable for 0 < α < 0.5.To be precise, z * is given by

Case α < 0
In this case we have a 21 , b 2 < 0 and b 1 > 0. Hence, terms multiplied by a 21 or b 2 change their role from the continuous to the discrete level.As a result, MPRK22(α) schemes (5) applied to (3) become + b 1 diag(y n+1 /σ) diag(A)y n + b 2 diag(A) diag(y n+1 /σ)y (2) . (7b) Again we summarize the results of [26] here.But as this is the important case for the purpose of this paper, we also give proofs in Appendix A, at least for specific cases.As for 0 < α < 0.5 the Jacobian Dg(y * ) in general depends on A T and the steady state y * itself.But for the purpose of this paper it is sufficient to restrict ourselves to system (4), in which case it can be seen that the dependence on A T and y * within (A4) disappears due to the property (A5) and we have as well as The nonzero eigenvalue of A in ( 4) is λ = −(a + b) < 0. Hence, we only need to consider z ∈ R with z < 0. For every α < 0 the function R is monotonically increasing on (−∞, 0) with < R(z) < 1 for z < 0 and α < 0.
If we restrict α additionally to α ≤ −0.5, then |R(z)| < 1 for z < 0 and it follows that MPRK22(α) schemes with α ≤ −0.5 are unconditionally stable, i. e. the stability is independent of ∆t.Moreover, convergence to the steady state of the corresponding initial value problem is guaranteed as well, whenever the initial values is sufficiently close to the fixed point.The situation is different for −0.5 < α < 0, for which the schemes are only conditionally stable.For −0.5 < α < 0 stability requires z * < z < 0 with As mentioned before, see [26] or Appendix A for details.Now we know that MPRK22(α) schemes with α ≤ −0.5 are unconditionally stable and converge to the steady state of the corresponding initial value problem, at least when applied to the linear system (4) and the initial value is not to far away from the steady state.As mentioned above, the numerical results in [19][20][21] suggested that this local behavior might actually be the global behavior.But it was observed in [26] and [22,25] that MPRK22(-0.5)and MPDeC (8) with equidistant nodes are indeed only locally stable.Such a linear stability behavior is unknown from general linear methods.
The conjecture that this can only occur in the presence of negative RK parameters is supported in Section 3, where we also show that there are further values of α with α < −0.5 for which the same stability behavior can be observed.Moreover, the experiments with α ≥ 0.5, corresponding to non-negative RK parameters, do not show a case for which the stability properties are local.

Numerical results
In this section we demonstrate numerically that for some MPRK22(α) schemes with α < 0 the linear stability and convergence to the correct steady state is indeed only given locally, i. e. the initial value must be sufficiently close to the fixed point.To show this, it is already sufficient to consider the linear system with initial value where the restriction on δ is necessary to ensure positive initial data.The solution of this initial value problem is and its steady state is The following computations were performed with MATLAB 2023b.The MATLAB code is available from https://github.com/SKopecz/locstabMPRK22.git.
First, we consider system (8) with a = 20 for which the system is nonstiff.In Figure 1 we see numerical solutions computed by MPRK22(−0.5)for two slightly different values of δ.The rather large time step size ∆t = 1, which is of course unsuitable for accuracy, was chosen to demonstrate the stability behavior of the scheme.As described in Section 2.3 the MPRK22(−0.5)scheme is stable and the iterates tend to the steady state of the corresponding initial value problem for all time step sizes, as long as the initial value is close enough to the steady state.This behavior is depicted in Figure 1a.Now, in Figure 1b the value of δ, i. e. the distance to the steady state, is slightly increased and as a result the numerical approximations tend to a spurious steady state.As discussed above, we would like to have a linear stability behavior similar to A-stability, i. e. the stability should not depend on the initial values.But here it is demonstrated that MPRK schemes exist for which stability crucially depends on the initial value.
To gain a better insight into the dependence of the stability on δ, we compute the steady states for N = 200 equidistantly spaced samples of δ between 0 and 0.5.In each case the steady state is computed by performing M = 10 4 steps of the MPRK22(−0.5)scheme to obtain y M .To check if the correct or a spurious steady state is approached we compute the distance d(α, δ) = ∥y * − y M ∥ ∞ .
We expect d(α, δ) ≪ 1 if the scheme is stable and tends to the steady state of the corresponding initial value problem and d(α, δ) ≫ 0 if a spurious steady state is reached.Figure 2a shows a plot of d(−0.5, δ) for 0 ≤ δ < 0.5.We see that the stability behavior changes abruptly from stable to unstable for some δ between 0.23 and 0.24 as suggested by Figure 1.If we increase the stiffness of the system by choosing a = 200 the critical value of δ, for which the change from stable to unstable behavior occurs, decreases, as can be observed from Figure 2b.In this case the stability behavior changes from stable to unstable for δ ≈ 0.06.Next, we want to answer the question if α = −0.5 is the only value of α for which the linear stability of MPRK22(α) crucially depends on the initial value.To answer this question, we again compute the distance d(α, δ), but this time we vary both δ and α.We use an equidistant grid with 160 samples of δ and 241 samples of α for 0 < δ < 0.5 and 0 < |α| ≤ 2. Again, we use a = 200 in (8), ∆t = 1 and M = 10 4 steps to compute the steady state.The result of this computation can be seen in Figure 3a and we note that the MPRK22(α) scheme is not defined for α = 0, which is not indicated in the plot.Based on Figure 3a, we can make the following statements.First, the expected unstable behavior for 0 < |α| < 0.5, which was discussed in Sections 2.2 and 2.3, is clearly visible.Second, for α ≥ 0.5 we observe stable behavior independent of δ.Therefore, local stability could actually be global stability for α ≥ 0.5.Third, for α ≤ −0.5 there exists an α * < −0.5 such that stability of MPRK22(α) is indeed only a local stability and for α < α * no unstable behavior can be observed.To see this more clearly, Figure 3b shows a zoom of Figure 3a with higher resolution for −0.6 ≤ α ≤ −0.5.As we can see, α = −0.5 is not the only value of α for which the stability depends on the initial value.The same is true for all α between α * ≈ −0.56 and −0.5.

Summary and outlook
We have discussed the linear stability behavior of MPRK22(α) schemes for all α ̸ = 0 based on the local stability theory of [19,20].Moreover, we have confirmed that negative RK parameters in MPRK schemes can lead to linear stability actually being a local property.In particular, we have shown that there are MPRK22(α) schemes with α < −0.5 that are linearly stable but converge to a spurious fixed point instead of the steady state of the corresponding linear system if the initial value is not sufficiently close to the steady state.Such a linear stability behavior is unknown from general linear methods.We further note that we cannot recommend to use MPRK22(α) schemes with α < 0.5 since apart from the time step size restrictions for 0 < α < 0.5, the stability of these numerical methods shows an artificial dependence on A T and y * that is not present in the continuous case.In addition, the convergence towards the steady state of the corresponding initial value problem cannot be guaranteed in the case of general linear systems by means of the stability theory from [19,20], see [22].Nevertheless, for MPRK22(α) schemes relevant in practice, i. e. α ≈ 1, the parameters of the underlying RK scheme are nonnegative and no restrictions on the initial value could be found numerically.Therefore, for these methods, it still seems that the local linear stability is actually a global linear stability.
Since we would like to have MPRK schemes with properties similar to A-stable general linear methods, it is therefore essential that the convergence to the steady state of the corresponding initial value problem is given for all initial values which satisfy n T i y 0 = n T i y * , where the vectors n i define the independent linear invariants of the system matrix.For this reason, there is a great need for a theory that makes statements not only about local stability but also about global stability.
By differentiation of (A3a) we obtain the Jacobians By differentiation of (A3b) one can write which shows that Dg(y * ) as given in (A2) in general depends on A T and y * .However, hereafter we only consider the linear system (4) for which y * = s(b, a) T for some s ∈ For z < 0 and α < 0 every term in the nominator is positive and the denominator is positive as well.Hence, R ′ (z) > 0 for all z < 0 and α < 0, which implies that R is monotonically increasing on (−∞, 0).Furthermore, we find R(0) = 1 and .
Together with the monotonicity of R this results in < R(z) < 1 for z < 0 and α < 0. According to the stability theory of [19,20] it follows that MPRK22(α) schemes with α ≤ −0.5 are unconditionally stable and converge to the steady state of the corresponding initial value problem independent of the time step size ∆t, if the initial value is close enough to the fixed point.On the other hand, for −0.Hence, stability for −0.5 < α < 0 requires z * < z < 0.