CONDITIONING AND RELATIVE ERROR PROPAGATION IN LINEAR AUTONOMOUS ORDINARY DIFFERENTIAL EQUATIONS

. In this paper, we study the relative error propagation in the solution of linear autonomous ordinary diﬀerential equations with respect to per- turbations in the initial value. We also consider equations with a constant forcing term and a nonzero equilibrium. The study is carried out for equations deﬁned by normal matrices.


Introduction. Consider the scalar real Ordinary Differential Equation (ODE)
y (t) = ay (t) , t ≥ 0, y (0) = y 0 where a is a given real number. If the initial value y 0 = 0 is perturbed to y 0 with relative error ε = | y 0 − y 0 | |y 0 | , the value y (t) = e at y 0 of the solution at the time t is perturbed to y (t) = e at y 0 with relative error Whereas the relation between the absolute errors | y (t) − y (t)| and | y 0 − y 0 | is | y (t) − y (t)| = e at | y 0 − y 0 | , t ≥ 0, the relation between the relative errors δ(t) and ε is So, whereas the absolute error of the perturbation grows or decays with time in dependence of the sign of a, the relative error does not change with time, for all a, y 0 = 0 and y 0 . This shows that the absolute and relative errors of the perturbed solution behave very differently.

STEFANO MASET
Now, we want to consider the same question for the n-dimensional real linear autonomous ODE y (t) = Ay (t) , t ≥ 0, where A is a given n × n real matrix. Let · be a norm on R n . If the initial value y 0 ∈ R n \ {0} is perturbed to y 0 with relative error ε = y 0 − y 0 y 0 , the value y (t) = e tA y 0 of the solution at the time t is perturbed to y (t) = e tA y 0 with relative error Observe that y 0 = 0 implies, for all t ≥ 0, y (t) = 0 and so δ (t) is defined. As for the absolute errors y (t) − y (t) and y 0 − y 0 , the following two inequalities are well known (see e.g. [8]): where M (t) grows polynomially with t and α (A) is the spectral abscissa of A defined as the maximum real part of the eigenvalues of A, and There are specific ways for computing µ(A) in the usual cases where the norm . is the 1-norm, the 2-norm or the ∞-norm (see e.g. [5]). For example, in case of the 2−norm, the relevant logarithmic norm µ 2 (A) is the maximum eigenvalue of the symmetric part 1 2 A T + A of A. Despite its basic importance, it seems that the relation between the relative errors δ (t) and ε has not been yet investigated in literature. Indeed, any perturbation analysis for the ODE (1) found in literature considers the relative error of the perturbed matrix exponential e t A with respect to e tA , when A is perturbed to A: see [8], [4], [6], [1] and [7]. Moreover, in these papers the role of the particular initial value y 0 is ignored since they deal with the matrix e tA rather than the vector e tA y 0 .
In Section 2 below, we study the relation between δ (t) and ε by using the 2−norm as norm · and by assuming that A is a normal matrix. In the successive Section 3, the same analysis is accomplished for a linear autonomous ODE with a forcing term and a nonzero equilibrium.
Observe that the relative errors δ (t) and ε are measured on vectors, not on their components. However, for a vector u ∈ R n with nonzero components and a perturbed vector u ∈ R n , we have whenever the norm · is such that |v i | ≤ v , v ∈ R n and i = 1, . . . , n, v = (|v 1 | , . . . , |v n |) , v ∈ R n , v ≤ w , v, w ∈ R n with |v i | ≤ |w i | for i = 1, . . . , n.
This holds for a p-norm. Thus, relative error on vectors and relative errors on components are strictly related. For an initial value y 0 with nonzero components y 0i perturbed to y 0 with components y 0i : if i.e. all components of y 0 are perturbed with relative errors within a tolerance TOL, then ε ≤ TOL; if (an estimate of) δ(t) is known, then we have information on the relative errors of the components of y(t) by %bigskip In numerical analysis, the term "conditioning" refers to how a perturbation in the data of a problem influences the result of the problem. Conditioning studies consider relative errors, rather than absolute errors, because relative errors quantify the loss of accuracy, in terms of significant figures, introduced by a perturbation. The magnification of the relative error in the result with respect to the relative error in the data is ruled by the so-called "condition numbers".
In the next sections, we present a conditioning analysis with respect to the initial value of linear autonomous ODEs.
2. Condition numbers with respect to the initial value. We begin by presenting an illustrating example.
In Figure 1 left, we see the graph of the absolute error y (t) − y (t) 2 for t ∈ [0, 3]. The values y (t) and y (t) are computed by using the MATLAB function expm for the matrix exponential. We observe a fast decrease of the absolute error in time. The decrease can be explained by noticing that the eigenvalues λ 1 and λ 2 of A are negative: λ 1 = −0.01 and λ 2 = −5. Observe that we have  since µ 2 (A) = λ 1 , but (3) is not sharp in this situation.
In Figure 1 right, we see the graph of the relative error δ (t) = for t ∈ [0, 3]. Also the relative error has a fast decrease in time.
In Figure 2, we see the graphs of the absolute and relative errors on the interval  The left side shows a slow decrease of the absolute error, more coherent with the inequality (3). On the other hand, the right side shows a fast grow of one order of magnitude for the relative error. Indeed, at t = 3 we have y (t) = (0.0970, 0.0970) .  So, whereas the absolute error is more or less equal to the initial one as it is shown in Figure 2, y(t) is almost 10 times smaller than the initial value y 0 . This explains the growth of one order of magnitude for the relative error. The relative error can behave even worse. Consider the new initial value The graphs of the absolute and relative errors on the interval [0, 3] are presented in Figure 3. Along with a very slow decrease of the absolute error very similar to that in Figure 2, we have an explosion of six order of magnitude of the relative error. Indeed, at t = 3, we have y (t) = 10 −6 · (0.3059, 0.3059) with an absolute error more or less equal to the initial one.
In conclusion, we could say that although the matrix (2) is universally considered "stable", in the sense that it dumps a perturbation in the initial value, at the same time it can be considered "unstable" because it can magnify the same perturbation when we look at the relative error rather than the absolute error.
The previous example shows that the time propagation of the relative error in the initial value of the ODE (1) is an interesting phenomenon, which merits to be investigated. To this aim, now we introduce some more or less known notions about the conditioning of linear problems (see e.g. [2]).
Consider the problem of computing the result v = Bu for the data u ∈ R n \ {0}, where B ∈ R n×n is given and non-singular. Let · be a fixed vector norm on R n and with the same symbol we denote the matrix norm on R n×n induced by this vector norm. If the data u is perturbed to u with relative error By writing u = u + ε u z, where z = u − u u − u is the direction of the perturbation ( z is an arbitrary unit vector when u = u, i.e. ε = 0), we have u is called the condition number of B at the data u with direction of the perturbation z.
In general, we do not know the specific direction of the perturbation z, but we have information only on ε. So, we can only write where κ (B, u) := max is called the condition number of B at the data u. If the direction of the perturbation z is such that we have κ (B, u) = κ (B, u, z) and equality in (4) holds.
If we do not know the specific data u, we can only write where κ (B) := max we have κ (B) = κ (B, u) and if, in addition, the direction of the perturbation z satisfies (5), we have κ (B) = κ (B, u) = κ (B, u, z) and equality in (6) holds. We define the following condition numbers of the ODE (1) with respect to the initial value by looking at the result y (t) = e tA y 0 relevant to the data y 0 : K (t, A, y 0 ) := κ e tA , y 0 = max z0∈R n z0 =1 K (t, A, y 0 , z 0 ) = e tA e tA y 0 (8) where z 0 = y 0 − y 0 y 0 − y 0 and y 0 = y 0 y 0 .
Apart from the well-known condition number (9) of the matrix exponential, which turns out to be not too much important since it describes the propagation of the relative error in a nongeneric situation for the initial value, at the best of our knowledge conditions numbers for ODEs like (7) and (8) have not been previously considered in literature.

2-norm and normal matrices.
In the following, we use the 2-norm and we assume that A is a normal matrix. When the 2-norm is used, the previous condition numbers K ( · ) are written K 2 ( · ).
Under the assumption that A is diagonalizable, A has the spectral decomposition and, for an analytic complex function g and u ∈ R n , we have where λ 1 , . . . , λ p are the distinct complex eigenvalues of A and, for i = 1, . . . , p, P i is the projection on the eigenspace of the eigenvalue λ i . Under the stronger assumption that A is normal, the eigenspaces are orthogonal and so we have The big advantage of considering normal matrices is that then we have (11) at our disposal. On the other hand, by confining our attention to normal matrices is not too much restrictive since the class of normal matrices includes important families of matrices as symmetric matrices, skew-symmetric matrices and orthogonal matrices. Also observe that the choice of considering only normal matrices is adequate for a first paper on the subject as the present one.
The following theorem gives expressions for the three condition numbers K 2 (t, A, y 0, z 0 ), K 2 (t, A, y 0 ) and K 2 (t, A) in case of a normal matrix A.
Theorem 2.2. Let A ∈ R n×n be a normal matrix. Assume that the spectrum {λ 1 , . . . , λ p } of A is partitioned in the subsets Figure 4). For i = 1, . . . , p, let P i be the orthogonal projection on the eigenspace of the eigenvalue λ i . Moreover, for j = 1, . . . , q, let be the orthogonal projection on the sum of the eigenspaces of λ ij−1+1 , λ ij−1+2 , . . . , λ ij . Finally, let y 0 ∈ R n \ {0} be an initial value of (1) and let y 0 ∈ R n be a perturbed initial value.

STEFANO MASET
For any t ≥ 0, we have Moreover: • if z 0 lies in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , . . . , λ i1 with maximum real part, then, for any t ≥ 0, • if y 0 lies in the sum of the eigenspaces of the eigenvalues λ iq−1+1 , λ iq−1+2 , . . . , λ iq = λ p with minimum real part, then, for any t ≥ 0, and if, in addition, z 0 lies in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , . . . , λ i1 , then, for any t ≥ 0, Proof. By (11) we have, for u ∈ R n , Thus, for t ≥ 0, we have . Now, we give the expression for K 2 (t, A, y 0 ). Observe that Moreover, if Q j z 0 = 0, j = 2, . . . , q, i.e. z 0 lies in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , . . . , λ i1 , then for such z 0 . Finally, we give the expression for K 2 (t, A). Observe that Moreover, if Q j y 0 = 0, j = 1, . . . , q − 1, i.e. y 0 lies in the sum of the eigenspaces of the eigenvalues λ iq−1+1 , λ iq−1+2 , . . . , λ iq = λ p , then Thus For a normal matrix A, the previous theorem says that if the initial value y 0 ∈ R n \ {0} of (1) is perturbed to y 0 with direction of the perturbation z 0 and relative error ε, then, regarding the relative error δ of the perturbed solution, we have the following three facts.
• For any t ≥ 0, we have • For any t ≥ 0, independently of z 0 , we have with ≤ replaced by = for any t ≥ 0 if z 0 lies in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , . . . , λ i1 . • For any t ≥ 0, independently of y 0 and z 0 , we have with ≤ replaced by = for any t ≥ 0 if y 0 lies in the sum of the eigenspaces of the eigenvalues λ iq−1+1 , λ iq−1+2 , . . . , λ iq = λ p and z 0 lies in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , . . . , λ i1 . Regarding the condition number K 2 (t, A, y 0 ) given in (13), we can say that it is an increasing function of t with K 2 (0, A, y 0 ) = 1 whose asymptotic behavior is given by where j * is the minimum index j = 1, . . . , q such that Q j y 0 = 0 and ∼ means limit 1 as t → +∞. Regarding the condition number K 2 (t, A, y 0 , z 0 ) given in (12), we can write R 2 (t, A, z 0 ) is a decreasing function of t with R 2 (0, A, z 0 ) = 1 whose asymptotic behavior is given by where j * * is the minimum index j = 1, . . . , q such that Q j z 0 = 0. So, the asymptotic behavior of K 2 (t, A, y 0 , z 0 ) is given by i.e. q = 1, then for all y 0 = 0 and y 0 as in the scalar case. Observe that skew-symmetric matrices fulfill (16). It is not a surprise that (17) holds for skew-symmetric matrices: if A is skew-symmetric, then , for t ≥ 0, e tA is orthogonal and so δ(t) = e tA ( y 0 − y 0 ) 2 e tA y 0 2 = y 0 − y 0 2 y 0 2 = ε.
Thus j * = 1 with Q 1 y 0 2 = 0.2 √ 2.60 and j * * = 2 with Q 2 z 0 2 = 1. We conclude that and this explains the fast decrease of the relative error in the right side of Figure 1.
In Figure 2, we have the same y 0 of Figure 1 and Thus j * and Q 1 y 0 2 are the same as above and now j * * = 1 with Q 1 z 0 2 = 1. We conclude that and this explains the growth of one order of magnitude for the relative error in the right side of Figure 2. Observe that in this case we have K 2 (t, A, y 0 , z 0 ) = K 2 (t, A, y 0 ), t ≥ 0, since z 0 lies in the eigenspace of λ 1 . Finally, in Figure 3 we have y 0 = (1, −1) and y 0 = y 0 + (0.01, 0.01) .
Thus j * = 2 with Q 1 y 0 2 = 1 and j * * = 1 with Q 1 z 0 2 = 1. We conclude that and this explains the explosion of the relative error in the right side of Figure 3.
Observe that in this case we have since y 0 lies in the eigenspace of λ 2 and z 0 lies in the eigenspace of λ 1 and then δ (t) = e 4.99t · ε, t ≥ 0.
What happens in this last case is that the unperturbed solution includes only the fast decreasing term e λ2t = e −5t , but any small perturbation in the initial value moving it outside the eigenspace of λ 2 gives rise to the slow decreasing term e λ1t = e −0.01t in the perturbed solution that soon dominates completely the fast decreasing term.
By concluding this section, we remark that, for a normal matrix A ∈ R n×n , "all the eigenvalues of A have real part ≤ 0" is a necessary and sufficient condition for not amplifying the absolute error of perturbations in the initial value, i.e.
On the other hand, above we have seen that "all the eigenvalues of A have the same real part" is a necessary and sufficient condition for not amplifying the relative error of perturbations in the initial value, i.e.
Observe that the condition (18) for the absolute error is a condition on the absolute position of the eigenvalues in the complex plane, whereas the condition (19) for the relative error is a condition on their relative position.
Since the real parts of the eigenvalues of the normal matrix A are the eigenvalues of the symmetric part 1 2 A T + A of A (see [3]), we obtain that (19) holds if and only if 1 2 A T + A = αI for some α ∈ R, i.e. a ii = α, i = 1, . . . , n, a ji = −a ij , i, j = 1, . . . , n with i = j.
So, a normal matrix does not amplify the relative error of perturbations in the initial value if and only if it is the sum of a skew-symmetrix matrix and a multiple of the identity. We also remark that all our previous results are not confined to normal real matrices but they are also valid for normal complex matrices.
3. Condition numbers for equations with nonzero equilibrium. Many mathematical models based on linear autonomous ordinary differential equations have, unlike (1), a nonzero equilibrium. Therefore, now we consider the ODE where A ∈ R n×n is non-singular and b ∈ R n \ {0}. The solution is given by For y 0 = x, the ODE (20) has the equilibrium solution which is the unique equilibrium.
A relative error analysis similar to (1) can be done. We consider an initial value y 0 = 0 and a perturbed initial value y 0 . Let ε = y 0 − y 0 y 0 .

CONDITIONING IN LINEAR AUTONOMOUS ORDINARY DIFF. EQUATIONS 2893
Let y be the solution relevant to the pertubed initial value y 0 and let We study the relation between the relative errors δ (t) and ε.
Unlike (1), it is not true that y (t) = 0, t ≥ 0, for any y 0 = 0. We consider the relative error δ (t) not defined for t ≥ 0 such that y (t) = 0. Now, we introduce the two condition numbers J (t, A, b, y 0 , z 0 ) and J (t, A, b, y 0 ), whose role is similar to K (t, A, y 0 , z 0 ) and K (t, A, y 0 ), respectively. We have with z 0 = y 0 − y 0 y 0 − y 0 being the direction of perturbation ( z 0 is an arbitrary unit vector when y 0 = y 0 , i.e. ε = 0). Moreover, we set On the other hand, the relative error analysis for (20) is different from the relative error analysis for (1): in both cases we have where z 0 is the direction of the perturbation, but y(t) = e tA y 0 for (1) and y(t) = x + e tA (y 0 − x) for (20).

Let us introduce the unit vector
x = x x and then we write an initial value as where the unit vector d 0 and c ≥ 0 are given by ( d 0 is an arbitrary unit vector when y 0 = x, i.e. c = 0). Since y 0 = 0, the situation d 0 = − x and c = 1 is excluded.
We have and so and 3.1. 2-norm and normal matrices. Now, as in the relative error analysis for (1), we use the 2-norm (by writing the condition numbers J( · ) as J 2 ( · )) and we assume that A is a normal matrix. Next theorem gives an expression for the two condition numbers J 2 (t, A, b, y 0 , z 0 ) and J 2 (t, A, b, y 0 ). In the theorem and in the following, ·, · denotes the standard scalar product on R n . Theorem 3.1. Let A ∈ R n×n be a non-singular normal matrix. Assume that the spectrum of A is partitioned by decreasing real parts as in Theorem 2.2. We use the same notations of that theorem. Let y 0 ∈ R n \ {0} be an initial value of (20) and let y 0 ∈ R n be a perturbed initial value.
Regarding the function µ j , j = 1, . . . , q, in the previous theorem, the following bound holds.
Proof. Since x − Q j x is orthogonal to the sum of the eigenspaces of λ ij−1+1 , . . . , λ ij , we have 3.2. Asymptotic behavior of J 2 (t, A, b, y 0 , z 0 ) and J 2 (t, A, b, y 0 ). By (25) and (15), we obtain where j * * is the minimum index j = 1, . . . , q such that Q j z 0 = 0. Observe that (33) describes the asymptotic behavior of the condition number J 2 (t, A, b, y 0 , z 0 ) in terms of the asymptotic behavior of the condition number J 2 (t, A, b, y 0 ).
In the case y 0 = x, i.e. c = 0, we have This is in accordance with the fact that the equilibrium x is unstable, stable and asymptotically stable, when r 1 > 0, r 1 = 0 and r 1 < 0, respectively. Moreover, in this case we have In the following, we assume y 0 = x, i.e. c = 0. Moreover, let j * be the minimum index j = 1, . . . , q such that Q j d 0 = 0.
By (26), we obtain the asymptotic behavior of the condition number J 2 (t, A, b, y 0 ) and then, by (33), the asymptotic behavior of the condition number J 2 (t, A, b, y 0 , z 0 ) in the three cases r j * > 0, r j * < 0 and r j * = 0.
Observe that by Theorem 3.2.
Now, we consider some questions concerning the case r j * = 0. Then, we analyze the asymptotic behavior in the cases r 1 > 0, r 1 < 0 and r 1 = 0.
Theorem 3.3. Assume Q j * x = 0. Let d j * be the dimension of the sum the eigenspaces of λ i j * −1 +1 , . . . , λ ii j * . If d j * = 2, i.e. there is only one pair of conjugate pure imaginary eigenvalues and they are simple eigenvalues, then holds for infinitely many t ≥ 0 which repeat themselves periodically. If d j * > 2, the condition "(40) holds for some t ≥ 0" is a nongeneric condition on x and d 0 .
Proof. Let t ≥ 0. Observe that (see (32)) and then Moreover, we have orthogonal vectors of R n . Here, Re (w) and Im (w) denote the vectors whose components are the real parts and the imaginary parts, respectively, of the components of the vector w. Now, by looking at the Cauchy-Schwarz inequality used in the proof of Theorem 3.2, we see that Let V j * be the sum the eigenspaces of λ i j * −1 +1 , . . . , λ i j * and let d j * be its dimension. If d j * = 2, then (42) holds for infinitely many t ≥ 0 which repeat themselves periodically. In fact, Q j * x ∈ V j * and, in case of d j * = 2, we have i j * − i j * −1 = 2 and Re(P i j * d 0 ) and Im(P i j * d 0 ) constitute a base for V j * .
On the other hand, observe that is a curve on the intersection W of V j * and the sphere of radius Q j * d 0 2 (see (41)). Thus, we have "(42) holds for some t ≥ 0" if and only if the point of W is on the curve (43). If d j * > 2, since W is a manifold of dimension d j * − 1 > 1 and (44) can be any point in W by varying the equilibrium x, we conclude that the condition "the point (44) is on the curve (43)" is a nongeneric condition on x, fixed d 0 .
By the previous theorem, we obtain that if there is only one pair of conjugate pure imaginary eigenvalues and they are simple eigenvalues, then The quantity (see (22) and (30)) appears at the denominator in (38) and (39). It is a measure of how much the solution y is asymptotically away from zero. Next theorem says when it is nonzero and when it becomes zero.
Theorem 3.4. If Q j * x 2 < 1 and c Q j * d 0 2 = 1, then If Q j * x 2 = 1 and c Q j * d 0 2 = 1, then if and only if Proof. We have with (see (39)) Observe that if Q j * x 2 = 1 and c Q j * d 0 2 = 1, then the condition (46) is satisfied when there is only one pair of conjugate pure imaginary eigenvalues and they are simple eigenvalues. So, when there is only one pair of conjugate pure imaginary eigenvalues and they are simple eigenvalues, we have The asymptotic behavior of J 2 (t, A, b, y 0 , z 0 ) when j * = 1 is described by (35) and we have In the case Q j * d 0 = 0, i.e. j * > 1, by (34) or (36) or (38) we have lim t→+∞ J 2 (t, A, b, y 0 ) = +∞.
The asymptotic behavior of J 2 (t, A, b, y 0 , z 0 ) when j * > 1 is described by (35) or (37) or (39). For example, in the case r j * > 0 we have
The relative error in the initial value is asymptotically vanishing and this is in accordance with the fact that the equilibrium x is globally asymptotically stable when r 1 < 0: for any initial value, the absolute error is asymptotically vanishing and the solution approaches to x = 0 as t → +∞.

3.2.4.
Asymptotic behavior in the case r 1 = 0. Suppose r 1 = 0. Since A is nonsingular, the eigenvalues of A with zero real part are conjugate pure imaginary pairs of eigenvalues and i 1 , the number of such eigenvalues, is even. Let where ω 2 , ω 4 , . . . , ω i1 > 0, be such conjugate pure imaginary pairs of eigenvalues.
For an initial value y 0 ∈ R n such that y 0 − x is in the sum of the eigenspaces of λ k−1 and λ k , k = 2, 4, . . . , i 1 , we obtain as solution the equilibrium plus a periodic function of period 2π ω k : For an initial value y 0 ∈ R n such that y 0 − x is in the sum of the eigenspaces of the eigenvalues λ 1 , λ 2 , λ 3 , λ 4 , . . . , λ i1−1 , λ i1 , we obtain as solution the equilibrium plus a sum of i1 2 periodic functions of periods 2π ω1 , 2π ω4 , . . . , 2π : For a generic initial value y 0 ∈ R n , we obtain as solution the equilibrium point plus a sum of i1 2 periodic functions of periods 2π ω2 , 2π ω4 , . . . , 2π ωi 1 plus a function asymptotically vanishing.
Observe that: 1. If c → 0, i.e. y 0 → x, then The asymptotic behavior of J 2 (t, A, b, y 0 , z 0 ) when j * > 1 is described by (37) and we have The asymptotic behavior of J 2 (t, A, b, y 0 , z 0 ) when j * = 1 is described by (39) (28)). We see this in the right side of Figure 5, where the graph of J 2 (t, A, b, y 0 ) for t ∈ [0, 3T ] is shown. We have lim sup  4. Conclusions. In this paper, we have studied how a relative error in the initial value of the ODE (1) with zero equilibrium or the ODE (20) with a nonzero equilibrium is propagated along the solution. In particular, given a norm · on R n ,   where y is the solution relevant to the initial value y 0 = 0 and y is the solution relevant to the perturbed initial value y 0 .
We have confined our study to the 2-norm and this is not restrictive. In fact, given an arbitrary norm · different from the 2-norm, the norm equivalence says that there exists constants c, C > 0 depending on the dimension n of the ODE such