Relative error long-time behavior in matrix exponential approximations for numerical integration: the stiff situation

In the stiff situation, we consider the long-time behavior of the relative error γn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _n$$\end{document} in the numerical integration of a linear ordinary differential equation y′(t)=Ay(t),t≥0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y^{\prime }(t)=Ay(t),\quad t\ge 0$$\end{document}, where A is a normal matrix. The numerical solution is obtained by using at any step an approximation of the matrix exponential, e.g. a polynomial or a rational approximation. We study the long-time behavior of γn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _n$$\end{document} by comparing it to the relative error γnlong\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _n^{\mathrm{long}}$$\end{document} in the numerical integration of the long-time solution, i.e. the projection of the solution on the eigenspace of the rightmost eigenvalues. The error γnlong\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma _n^{\mathrm{long}}$$\end{document} grows linearly in time, it is small and it remains small in the long-time. We give a condition under which γn≈γnlong\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _n\approx \gamma _n^{\mathrm{long}}$$\end{document}, i.e. γnγnlong≈1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\gamma _n}{\gamma _n^{\mathrm{long}}}\approx 1$$\end{document}, in the long-time. When this condition does not hold, the ratio γnγnlong\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\gamma _n}{\gamma _n^{\mathrm{long}}}$$\end{document} is large for all time. These results describe the long-time behavior of the relative error γn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _n$$\end{document} in the stiff situation.


Introduction
Consider the ordinary differential equation (ODE) where A ∈ ℝ d×d and y(t) ∈ ℝ d , and consider, over the mesh of constant stepsize h > 0 , a numerical solution of (1.1) given by where R ∶ D ⊆ ℂ → ℂ is a analytic approximant of the exponential e z , z ∈ ℂ . When the numerical solution is obtained by a Runge-Kutta (RK) method, the approximant R is the stability function of the RK method and it is a polynomial or a rational function.
The paper [9] analyzed in the non-stiff situation the time behavior of the normwise relative error in case of a normal matrix A. It seems to be the first paper in literature dealing in detail with the relative error time behavior of numerical solutions of ODEs. This is quite surprising because relative errors are generally considered better than absolute errors as quality measures of approximations. Indeed, componentwise relative errors are involved in the stepsize control mechanism (see [12]).
The present paper continues to analyze, in case of A normal, the error n by considering its long-time behavior in the stiff situation. Next subsection, with all its subsubsections, contains the basic material for facing such an analysis. Part of this material was introduced in [9].

Small and large
We set that, for a ≥ 0 , "a is small" is the same as " a ≪ 1 " and "a is large" is the same as " a ≫ 1".
For b ≥ 0 and c > 0 , b ≪ c means b c ≪ 1.

The meaning of "it is expected"
In the paper, we often say "it is expected S", where S is a statement, with the meaning that the statement not S is "unlikely" or "unusual" or "extreme". Sentences of this form can seem vague, although they are able to convey significant information. However, they are never used in definitions or theorems, which are stated in a precise manner without any such type of vagueness. The sentences are used for a better understanding of technical notions and results.
By introducing probability measures on data, we could made "it is expected S" mathematically precise, but this is out of the scope of the present paper.

The spectrum of A
The spectrum of the normal matrix A, where 1 , 2 , … , p are the distinct eigenvalues of A, is partitioned by decreasing real part in the subsets 1 , 2 , … , q (see Fig. 1 with 0 = i 0 < i 1 < ⋯ < i q = p , and For i = 1, … , p , let P i be the orthogonal projection on the eigenspace of i . For j = 1, … , q , let For a nonempty subset of , let

The initial value y 0
We assume y 0 ≠ 0 . Thus y(t) = e tA y 0 ≠ 0 , for any t, and the relative error n is defined for any n. Let be the normalized initial value. Let The generic situation for the initial value y 0 is * = . In order to use simpler notations, we assume this generic situation. If it does not hold, then below we have to see 1 , … , q as * 1 , … , * q without the sets * j that are empty. In other words, we see 1 as j * 1 where 2 as j * 2 where j = { i j−1 +1 , i j−1 +2 , … , i j } Re i j−1 +1 = Re i j−1 +2 = ⋯ = Re i j = r j j = 1, 2, … , q, Q j ∶= ∑ and so on. Of course, when we do this, the number q of sets in 1 , … , q is no longer equal to the number of possible real parts in the spectrum , but it is equal to the number of possible real parts in * .

Rightmost and non-rightmost eigenvalues
The set 1 is the set of the rightmost eigenvalues. The set is the set of the non-rightmost eigenvalues. We assume q > 1 in order to have − ≠ �. By recalling the definitions (1.4), we set

The numbers ˇj
For j = 2, … , q , i.e. for any non-rightmost real part, let Observe that It is expected | 2 | non-small.

Dimensionless quantities
We use the dimensionless stepsize h 1 , or h , and the dimensionless time t 1 , or t , rather than the stepsize h and the time t, respectively, because they are small or large independently of the unit used for time.
In this paper, when we say that a certain quantity is small or large, this quantity is always dimensionless.
The numbers j defined above are dimensionless, as well as the errors i now introduced.

The errors i
We assume that the approximant R has order l, where l is a positive integer. This means with C ≠ 0 . It is assumed that the domain D of R includes a neighborhood of zero. Moreover, we assume h i ∈ D , i = 1, … , p.
We introduce the complex numbers where is the relative approximant. The numbers i are logarithmic errors of R as an approximant of the exponential, since For a nonempty subset of , we have as h → 0 , where and are defined in (1.4).

Local relative errors and global relative errors
As particular cases of (1.9), we obtain . This result is useful in the base situation, where it is expected E 1 ≪ 1.

Remark 1.2
1. By taking a small c in the previous theorem, we have This holds for times t n 1 ≤ x , where x > 0 is such that xE 1 ≪ 1 . If 1 is constituted by a real eigenvalue or by a complex conjugate pair of eigenvalues (the generic situation for the matrix A), we have K 1 = 1 and then = 0 for any n. For the case 0 ∈ 1 and 1 ≠ {0} , see point 5 of Remark 5.2 in [9].

Long-time behavior of n
We want to study the long-time behavior of the error n . This is done by comparing it to the error long n . Since in the long-time the solution y becomes the solution y long whose error n is just long n , it is quite reasonable to have n ≈ long n in the long-time. Indeed, at point 4 of Remark 5.3 in [9], it is stated the following result.
For any > 0 , there exist H 0 > 0 (independent of ) and s ≥ 0 (dependent on ) such that, for h ≤ H 0 and s ≤ t n ≤ c E , we have n ≈ long n with degree .

Remark 1.3
The theorem assumes q > 1 . If q = 1 , then n = long n for any n. In addition, it assumes 0 ∉ 1 . If q > 1 and 1 = {0} , then long n = 0 for any n and it does not make sense look at n ≈ long n , since this implies n = 0 . For the case q > 1 , 0 ∈ 1 and 1 ≠ {0} , see point 6 of Remark 5.3 in [9].
The previous theorem is of interest in the non-stiff situation, where the condition h ≤ H 0 is not restrictive. In fact, in the non-stiff situation it is expected h non-large.
On the other hand, the result is not useful in the stiff situation, since the condition h ≤ H 0 is restrictive. In fact, in the stiff situation it is expected h non-small.

The contents of this paper
The present paper wants to study the long-time behavior of the relative error n in the stiff situation. As above, this is done by comparing it to long n . In the stiff situation, it is important to have n ≈ long n in the long-time. In fact, if this happens, since long n is small up to large times t n 1 , we have the very surprising fact that the error n is small in the long-time, although the stepsize h is 23 Page 12 of 62 tuned only for having a small local relative error of the long-time solution and, because of this, the local relative error of the solution is not small.
In other words, when we are interested in the numerical integration of the solution in the long-time, we can start from the beginning with a stepsize suitable for integrating with a small local relative error the long-time solution, larger than the stepsize suitable for integrating with a small local relative error the solution, and in the long-time we will have a small error n .
As in [9], we confine our attention to normal matrices. This should be not considered as a limitation, since the class of the normal matrices is sufficiently large to include important types of matrices and, moreover, the test problem (1.1) with A normal shows unexplored and interesting situations in numerical ODEs.
The plan of the paper is as follows.

Replies to general questions or criticisms
This final subsection includes replies to general questions or criticisms which could be issued about the contents of this paper.
-Question. What is the motivation of this paper?
Reply. This paper studies the relative error of numerical approximations of ODEs, although confined to linear systems with normal matrix. Of course, the absolute error and the relative error of the numerical approximations have the same order of convergence with respect to the stepsize h, but they have a different time behavior in the numerical integration of a solution spanning over various orders of magnitude.
The motivation for studying the relative error time behavior in numerical ODEs, as the present paper is doing, comes from the following two facts: -as it is widely recognized, the relative error is an important measure for the quality of an approximation, often better than the absolute error; -there has been no attention in the numerical ODEs community about the relative error time behavior of numerical approximations.
Anyway, the fact that in the numerical ODEs field the relative error is considered important is attested by the numerical solvers, which accept as an input argument a tolerance on the componentwise relative error. Thus, this paper (similarly to [9] and [10]) try to fill this gap between theory, where there are not studies on the relative error, and practice, where the relative error is used.
-Question. What is the relevance of the results achieved?
Reply. For the numerical ODEs community, it should be of interest to know the relative error time behavior of numerical approximations of the ODE (1.1) with A normal. The results achieved describes this time behavior and their relevance is that they give a new prospective on the numerical integration errors. We can summarize this new perspective in the following points.
-In the non-stiff situation, the relative error is small and it grows linearly in time. Moreover, this linear growth is determined in the long-time only by the rightmost eigenvalues. -In the stiff situation, the relative error is not small at the beginning of the numerical integration and it is not guaranteed that in the long-time it will become small, with a linear growth determined only by the rightmost eigenvalues. This happens if and only if a certain condition is satisfied and this condition is a novelty in the numerical ODE theory. -Gauss RK methods, despite they are considered stable in the classic numerical linear stability theory (they are A-stable methods), are not suitable to have the above condition satisfied. On the other hand, Radau and Lobatto IIIC RK methods are suitable to have this condition satisfied.
-Criticism. Componentwise relative errors where y n,i and y i (t n ) , i = 1, … , d , are the components of y n and y(t n ) , should be considered (as in the numerical ODE solvers), not the normwise relative error (1.3).
Reply. In literature both normwise relative errors and componentwise relative errors are considered as quality measures of vector approximations (see [2]). The componentwise approach has the advantage that it gives information on the precision of the components, but it has the drawback that the components must be nonzero (when some component becomes zero, we need to switch to the absolute error). On the other hand, the normwise approach can give anyway information about the componentwise relative errors (for example, a large normwise relative error implies that some component has a large relative error) and it works also when some component becomes zero. Reply. In mathematical modeling and numerical analysis there is a threshold in the order of magnitude of quantities (scalars or vectors) under which they are considered zero. Under the threshold, it is important to use the absolute error for approximations, since they are considered approximations of zero. But, in case of a solution of (1.1) which is going to zero, and so it is spanning over several orders of magnitude, it could be of interest to compute with a good precision this solution for the orders of magnitude larger than the threshold. In this situation, the relative error is important.
Of course, the numerical analyst's point of view is that the threshold is the order of magnitude of the machine epsilon, but in applications this threshold can be larger.
As an example, we can consider the radioactivity decay of radionuclides, where the activity a(t) (measured in becquerel (Bq) by a Geiger counter) of a given amount of radionuclide satisfies a � (t) = − a(t) with > 0 . For a decay chain, we have a � (t) = Aa(t) , where A is a lower bi-diagonal matrix, the so-called Bateman equation. The threshold could be the order of magnitude 10 2 Bq∕kg of the background radiation. Of course, this threshold becomes a much smaller ten power by using an unit larger than the becquerel, e.g. the curie. It could be interesting to numerically compute with a good precision a solution a(t) whose initial value has order of magnitude 10 6 Bq∕kg (like in a nuclear plant accident). Since the solution becomes small compared to the initial value, using the relative error for the approximations of the solution is better than using the absolute error when the solution is not yet considered as zero.
Another example could be a space discretization of the heat equation, with homogeneous Dirichlet boundary condition, by the method of lines. In this case, the space discrete temperature approaches zero (the border temperature) and under a given threshold in the order of magnitude, say 10 −2 • C , it can be considered zero. But, over this order of magnitude, the temperature is not zero and it becomes important to use the relative error for time-space approximations, especially when the solution spans over several orders of magnitude due to an initial value with order of magnitude larger than the threshold, for example 10 2 • C.
We remark that the analysis in this paper also consider the situation where the solution, instead of approaching to zero, grows up to large values with respect to the initial value. Also in this situation the relative error is important. -Criticism. The paper considers ODEs (1.1) with matrix A normal. Such problems can be diagonalized with a unitary transformation and then one can assume without loss of generality that A is diagonal. Reply. In the paper, we do not assume from the beginning that A is diagonal because this does not simplify the exposition. In fact, the analysis presented starts from the fundamental relation (4.6) given below for the relative error, which maintains the same form when A is diagonal. We have such a net expression for the relative error precisely for the possibility to reduce to the diagonal case by a unitary transformation. Hence, the assumption that A is diagonal is already implicitly done when one decides to deal with a normal matrix. -Criticism. Since it is possible to reduce to the diagonal case, it would be sufficient to study the behavior of the numerical scheme at a scalar problem, which is really trivial.
Reply. Although we can reduce to a linear systems of uncoupled scalar differential equations, this does not mean that they are fully uncoupled in the numerical scheme, since we are using the same stepsize h in all scalar equations. This reflects the fact that the numerical scheme is applied to an ODE (1.1) with a matrix A in general non-diagonal, without thinking to diagonalize it in advance. Moreover, the analysis of the present paper requires to have rightmost and nonrightmost eigenvalues. In other words, we need eigenvalues with different real parts, i.e. an ODE (1.1) with different time scales. The case of a sole scalar equation is not considered. Anyway, we can observe that in the base situation for a scalar equation, the relative error n = long n is expected to be small and linearly growing in time up to large times.

Examples
In this section, we give two examples of stiff situations where the error n is not small from the beginning of the numerical integration and it grows without to approach in the long-time to the small error long n . We remind that the stability region of the approximant R (see [5]) is the set and the order star of R (see [5][6][7]13]) is the set where S is the relative approximant given in (1.8). The complementary set of S is

Same approximant with different ODEs
As first example, we consider the ODE (1.1) with the symmetric matrix whose eigenvalues are a and b with relevant eigenvectors (1, 1) and (1, −1) , respectively. We consider a = −1 and the following three possibilities for b: The initial value is y 0 = (2, −1) , for which we have The solution y quickly approaches to the long-time solution y long : we have . For the numerical integration, we use the Taylor approximant of order five with stepsize h = 1 5 over N = 25 steps up to t N = Nh = 5. We have:   Starting from a non-small 1 (remind that 0 = 0 ), the error n goes down to the small error long n . In the long-time, we have small errors n although the stepsize is tuned only for having a small 1 , without any concern about 2 .
For the possibility P2), we see in Fig. 3 the same as in Fig. 2. As in P1), starting from a non-small 1 , the error n goes down to long n , although long n is reached at a larger time with respect to P1).
Finally, for the possibility P3), we see in Fig. 4 the same as in Figs. 2 and 3. Unlike P1) and P2), the error n does not go down to long n , but it continues to grow.

Order star and stability region
Fixed a = −1 , we are interested in understanding for which b, with b < a , we have n ≈ long n in the long-time. This happens in P1) and P2), but not in P3). Order star and stability region for the Taylor approximant of order five are depicted in Fig. 5.
The values of |S(hb)| and |R(hb)| are: Observe that hb ∈ R for all three possibilities and hb ∈ S c only in P1). In other words, by looking at the negative real axis of Fig 5, hb lies in the red region for all three possibilities and hb lies in the white finger only in P1).
In the next Sect. 4, we will see a condition on hb for having n ≈ long n in the long-time. When it does not hold, we have n long n ≫ 1 for all time. The condition is something between hb ∈ S c (i.e. to stay in the white finger) and hb ∈ R (i.e. to stay in the red region). Indeed, to have hb ∈ S c is sufficient, but not necessary, for this condition on hb and to have hb ∈ R is necessary, but not sufficient.

Same ODE with different approximants
As second example, we consider the ODE (1.1) with the normal matrix with and Q with orthonormal columns u 1 , Consider the initial value y 0 = (3, 3, 3, −2) for which we have The solution y consists of two decaying oscillations: the fast oscillation y − y long decays faster than the slow oscillation y long and, in the long-time, only the slow oscillation is present. We have y(t) ≈ y long (t) if (look at (1.17)).
Assume that the numerical integration of the ODE is accomplished by the fourth order two-stage Gauss RK method, corresponding to the (2, 2)−Padé approximant and by the third order two-stage Radau RK method, corresponding to the (1, 2)− Padé approximant Both methods are applied with stepsize h = 1 10 over N = 100 steps up to t N = Nh = 10 . Observe that such a stepsize is not suitable for approximating the fast oscillation.
We are in the stiff situation: Since In the upper part of Fig. 6, we see the trajectory t n ↦ y 1 t n , y 2 t n in the plane ℝ 2 for the first two components of the exact solution y(t n ) , when t n ∈ [8,10] . In the middle and lower parts, we see the trajectory t n ↦ y n,1 , y n,2 for the first two components of the numerical solution y n , when t n ∈ [8,10] . Middle part for the Gauss RK method and lower part for the Radau RK method.
For the long-time t n ∈ [8,10] , where only the slow oscillation y long is present, the exact components y 1 (t n ) and y 2 (t n ) are equal and have order of magnitude 10 −4 . The Gauss RK method exhibits numerical components y n,1 and y n,2 of order of magnitude 10 0 . On the other hand, the Radau RK method exhibits accurate numerical components y n,1 and y n,2 , although the stepsize is not suitable for approximating the fast oscillation.
In Fig. 7, we see the error n , for n = 0, 1, … , N , for both approximants: for the Gauss RK method the error continues to grow and for the Radau RK method it goes down to long n .

Order star and stability region
Fixed 1 = −1 + i and 3 = −3 + 1000i , we are interested in understanding for which approximants we have n ≈ long n in the long-time. In our situation, this happens for the Radau RK method, but not for the Gauss RK method.
Order star and stability region for such approximants are shown in Fig. 8. We have h 3 ∈ R for both methods, since they are A-stable. On the other hand, we have h 3 ∈ S c only for the Radau RK method: In Sect. 4, we will see a condition on the approximant for having n ≈ 3 The appropriate definition of n ≈ long n in the long-time In the following, we assume to be in the base situation. Then, it is expected E 1 small and h 1 non-large. To make easier the exposition, we assume E 1 small and h 1 non-large. Since E 1 is small, the error long n grows linearly in time and it is small up to large times t n 1 .
We also fix a number c > 0 and let (The number c plays a role similar the number c appearing in Theorems 1.1, 1.2 and 1.3). As a reference value for c, one can take c = 1 . As a matter of generality, we do not confine c only to this value. In all theorems below, it is stated for which c > 0 3494 for the Gauss RK method 0.0270 for the Radau RK method.
they are valid. However, when the theorems are applied, c is considered non-small, so we have and such that g(c) < 1 , i.e. c < 1.2564 , with 1 − g(c) non-small. In order to describe the long-time behavior of the error n , we compare it to long n and we are interested in whether or not n ≈ long n in the long-time. Here, "in the long-time" does not mean t n 1 → +∞ . In fact, it is not of great interest to consider what happens for t n 1 → +∞ , since long n becomes non-small for a sufficiently large t n 1 . It is of interest to have n ≈ long n starting from times t n 1 such that long n is still small. So, we introduce the following definition. In the definition, we consider times t n 1 up to . Observe that if K 1 is not large (remind (1.12) and remind that K 1 = 1 is the generic situation for the matrix A), then the error long n is not small for t n 1 at the end of the interval [0, ]. In fact, for t n 1 ∈ [ , ] , where ∈ (0, 1] is not small, by Theorem 1.2 we have and the right-hand side in this inequality is not small.

The definition of n ≈ long n in the long-time with monitor function
We can make the previous definition more precise by a monitor function.
We say that n ≈ long n in the long-time with monitor function s if, for any > 0 , we have where e n is such that n = long n (1 + e n ).

Remark 3.1
In the previous definition, we also allow monitor functions In this case, we have to specify that (3.1) holds for ∈ (0, a] and (3.2) holds for ∈ (0, a] and ≥ b.  In summary: Regarding the satisfiability of s( , ) ≪ 1 , observe that s satisfies (3.1) and we have ≫ 1.

Analysis of the long-time behavior of n
In the paper [9], it was presented an analysis of the long-time behavior of the error n important for the non-stiff situation. In the present paper, it is developed another type of analysis important for the stiff situation. In this new analysis, the complex numbers w i and i introduced below are important.

The numbers w i
For any i ∈ − , i.e. for any non-rightmost eigenvalue, we introduce the complex number where j = 2, … , q is such that i ∈ j . We set

The numbers ˛i
For any i ∈ − , we introduce the complex number If n ≈ long n in the long-time with monitor function s and for some small > 0, then n ≈ long n in the long-time. In particular, We set Since we have As a consequence we obtain Here are some observations about .
-It is expected | | non-small. In fact, let i ∈ j , with j = 2, … , q , be a non-rightmost eigenvalue such that The case where | | is small, i.e.
with |e| ≪ 1 , is "unlikely". Observe that it is expected | j | non-small. -In the non-stiff situation, it is expected negative non-small. In fact, it is expected | 2 | non-small and, in the non-stiff situation, it is expected that the righthand side of (4.4) is small and then it is expected | − 2 | small.

The basic theorem
The next theorem is, in our new analysis, the analog of Theorem 5.3 in [9] (which was suitable for studying the long-time behavior of n in the non-stiff situation).
In the following, we continue to assume 0 ∉ 1 , but all our conclusions are valid (with easy adaptations) also for the case 0 ∈ 1 and 1 ≠ {0}.

A first result
We give a first theorem about n ≈ long n in the long-time with a monitor function. (4.10) where the last equality follows by (4.3) (remind (1.13)). So, for t n 1 ∈ [s, ] , we obtain By inverting the function f, we obtain the monitor function (4.10). ◻   Proof Use the previous theorem with = 1 = E 1 ≪ 1 and observe that (4.12) ◻ It is expected that the last three conditions in (4.13) are satisfied. In fact, since E 1 ≪ 1 , they are not satisfied only in "extreme" cases. Moreover, in the non-stiff situation, it is expected that the first condition is satisfied. In fact, since in the non-stiff situation it is expected max i ∈ − | i| h 1 ≪ 1 , the first condition is not satisfied only in "extreme" cases.
So, we can state the following important conclusion.

Conclusion 4.5
Suppose to be in the non-stiff situation. It is expected n ≈ long n in the long-time.

The condition A
We introduce the condition where W and are defined in (4.1) and (4.2), respectively.
Next theorem shows that, under the condition A, we have n ≈ long n in the longtime with a new monitor function different from (4.10).
If A holds, then n ≈ long n in the long-time with monitor function defined for x ≥ 1 and for > 0 such that the right-hand side of (4.14) with x = 1 is greater than or equal to 1, so we have s( , x) ≥ 1 for x ≥ 1 and such .
the advantage that it is no longer necessary to suppose ≤ a (where is a given at point 1) above) and ≥ 1 . However, we prefer to use the old monitor function (4.14) because it has an explicit expression.
The previous theorem with c = 1 gives the next important results.   Proof Use the previous theorem with It is expected that if A holds, then the three conditions in (4.16) are satisfied. In fact, since E 1 ≪ 1 , they are not satisfied only in "extreme" cases. So, it is expected that if A holds then n ≈ long n in the long-time. Of course, we already know that in the non-stiff situation it is expected n ≈ long n in the long-time, independently of the condition A, as well as we know that it expected that A holds in the non-stiff situation.
So, what is really important is the following conclusion.

Conclusion 4.9
Suppose to be in the stiff situation. It is expected that if A holds, then n ≈ long n in the long-time.

Order star and stability region
The condition A can be related to the order star and the stability region of the approximant (recall the beginning of Sect. 2).

The region R x
For any x ∈ ℝ , let

R.
For i ∈ − , we have Thus, the condition A can be restated as In the case r 1 = 0 , it becomes according to Theorem 4.11 about the cases r 1 ≤ 0 and r 1 ≥ 0.

The conditions B and C
When A does not hold, we have B or C, where and (4.17) In the next section we study, what happens when A does not holds. Of course, when A does not hold, it is expected that B holds.

When the condition A does not hold
Next theorem helps to say when n long n ≫ 1 . We exclude the time t n 1 = 0 , i.e. the index n = 0 , since the ratio n long n for n = 0 is indeterminate 0 0 .

n ≫ 1 for all time
Here, "for all time" we do not mean for all times t n 1 , since in our analysis we consider t n 1 up to . So, we introduce the following definition This definition is made more precise by using a monitor function.

The condition B
The next theorem explains what happens under the condition B.
The previous theorem with c = 1 gives the following important results.  In the stiff situation, it is expected that if B holds, then (5.5) holds. In fact, E 1 ≪ 1 and it is expected | | non-small. So, we can state the following important conclusion.

The condition C
Although it is expected that C does not hold, we study anyway the condition C since it characterizes the transition between n ≈ long n in the long-time and n long n ≫ 1 for all time.
For the condition C, we need a weak form of n long n ≫ 1 for all time. Here is the definition with a monitor function. .

S(x)
e Re( i )t n 1 − e j t n 1 t n 1 ≤ 1 − e j t n 1 t n 1 .
e Re( i )t n 1 − e j t n 1 t n 1 = 1 − e j t n 1 t n 1 .
for t n 1 ∈ (0, S] . In particular, for S = v , in (5.7) we have whenever ≥ 1 . Then The previous theorem with c = 1 gives the following results. Suppose to be in the stiff situation and suppose that C holds and E 1−v 1 ≪ 1 . Then it is expected that (5.8) holds. So, we can state the following conclusion.
Suppose to be in the stiff situation and suppose that C holds. It is expected S-weakly n long n ≫ 1 for all time.

Examples revisited
Now, we look at the two examples of Sect. 2 in the light of the results of Sects. 4 and 5.

Same approximant with different ODEs
The conditions A, B and C, are W < 1 , W > 1 and W = 1 , respectively, where − 15

The region R hr 1
The condition A can be written as The region R −0.1 for the two methods is shown in Fig. 10. In the left part of the figure, we see that the region for the Gauss RK method does not cover points with large imaginary part on the line On the other hand, in the right part, we see that the region for the Radau RK method completely includes this line.

Independence of the non-rightmost spectrum
In this section, we study when the condition A holds independently of the particular non-rightmost spectrum − . Here, we consider an analytic approximant R with domain D such that {z ∈ ℂ ∶ Re(z) < R } ⊆ D for some R ∈ (0, +∞] , i.e. D includes a left half-plane.

The property A(x)
We introduce the property A(x) of the approximant R.
where def ⟺ has the meaning of "if and only if" by definition.
The property A(x) can be also written as where R x is the region defined in (4.18). Observe that A(0) is the A-stability property.
The property A(x) is important because A(hr 1 ) implies the condition A for all nonrightmost spectra − .
It is of interest to consider the property A(x) for |x| non-large. In fact, |hr 1 | ≤ h 1 and we are assuming h 1 non-large.
We have the following negative result.

Theorem 7.1
There exists x 0 > 0 such that, for x < R with |x| ≤ x 0 and x ≠ 0 , A(x) is not true.
Proof Remind that l is the order of the approximant R. In the complex plane, there exists a small disk centered at the origin which consists of l + 1 sectors of width l+1 included in the order star S , intercalated with l + 1 sectors of width l+1 included in S c . Thus, there exists x 0 > 0 such that, for x < R with |x| ≤ x 0 and x ≠ 0 , the line Re(z) = x has a non-empty intersection with the order star S . Let w be a point in this intersection. We have Re(w) = x and Then, due to the continuity of R, there exists > 0 such that, for any z ∈ ℂ with x − ≤ Re(z) ≤ x and Im(z) = Im(w) , we have e −x |R(z)| > 1. The properties A(x, a) and B(x, a) can be also written as The properties A(x, a) and B(x, a) are important because A(hr 1 , a) implies the condition A for all non-rightmost spectra − such that h − ≥ a and B(hr 1 , a) implies the condition B for all non-rightmost spectra − such that h − ≥ a . Remind that − and − are defined in (1.5).

The limit L
Now, we assume that exists. In addition, we also assume the following.
-When L < +∞ : where k > 0 , and, for any x < R and D ≥ 0 , where k > 0 , and, for any x < R and D ≥ 0 , where C = C(x, D) > 0. The next two subsections consider, for x < R , the cases L > e x and L < e x .
A(x, a) def ⟺ e −x |R(z)| < 1 for all z ∈ ℂ such that Re(z) < x and |z| ≥ a B(x, a) def ⟺ e −x |R(z)| > 1 for all z ∈ ℂ such that Re(z) < x and |z| ≥ a.
(7.5) e −x |R(z)| ≤ for all z ∈ ℂ such that Re(z) < x and |z| ≥ a.  If, in addition, L = L inf , then is not smaller than this negative number and can be arbitrarily close to it.

Approximants with L = +∞
Consider approximants with L = +∞ . Examples of such approximants are Taylor approximants and superdiagonal Padé approximants. The results in Sect. 7.5 say that the condition B holds for h − sufficiently away from zero, as confirmed in the first example of Sect. 2. In particular, B holds for As h − → +∞ , B holds with → +∞.

Approximants with L = 0
Consider approximants with L = 0 . Examples of such approximants are subdiagonal Padé approximants. Radau e Lobatto IIIC RK methods correspond to the first and second subdiagonal Padé approximants, respectively.
The results in Sect. 7.6 say that the condition A holds for h − sufficiently away from zero. In particular, A holds for As h − → +∞ , A holds with → −∞ . Moreover, Theorem 7.1 says that, for any rightmost real part r 1 ≠ 0 with |hr 1 | ≤ x 0 , we cannot have that A holds for all h − .
A-stable approximants with L = 0 are called L-stable (see [5]) and they are considered particularly suitable for integrating very stiff ODEs (see [1,3,8,11]). Observe that here we are also considering approximants with L = 0 which are not A-stable. Indeed, the A-stability property does not play a crucial role in this context. Among subdiagonal Padé approximants, only the first and second subdiagonal Padé approximants (Radau and Lobatto IIIC methods) are A-stable.

Approximants with L = 1
Consider approximants with L = 1 . Examples of approximants with L = 1 are diagonal Pad é approximants, which are also A-stable. Gauss methods correspond to the diagonal Padé approximants.
Suppose r 1 < 0 . The results in Sect. 7.5 say that the condition B holds for h − sufficiently away from zero, as confirmed in the second example of Sect. 2. In particular, B holds for For an A-stable approximant, B holds with ≤ − r 1 1 and, as h − → +∞ , → − Suppose r 1 > 0 . The results in Sect. 7.6 say that the condition A holds for h − sufficiently away from zero. In particular, A holds for For an A-stable approximant, A holds with ≥ −

Non-significant eigenvalues II
In this subsection we study when any non-rightmost eigenvalue i with |w i | ≥ 1 is non-significant (see Sect. 7.2). The importance of the region P x is due to the fact that, for a non-rightmost eigenvalue i , we have |w i | ≥ 1 if and only if h i ∈ P hr 1 .

The number a(x)
In other words, a(x) is the infimimum of the radii of open disks centered at the origin and including P x .
The importance of the number a(x) is given by the following theorem.

Theorem 7.8 For a non-rightmost eigenvalue
Proof The closed disk of radius a(x) centered at the origin includes the region P x . The theorem follows by reminding that P x contains the non-rightmost eigenvalues i such that |w i | ≥ 1 . ◻

The theorem on the non-significant eigenvalues
Next theorem says when any non-rightmost eigenvalue i with |w i | ≥ 1 is non-significant. It involves the behavior of a(x) as x → 0.

and then
The theorem now follows by reminding the definition of non-significant eigenvalue. By the previous theorem we obtain the following important conclusion.

Conclusion 7.10
Suppose that the approximant satisfies (7.7). It is expected that A holds.
In fact, suppose A does not hold, i.e. there is a non-rightmost eigenvalue i with |w i | ≥ 1 . It is expected that this eigenvalue is significant. On the other hand, if it is significant, then, by the previous theorem, we obtain that (7.8) does not hold and this is "unlikely".
In the next subsection, we show that the implicit Euler method satisfies (7.7).

The implicit Euler method
We examine the property A(x) and determine the number a(x) for the the implicit Euler method, corresponding to the (0, 1)-Padé approximant This approximant has R = 1.
The region R x , x < 1 , for this approximant is the exterior of the disk of center 1 and radius e −x and the region P x is the part of the closed disk at the left of the line Re(z) = x (see Fig. 11).
Theorem 7.11 Let x < 1 . For the implicit Euler method, we have A(x) if and only x = 0 . Moreover, we have Proof When x = 0 , A(x) is the A-stability. When x ≠ 0 , we have 1 − e −x < x and then since the complementary set R c x of R x is the closed disk of center 1 and radius e −x (see Fig. 10). Thus A(x) is not true.
For the second part, let b ≥ 0 . An easy computation shows that, for z ∈ ℂ such that |z| = b , we have Hence if and only if For x > 0 , (7.10) is equivalent to For x < 0 , (7.10) is equivalent to Now, equation (7.9) follows. ◻ By (7.9), we obtain We can conclude that it is expected that A holds for the implicit Euler method.

Conclusions
In the stiff situation, we have studied the long-time behavior of the relative error in the numerical integration of the ODE (1.1) with A normal. The numerical integration is accomplished over a mesh of constant stepsize h, by using at any step of an analytic approximant R of the exponential: see (1.2). The relative error n of the numerical integration is given in (1.3).
� ≠ z ∈ ℂ ∶ z ∈ P x and |z| = b = z ∈ ℂ ∶ Re(z) < x and |z| = b and z ∈ R c x = z ∈ ℂ ∶ |z| = b and 1 2 b 2 + 1 − e −x ≤ Re(z) < x (7.10) 1 2 b 2 + 1 − e −x < x and x > −b and We have defined the long-time solution y long as the solution of (1.1) projected on the eigenspace of the rightmost eigenvalues and we have considered the relative error long n of the numerical integration of y long . The error long n grows linearly in time, it is small and it remains small in the long-time.
We have introduced the condition where r 1 is the real part of the rightmost eigenvalues of A. When A holds, we have n ≈ long n in the long-time. When A does not hold, we have n long n ≫ 1 for all time. Let L = lim z→∞ |R(z)| . In order to have the condition A satisfied, it is better to use approximant with L = 0 (for example Radau and Lobatto IIIC methods). Approximants with L = 1 (for example Gauss methods) does not work well when r 1 < 0.
The paper [10] analyzes the numerical integration in the stiff situation by looking to a different question. In [10], the interest is about numerical approximations (1.2) of the long-time solution starting with a perturbed initial value. The approximants are analyzed by means of their error growth function R (see [4,5]) in order to study how they propagate the initial perturbation from the relative error point of view. In this other context, we have a non-large propagation of the initial perturbation if and only if We have considered the case of A normal. Some numerical experiments, not included here, suggest that also for non-normal matrices we have n ≈ long n in the long-time when the condition A holds and n long n ≫ 1 for all time when A does not hold. In light of this, the results of Sect. 7 becomes more important, since they are about the condition A.
We conclude by remarking that the findings of this paper are interesting in applications involving differential models described by linear ODEs with r 1 ≠ 0 . In particular, they are interesting when we are integrating an ODE whose solution decreases to small orders of magnitude (case r 1 < 0 ), but it is not yet considered as zero, or grows up to a large orders of magnitude (case r 1 > 0 ), but it is not yet considered as infinite.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.
Funding Open access funding provided by Università degli Studi di Trieste within the CRUI-CARE Agreement.
A: |R(h )| < e hr 1 for any non-rightmost eigenvalue of A,