Optimal subsets in the stability regions of multistep methods

In this work we study the stability regions of linear multistep or multiderivative multistep methods for initial-value problems by using techniques that are straightforward to implement in modern computer algebra systems. In many applications, one is interested in (i) checking whether a given subset of the complex plane (e.g. a sector, disk, or parabola) is included in the stability region of the numerical method, (ii) finding the largest subset of a certain shape contained in the stability region of a given method, or (iii) finding the numerical method in a parametric family of multistep methods whose stability region contains the largest subset of a given shape. First we describe a simple procedure to exactly calculate the stability angle $\alpha$ in the definition of $A(\alpha)$-stability. As an illustration, we consider two finite families of implicit multistep methods: we exactly compute the stability angles for the $k$-step BDF methods ($3\le k\le 6$) and for the $k$-step second-derivative multistep methods of Enright ($3\le k\le 7$). Next we determine the exact value of the stability radius in the BDF family for each $3\le k\le 6$, that is, the radius of the largest disk in the left half of the complex plane, symmetric with respect to the real axis, touching the imaginary axis and lying in the stability region of the corresponding method. Finally, we demonstrate how some Schur--Cohn-type theorems of recursive nature and not relying on the RLC method can be used to exactly solve some optimization problems within infinite parametric families of multistep methods. As an example, we choose a two-parameter family of implicit-explicit (IMEX) methods: we identify the unique method having the largest stability angle in the family, then we find the unique method in the same family whose stability region contains the largest parabola.

First we describe a simple procedure to exactly calculate the stability angle α in the definition of A(α)-stability by representing the root locus curve (RLC) of the multistep method as an implicit algebraic curve.As an illustration, we consider two finite families of implicit multistep methods.We exactly compute the stability angles for the k-step BDF methods (3 ≤ k ≤ 6) and discover that the values of tan(α) are surprisingly simple algebraic numbers of degree 2, 2, 4 and 2, respectively.In contrast, the corresponding values of tan(α) for the k-step second-derivative multistep methods of Enright (3 ≤ k ≤ 7) are much more complicated; the smallest algebraic degree here is 22.
Next we determine the exact value of the stability radius in the BDF family for each 3 ≤ k ≤ 6, that is, the radius of the largest disk in the left half of the complex plane, symmetric with respect to the real axis, touching the imaginary axis and lying in the stability region of the corresponding method.These radii turn out to be algebraic numbers of degree 2, 3, 5 and 5, respectively.
Finally, we demonstrate how some Schur-Cohn-type theorems of recursive nature and not relying on the RLC method can be used to exactly solve some optimization problems within infinite parametric families of multistep methods.As an example, we choose a two-parameter family of implicit-explicit (IMEX) methods: we identify the unique method having the largest stability angle in the family, then we find the unique method in the same family whose stability region contains the largest parabola.

Introduction
In the stability theory of one-step or multistep methods for initial-value problems, one is often interested in various geometric properties of the stability region S ⊂ C of the method.In this work we study the shape of the stability region of linear multistep methods (LMMs) or multiderivative multistep methods (also known as generalized LMMs) as follows.
Suppose we are given a) a stability region S, or b) a family of stability regions S β parametrized by some β ∈ R d , and a family of subsets of C, denoted by F. Due to their relevance in applications, we will consider the following three classes: is the family of infinite sectors in the left half of C, with vertex at the origin, symmetric about the negative real axis, and parametrized by the sector angle α ∈ (0, π/2); is the family of disks in the left half of C, symmetric with respect to the real axis, touching the imaginary axis, and parametrized by the disk radius r > 0; • F = F para m is the family of parabolas in the left half of C, symmetric with respect to the real axis, touching the imaginary axis, and parametrized by some m > 0.
Our goal is to find the set H ∈ F with the largest parameter (α, r, or m) such that • H ⊂ S in case a); • H ⊂ S βopt for some stability region in the family in case b), but H ⊂ S β for β = β opt .
We will present some tools to handle these shape optimization questions, and, as an illustration, exactly solve some of them by using Mathematica version 11 in the BDF and Enright families (as LMMs and multiderivative multistep methods, respectively), and in an infinite family of IMEX methods with d = 2 parameters.

Motivation and main results
A When solving stiff ordinary differential equations, one desirable property of the numerical method is A-stability: a method is A-stable if the closed left half-plane {z ∈ C : Re(z) ≤ 0} belongs to S. Many useful methods are not A-stable, still, S contains a sufficiently large infinite sector in the left half-plane with vertex at the origin and symmetric about the negative real axis.This leads to the notion of A(α)-stability: a method is A(α)-stable with some 0 where the argument of a non-zero complex number satisfies −π < arg ≤ π.The largest 0 < α < π/2 such that (1) holds is referred to as the stability angle of the method [17].Various other stability concepts-such as A(0)-stability, A 0 -stability, • A-stability, stiff stability, or asymptotic A(α)-stability-have also been defined, and theorems devised to test whether a given multistep method is stable in one of the above senses; see, for example, [34,9,20,21,22,23,26,11,3,40].There are various techniques to test A(α)-stability for a given α value.
In [3], for example, the sector on the left-hand side of ( 1) is decomposed into an infinite union of disks, and a bijection between each disk and the left half-plane is established via fractional linear transformations to employ a Routh-Hurwitz-type criterion.Another way of studying A(α)-stability is to consider the root locus curve (RLC) of the multistep method [17].Based on the RLC and some theorems from complex analysis, [38] presents a criterion for a LMM to be A(α)-stable for a given α; the stability angle is then obtained as the solution of an optimization problem involving Chebyshev polynomials.The procedure in [38] is formulated only for LMMs but not for multiderivative multistep methods.
The first goal of the present work is to describe an elementary approach to exactly determine the stability angle of a LMM or multiderivative multistep method: by eliminating the complex exponential function from the RLC and using a tangency condition, a system of polynomial equations in two variables is set up whose solution yields the stability angle.This process is easily implemented in computer algebra systems.As an illustration, we consider two finite families: the BDF methods [13,17,38] as LMMs, and the second-derivative multistep methods of Enright [10,6,17].With α BDF k denoting the stability angle of the k-step BDF method for 3 ≤ k ≤ 6, we show that tan α BDF k is an unexpectedly simple algebraic number, having degree 2 for k ∈ {3, 4, 6}, and degree 4 for k = 5; see Table 1.For the k-step Enright methods with 3 ≤ k ≤ 7, the corresponding constants tan α Enr k (with approximate values listed in Table 2) are much more complicated algebraic numbers of increasing degree (starting with 22).As far as we know, exact values α ∈ (0, π/2) for the stability angles of multistep methods were not presented earlier in the literature.
Remark 1.1.The k-step BDF methods for k ∈ {1, 2} are A-stable.For k ≥ 7 they are not zero-stable [8,7,16], therefore not interesting from a practical point of view.
Remark 1.2.In [38, Table 1] one finds some approximate values for the BDF stability angles, however, some of these values are not correct.The k = 3 value is wrong because the polynomial R 3 is not computed properly.The approximate values for k = 4 and k = 5 given in [38,Table 1] are correct (up to the given precision).The value for k = 6 is again incorrect because an error was committed in the minimization process.If the optimization in [38,Section 3] is carried out exactly with the correct R j polynomials, we recover the stability angle values in our Table 1.The errors in [38,Table 1] propagated in the literature, see, for example, [33, p. 242].As a consequence, some works that appeared in the current millennium also contain the erroneous angles.In [17, Chapter V.2, (2.7)] the correct approximate values are presented.
Remark 1.3.At the time of writing this document, we learned (through personal communication) that [1] also contains the exact stability angles for the BDF methods with 3 ≤ k ≤ 6 steps: although they use a different technique to derive the results and the arcsin function to express the final constants, the values given in [1] and our Table 1 are the same.Notice, however, that the stability angle for k = 5 given in [1] has a slightly more complicated structure than the value in our Table 1.
Remark 1.4.The k-step Enright methods are A-stable again for k ∈ {1, 2}, see [17], and unstable for k ≥ 8.More precisely, [11] proves that these methods are not A 0 -stable for k ≥ 8, hence they cannot be stiffly stable either, see [23,Theorem 3] (cf.[22,26]).However, in [17 Remark 1.5.For LMMs (and for more general methods as well), various other step-size coefficients have been introduced in the context of linear or non-linear problems.These coefficients govern the largest allowable step-size guaranteeing certain monotonicity or boundedness properties of the LMM, including the TVD and SSP properties [15].These properties are relevant, for example, in the time integration of method-of-lines semi-discretizations of hyperbolic conservation laws [41,19,35].
Remark 1.6.In [31] the largest inscribed and smallest circumscribed (semi)disks are computed for certain one-step methods.
The second goal of the present work is to compute the stability radius for some multistep methods.We will achieve this by using again the algebraic form of the RLCs.Table 3 contains the exact values in the BDF family for 3 ≤ k ≤ 6.
C The RLC, as the graph of a [0, 2π] → C function (or a union of such functions for generalized LMMs), yields information about the boundary of the stability region, ∂S.It is known, however, that in general the RLC does not coincide with ∂S (see Figure 3).This does not pose a problem when a fixed multistep method is considered-one can evaluate the roots of the characteristic polynomial at finitely many test points sampled from different components of C determined by the RLC to see which component belongs to S and which one to C \ S.But when working with parametric families of multistep methods, the precise identification of the stability region boundaries or components can become challenging with the RLC method.One can overcome this difficulty for example by invoking a reduction process, the Schur-Cohn reduction, formulated in e.g.[37].Instead of using auxiliary fractional linear transformations and applying Routh-Hurwitz-type criteria [36,28] as mentioned above, these Schur-Cohn-type theorems in [37] are directly tailored to the context of multistep methods to locate the roots of the characteristic polynomials with respect to the unit disk.
The third goal of the present work is to demonstrate the effectiveness of the Schur-Cohn reduction when we solve two optimization case studies in a family of implicit-explicit (IMEX) multistep methods taken from [18].On the one hand, we find the method in the IMEX family that has the largest stability angle, that is, the method whose stability region contains the largest sector (see our Theorem 5.3).On the other hand, we illustrate the versatility of the reduction technique by also finding the method whose stability region contains the largest parabola (see Theorem 6.1); the inclusion of a parabola-shaped region in S is relevant when studying semi-discretizations of certain partial differential equations (PDEs) of advection-reaction-diffusion type [28,5,18].The chosen IMEX family is described by two real parameters, and the corresponding characteristic polynomial is cubic.The Schur-Cohn reduction process recursively decreases the degree of the characteristic polynomial, so instead of analyzing the roots of high-degree polynomials, we finally need to check polynomial inequalities in the parameters present in the coefficients of the original polynomial.Besides the two real parameters, two complex variables are involved in our calculations-the nontrivial interplay between these six real variables determines the optimum in both cases.We emphasize that we solve the optimization problem exactly, and RLCs are not relied on in the rigorous part of the proofs (only when setting up conjectures about the optimal values).
Remark 1.7.The Schur-Cohn reduction is also used in [25] to explore certain properties of a discrete parametric family of multistep methods.Conditions for disk or segment inclusions in the stability regions of a two-parameter family of multistep methods are formulated in [39].Optimality questions about the size and shape of the stability regions of one-step or multistep methods are investigated in detail in [27].Properties of optimal stability polynomials and stability region optimization in parametric families of one-step methods are discussed, for example, in [29,30].

Structure of the paper
In Section 2.1, we introduce some notation.In Sctions 2.2-2.3,we review the Schur-Cohn reduction and the definition of the stability region of a multistep method.In Sections 2.4-2.5, the definition of the root locus curve is recalled in two special cases: for linear multistep methods and for second-derivative multistep methods.Here we consider the BDF and Enright families as concrete examples.
Regarding the new results, a simple algebraic technique is described in Section 3.1 to exactly compute the stability angle of a linear multistep or multiderivative multistep method.Stability angles for the BDF and Enright families are tabulated in Sections 3.2-3.3.In Section 4, we exactly compute the stability radii in the BDF family by using the same approach.In Section 5, we first describe a two-parameter family of IMEX multistep methods, in which we determine the unique method with the largest stability angle, then, in Section 6, the unique method whose stability region contains the largest parabola.The techniques in Sections 5-6 do not rely on root locus curves but use the Schur-Cohn reduction instead; the full proofs are deferred to Appendices A and B.

Notation
The set of natural numbers {0, 1, . ..} is denoted by N. For z ∈ C, Re(z), Im(z), and z denote the real and imaginary parts, and the conjugate of z, respectively, and i is the imaginary unit.The boundary of a (possibly unbounded) set H ⊂ C is ∂H ⊂ C. When describing certain algebraic numbers of higher degree, a polynomial n j=0 a j x j with a j ∈ Z, a n = 0 and n ≥ 3 will be represented simply by its coefficient list {a n , a n−1 , . . ., a 0 }.For a polynomial Q(z) = n j=0 a j z j with 0 ≤ n ∈ N, a j ∈ C (0 ≤ j ≤ n) and a n = 0, we denote its degree, leading coefficient and constant coefficient by deg Q = n, lc Q = a n and cc Q = a 0 .The acronyms RLC and LMM stand for root locus curve and linear multistep method, respectively.

The Schur-Cohn reduction
In the rest of this section we assume that Q is a univariate polynomial with deg Q ≥ 1, and follow the terminology of [37]-we have explicitly added the deg Q ≥ 1 condition, being implicit in [37].We say that The reduced polynomial of Q(z) = n j=0 a j z j is defined as When this reduction process is iterated, we write Q rr for (Q r ) r , for example.The following theorems from [37] use the notion of the reduced polynomial and the derivative to formulate necessary and sufficient conditions for a polynomial to be in the above classes.In all three theorems below it is assumed that lc Remark 2.6.Let us consider the following example when applying the theorems above, e.g.Theorem 2.4.For any λ > 0, we set Q λ (z) := z 2 + λiz + 1.Then the roots of . This shows that it can happen that the degree of the original polynomial is > 1, but its reduced polynomial is a non-zero constant, so the relation Q r ∈ vN is undefined.In these cases, when Q r is a non-zero constant, notice that neither |Q r | < 1, nor |Q r | = 1, nor |Q r | > 1 can help us in general to determine whether Q ∈ vN or not (of course, the other condition |lc Q| > |cc Q| is violated now); cf. the sentence above [37, Theorem 5.1].

The stability region of a multistep method
Stability properties of a broad class of numerical methods (including Runge-Kutta methods, linear multistep methods, or multiderivative multistep methods) for solving initial value problems of the form can be analyzed by studying the stability region of the method.When an s-stage k-step method (s ≥ 1, k ≥ 1 fixed positive integers; for k = 1 we have a one-step method, while for k ≥ 2 a multistep method) with constant step size h > 0 is applied to the linear test equation y = λy (λ ∈ C fixed, y(0) = y 0 given), the method yields a numerical solution (y n ) n∈N that approximates the exact solution y at time t n := t 0 + nh and satisfies a recurrence relation of the form The characteristic polynomial associated with the method takes the form With Φ(•, µ) abbreviating the polynomial ζ → Φ(ζ, µ), the stability region of the method is defined as Remark 2.7.Some other variations of the above definition of the stability region of a multistep method have also been proposed in the literature, see, e.g.[24].In [4, p. 344], the "open stability region" is defined as the set see also [44, p. 348], [12, p. 452] or [33].In e.g.[17,32], the stability region of the method (3) is defined as and multiple roots satisfy |ζ j (µ)| < 1}, that is, essentially, Φ(•, µ) ∈ svN.In [27, Formula (2.5)] the stability region is given by and if |ζ j | = 1, then it is a simple root}, with C denoting the extended complex plane.
We can regroup the terms in (4) as ζ with some suitable polynomials C .The inequality condition in (3) implies that the leading coefficient C k does not vanish identically; it may happen that for some exceptional µ values the leading coefficient is zero: For example, for the implicit Euler (IE) method (7)) is interpreted formally, we have for the IE method that E = {1} ⊂ S (because (6) is satisfied vacuously).Similarly, for the BDF2 method, E = {3/2} ⊂ S (because then the unique root of However, elements of E or E ∩ S can be problematic.(i) For µ ∈ E, the order of the recursion (3) decreases, thus, in general, the starting values y 0 , y 1 , . . ., y k−1 of the numerical method cannot be chosen arbitrarily.(ii) Some exceptional values µ ∈ E ∩ S can be located in the interior of the corresponding region of instability of the method-this is the case for example for both the IE and BDF2 methods.When the step size h > 0 is chosen in a way that µ ∈ E ∩ S is such an isolated value, the recursion (3) generated by the numerical method becomes practically useless (it quickly "blows up" for arbitrarily small perturbations of h).(iii) RLCs are often used to identify the boundary ∂S of the stability region (see Sections 2.4-2.5 below).In [27,Definition (2.21)], the RLC is given by It can happen that ∂S is a proper subset of the corresponding RLC (see, for example, our Figure 3), but in [27,Corollary 2.6] it is shown that for a numerical method satisfying Property C (see [27,Formula (2.9)] or [17,Definition 4.7]), the RLC coincides with ∂S.
According to [17, Section V.4], all one-step methods have Property C, so the IE method also has.And indeed, applying [27,Proposition 2.7] to the IE method we now have that and σ have no common root and /σ is univalent on the set {z ∈ ∈ Γ = ∂S.This apparent contradiction seems to indicate that the authors of [27] interpreted definition (7) intuitively: a root ζ = ∞ is tacitly introduced as soon as the leading coefficient C k (µ) becomes zero.So [27, Corollary 2.6], for example, actually relies on definition (5) rather than on definition (7) (or (6)).
Notice that, with the theorems cited in our Section 2.2, one can directly investigate the stability region of a numerical method, without constructing the corresponding RLC or without analyzing the relation between ∂S and the RLC (see Sections 5-6 below).
Finally we remark that the above considerations also play an important role, e.g. in control theory [2, Chapter 1], where a "degree invariance" (i.e., "no degree loss") condition is incorporated in the Boundary Crossing Theorem.[2, Chapter 1] also recalls several stability results for polynomials, e.g. the Routh-Hurwitz, Jury, or the recursive Schur(-Cohn) stability tests.

The RLC of a LMM
A linear multistep method for (2) has the form where f m := f (t m , y m ), and the numbers α j ∈ R and β j ∈ R (j = 0, . . ., k) are the suitably chosen method coefficients with α k = 0.The method is implicit, if β k = 0.By setting One way to study the stability region (5), or its boundary ∂S in the complex plane is to depict the RLC corresponding to the method [17]: observe that P 1 is linear in µ, so

RLCs for the BDF methods
Each member of the BDF family is a special case of ( 8).The k-step BDF method (having order k) is given by where ∇ denotes the backward difference operator ∇y n+1 := y n+1 − y n , and ∇ j y n+1 := ∇ j−1 y n+1 − ∇ j−1 y n (for j > 1).It is known [17] that the corresponding RLC is Figures 1-3 show the RLCs for some BDF methods.

The RLC of a multiderivative multistep method
A second-derivative multistep method is more general than (8) and can be written as where , and the method is determined by the coefficients α j , β j and γ j , see [17].Now the associated characteristic polynomial (4) becomes This time we have two RLCs: where µ 1,2 are the two solutions of P 2 e iϑ , µ = 0.For any choice of the method coefficients α j , β j and γ j , one can construct µ 1,2 explicitly, since P 2 is only quadratic in µ.

RLCs for the Enright methods
The Enright methods are special cases of (12), and for k ≥ 1 they are defined [17] as where with the usual extension of the binomial coefficients.From ( 14) one obtains the RLCs of the Enright methods, see Figures 4 and 5.The order of the k-step Enright method is k + 2. 3 Optimal sector inclusions

The RLC in implicit algebraic form
Computing the stability angle of a method with stability region S is equivalent to finding the slope of the unique line L that passes through the origin, touches ∂S at some point in the open upper left half-plane such that ∂S lies on the right-hand side of L (viewed from the origin) in this quadrant.This last requirement is necessary since ∂S ∩ L can consist of more points, even in the open upper left half-plane, see Figure 7. Assume now that ∂S can be represented by the RLC of the method (cf.Remark 2.7).As we have seen, the RLC is the image of the function µ(•) in (10) for LMMs, or the union of the images of the functions µ 1,2 (•) in (13) for second-derivative multistep methods.The function µ is given as a simple ratio, but to get the explicit forms of µ 1,2 , one should solve a quadratic equation.As the value of k gets larger, these explicit formulae for µ 1,2 corresponding to a k-step second-derivative multistep method become more and more complicated.Moreover, obtaining explicit and practically useful parametrized formulae for the RLCs associated with multistep methods based on higher-than-second order derivatives would almost be impossible.
To avoid these difficulties, we now describe a more general and effective technique which reduces the determination of the stability angles to the solution of a suitable system of we have e iϑ = (i − t)/(i + t); so instead of solving Φ(e iϑ , µ) = 0 for µ, we can solve without trigonometric functions.Notice that originally we have ϑ ∈ [0, 2π] in e iϑ , or equivalently, ϑ ∈ (−π, π], but π is not in the range of the function 2 arctan, therefore we define to restore the missing µ value(s) due to the reparametrization.Then, clearly, (15) can be brought to the form Q(t, µ)/R(t) = 0 with some (complex) polynomials Q and R. By writing µ = a + bi (a, b ∈ R) we get that there exist two real polynomials Q re : R 3 → R and Q im : R 3 → R such that the solutions of Q(t, µ) = 0 for any fixed t ∈ R are obtained as the solutions of the system Now we eliminate t by taking the resultant [14] of Q re and Q im with respect to this parameter, and get that there exists a real polynomial F : R 2 → R such that if (16) By taking again the resultant of the first two polynomial equations, one of the variables, say b 0 , is eliminated.The resulting univariate polynomial yields in the general case finitely many possible a 0 values to choose from.With α denoting the angle (in radians) between L and the negative half of the real axis, we get that tan(α) = −b 0 /a 0 .To select the appropriate solution (a 0 , b 0 ) (and hence the appropriate tangent line L), we verify in the concrete case that (a 0 , b 0 ) ∈ ∂S ⊂ C = R 2 , and determine whether ∂S lies on the right-hand side of L.
The appropriately chosen α angle then yields the desired stability angle.

Results for the BDF methods
The simplest non-trivial case illustrating the steps in Section 3.1 is the determination of the stability angle for the 3-step BDF method.Formula (11) with k = 3 yields the following trigonometric parametrization of the RLC in R 2 after a simplification: After eliminating the trigonometric functions, (15) can be written as Then Q re and Q im in (16) become We eliminate t from this system and obtain Remark 3.1.The above RLC for the 3-step BDF method can also be parametrized as Here M −1 = {(20/3, 0)} ⊂ R 2 , corresponding to the t → ±∞ limiting value of the parametrization.
The remaining stability angle values for 4 ≤ k ≤ 6 can be computed analogously, so Table 1 shows only the final exact results.9) has a double root with modulus equal to 1. Therefore µ ± ∈ ∂S \ S, hence this S is not closed.The right figure depicts the 6 roots of P 1 (•, µ + ), and the double root is located at 1  2 1 + i √ 3 (note that µ + ∈ C \ R, so these roots are not symmetric with respect to the real axis).

Results for the Enright methods
By applying the algorithm described in Section 3.1, we can exactly determine the stability angles for the Enright methods, see Table 2.But since the c E k values are much more complicated algebraic numbers than the corresponding c B k constants in Table 1, Table 2 contains only a numerical approximation to the exact stability angles.Remark 3.5.Besides the stability angle, there are other measures of stability for A(α)stable methods.One of these characteristics is the stiff stability abscissa, being the smallest constant D > 0 such that {z ∈ C : Re(z) ≤ −D} ⊂ S. For example, for the 3-step Enright method, Table 3.1 in [17, Chapter V.3] contains the approximate value D ≈ 0.103.By using our implicit representation of ∂S, it is straightforward to determine the exact value of D ≈ 0.10341810907195; it is an algebraic number of degree 12, and the total number of digits in the coefficients of its defining integer polynomial is 529.

It turns out that c
As for the k = 4 case, the algebraic degree of c Enr can be given as roots of increasingly more involved integer polynomials, so we do not reproduce these polynomials here.During the computations in the k = 7 case, for example, we had to manipulate intermediate polynomials of degree of a few hundred, or polynomials with a total number of coefficient digits of approximately 470000.We could describe the final defining polynomial for c Enr  The angle between the dashed red line and the negative half of the real axis has also been computed exactly; its approximate value is ≈ 89.9999527 • .

Optimal disk inclusions
As for the largest inscribed disk in the stability region S, we again expect-similarly to Section 3.1-that ∂S (or the RLC) and the optimal disk possess a common tangent line (with point of tangency different from the origin).By using • the implicit algebraic form F (a, b) = 0 of the RLC, • the implicit equation (a + r) 2 + b 2 − r 2 = 0 for the boundary of the inscribed disk, • and the condition for a common tangent line , we obtain a system of 3 polynomial equations in 3 unknowns (a, b, r).By taking resultants and successively eliminating the variables (a, b), we obtain a univariate polynomial in r whose positive root will yield the optimum value of the stability radius.The exact optimal stability radii r BDF k for the k-step BDF methods (3 ≤ k ≤ 6) are found in Table 3; see Figure 8 also.The degree of the algebraic number r BDF k is 2, 3, 5, 5 for 3 ≤ k ≤ 6, respectively.3.
Remark 4.1.It is quite surprising that the algebraic numbers listed in Table 3 have such a low degree for the following reasons.For the 3-step BDF method, the univariate polynomial r mentioned above has degree 28, but it can be split into several factors of lower degree, and has a unique positive root r BDF 3 ≈ 7.0497.For the 4-step BDF method, the corresponding rpolynomial has degree 52 and a unique positive root ≈ 2.7272.The r-polynomial for the 5-step BDF method has degree 88 and a unique positive root ≈ 1.3579.Finally, the r-polynomial for the 6-step BDF method has degree 128, and a unique positive root ≈ 0.5599.

Optimal stability angle in a family of multistep methods
In [18], ODEs of the form u (t) = F (u(t)) + G(u(t)), u(0) = u 0 are considered, with F and G representing non-stiff and stiff parts of the equation, respectively.To solve these equations numerically, the authors construct several implicit-explicit (IMEX) LMMs, and thoroughly analyze them from the viewpoint of numerical monotonicity, boundedness and stability.Their analysis involves finding optimal methods with respect to various criteria in certain families.
Here we take their simplest case study from [18, Section 3.2.1], a 2nd-order, 3-step explicit method augmented by an implicit method (note that we changed their notation from b j to β j ): The values of β 2 := −3β 0 − 2β 1 + 3 and β 3 := 2β 0 + β 1 − 3 2 are determined from the order conditions, so (18) becomes a 2-parameter family of methods, with real parameters β 1 and β 0 .The three figures in [18, Figure 1] then depict the A(α)-stability angles, the "damping factors" and the "absolute error constants", respectively, of members of the family (18).In what follows, we do not consider these last two categories but focus only on the leftmost figure in [18, Figure 1]-as the authors conclude in [18, Section 3.2.1], a method with large stability angle does not necessarily have a good damping factor or a small error constant, and vice versa; the different optimization criteria are often conflicting.In other words, our goal in this section is to find the IMEX method in the family (18) with the largest stability angle.
To begin the A(α)-stability investigation, the authors of [18] define the usual linear test functions F (u) := λu and G(u) := λu.They then assume that ∆t • λ = iη and ∆t • λ = ξ with η ∈ R and R ξ ≤ 0: this choice is relevant "for example, for advection-diffusion equations if central finite differences or spectral approximations are used in space".These assumptions lead to the following characteristic polynomial of the IMEX multistep family, see [18, (2.4)-(2.7)]: To create the leftmost figure in [18, Figure 1] approximately indicating the optimal stability angle within the family, the authors use (19) to construct the RLCs and study these curves "for ξ → −∞" to estimate the stability angles 1 .
In the rest of this section we confirm their numerical findings, but we solve the optimization problem rigorously and exactly.We have selected this family (18) because the final result-the optimal stability angle-has a particularly simple form (see our Theorem 5.3 below), and, at the same time, our straightforward approach based on the theorems cited in Section 2.2 is readily illustrated.We emphasize that our analysis avoids the construction of the RLCs: as we have seen (for example, in Figure 3), they may have complicated self-1 When the stability angle α of a method is defined in [18, Sections 2.3 and 3.2.1]notice that we should require that the sector ξ ≤ 0, |η/ξ| ≤ tan(α) with angle α ≤ π/2 be included in the stability region in the (ξ, η)-plane (with the ξ = 0 and α = π/2 cases interpreted appropriately).In other words, arctan(α) in [18] is to be replaced by tan(α), otherwise the sector would not "open wide enough" and A-stability would not be recovered in the α → π/2 − limit.See also Footnote 2.
intersections, and it is often not obvious a priori whether a particular segment of the RLC coincides with the stability region boundary or not.

Summary of the main steps and results
By rearranging (19) and inserting the values of β 2 and β 3 given below ( 18), we define where Our goal is to find the parameters (β 1 , β 0 ) such that the stability region contains the infinite sector with the largest m > 0 in the definition of A(α)-stability.In other words, we are to find holds with the largest possible m > 0. Note that for convenience we have identified C with R 2 , hence stability regions in this section are subsets of R 2 .As a first step, Lemma 5.1 below yields a necessary condition for the inclusion (22).In its proof-presented in Appendix A.1-we use the argument proposed in [18] and consider the ξ → −∞, η = 0 limiting values.At this point it is convenient to recall the notion of Then a method of the form (18) is not As a consequence, from now on we can assume (β 1 , β 0 ) ∈ W , see Figure 9.Note that the orientation of the axes in Figure 9 and in the leftmost figure in [18, Figure 1] is the same: the β 1 -axis is horizontal, while the β 0 -axis is vertical.Lemma 5.1 thus also proves that the wedge-like object in the parameter space in the leftmost figure in [18, Figure 1] is indeed a perfect (infinite) wedge given by W .

Optimal parabola inclusion in a family of multistep methods
In the previous section we demonstrated how one can find the optimal sector in a family of stability regions of multistep methods.Here we show that the same algebraic approach allows us to replace the sector with more general shapes: we use again the multistep family (18) as a test example and determine the optimal stability region that contains the largest parabola.
The motivation for considering the shape of a parabola comes from [18] ("for advectiondiffusion equations, stability within a parabola2 can be more relevant than for a wedge"), or from [5, Sections 3-4] (where linearly implicit Runge-Kutta methods are developed for the numerical integration of semidiscrete equations originating from spatial discretizations of PDEs of advection-reaction-diffusion type).With P β 1 ,β 0 and S β 1 ,β 0 defined in ( 20)-( 21), we are now looking for the largest possible m > 0 such that the stability region of a suitable member of the family (18) contains the parabola that is, the inclusion holds.Clearly, we need • A-stability again to have (25) with some m > 0, so from now on, by Lemma 5.1, we can assume that (β 1 , β 0 ) ∈ W (see Figure 9).
In Appendix B.1 we apply a simple geometric argument: we first formulate the RLCs for the members of the multistep family as implicit curves {(ξ, η) ∈ R 2 : F β 1 ,β 0 (ξ, η) = 0}, then invoke the notion of discriminant [14] to construct a polynomial in m (and depending on the parameters β 1 and β 0 ) whose suitable root can yield the optimal value m opt in (25).The simple observation is the same as the one used in Section 3.1 (or in Section 4): the optimal inscribed object (now a parabola) touches the boundary of the optimal stability region.
Based on this technique and by using Mathematica, we conjecture that the parameter values β 1 = 1/5 and β 0 = 37/40 give m opt = 6/5.In Appendix B.2 we use a uniqueness argument to rigorously prove this conjecture.We emphasize that, similarly to Appendix A.2, no RLCs are involved in this uniqueness proof; the RLCs are used only as auxiliary objects to conjecture the optimum.Given the complexity of intermediate calculations, it is again surprising that the final result m opt is a simple rational number.In summary, we have the following theorem.Theorem 6.1.Suppose that (β 1 , β 0 ) ∈ W . Then the largest m > 0 such that (25) holds is m ≡ m opt := 6/5.Remark 6.2.The authors of [18] observe that "for the methods considered in this paper, a large angle α will correspond to a large β" (with α and β interpreted in our Footnotes 1 and 2).According to our results, the optimal (β 1 , β 0 ) parameter pairs (3/8, 3/4) and (1/5, 37/40)-determining the stability regions with the largest inscribed sector and parabola, respectively-do not coincide, although they are both located on the right boundary of W in Figure 9 (see also Remark B.1).

A Appendix
A.1 The proof of Lemma 5.1 Clearly, if |ξ| is large enough, the coefficients of the RHS polynomial can be arbitrarily close to 0. So by the fact that the roots of a polynomial are continuous functions of its coefficients, we get that "the ζ j roots of LHS(ζ) = RHS(ζ) can be made arbitrarily close to those of LHS(ζ) = 0 by choosing |ξ| large".To make the previous "statement" precise, we distinguish two cases according to whether the leading coefficient of LHS vanishes or not: for β 0 = 0, the LHS polynomial has at most two roots, whereas the difference LHS − RHS has three.Case I: β 0 = 0.By the above statement we easily see that if LHS β 1 ,β 0 (•) / ∈ vN, then P β 1 ,β 0 (•, ξ, 0) / ∈ svN for |ξ| large enough.We now show that So let us suppose in the rest of Case I that (β 1 , β 0 ) / ∈ W and β 0 = 0. Case I a .First we check the case when cc and, since now 2β 1 − 3 = 0, we can apply Theorem 2.4 to the above polynomial in Case I b .The conditions cc LHS β 1 ,β 0 (•) = 0 = β 0 mean that we can apply Theorem 2.4 to LHS β 1 ,β 0 (•).It is easy to verify that (LHS β 1 ,β 0 (•)) r does not vanish identically, so LHS β 1 ,β 0 (•) ∈ vN if and only if We show in Cases I b1 and I b2 below that ( 27) never occurs.First we observe that the inequality constraint in (27) yields that lc (LHS Case I b2α : when (LHS β 1 ,β 0 (•)) rr ≡ 0 and [(LHS β 1 ,β 0 (•)) r ] ∈ vN.In this case, however, the unique root of the polynomial [(LHS β 1 ,β 0 (•)) r ] , , for which we again have |ζ 1 | > 1, completing Case I.
A.2 The proof of Theorem 5.3 Proof.In the proof we suppose m > 0 and, due to Lemma 5.1, that (β 1 , β 0 ) ∈ W .
Remark B.1.By setting (β 1 , β 0 ) = (3/8, 3/4) for example (corresponding to the "sectoroptimal" method in Section A.2), we have that the RLC in (36) is identical to 3/16 By studying the positive roots of ∆ β 1 ,β 0 (•) as (β 1 , β 0 ) is varied within W , we can conjecture that the value of m in (25) cannot be greater than 6/5 for the family (18).Moreover, m = 6/5 occurs only for β 1 = 1/5 and β 0 = 37/40, and in this case the RLC and the upper parabola branch touch each other at (ξ, η) = (−10/7, 2 3/7).Since the polynomial ∆ β 1 ,β 0 (m) is much more complicated than the corresponding polynomial Q β 1 ,β 0 (m) in Appendix A.2, this time Mathematica could not confirm in a reasonable amount of computing time that the value m = 6/5 is indeed the optimal one.B.2 The proof of optimality in Theorem 6.1 However, once the unique optimum has been conjectured properly, the proof of optimality becomes straightforward to complete.
But (P β 1 ,β 0 (•, ξ 0 , η 0 )) rr is a linear polynomial (it is easily checked that it cannot be a constant polynomial), so its unique (non-real complex) root can be directly expressed: one sees that the absolute value of this root is ≤ 1 if and only if The product of the first two factors is strictly negative in W , and a standard constrained optimization computation shows that the third factor is ≥ 0 in W if and only if (β 1 , β 0 ) = (1/5, 37/40), completing Step 1.
These mean that the stability region boundary in the optimal case coincides with the corresponding RLC (in the left half-plane), and the boundary of the optimal inscribed parabola touches the stability region boundary in the open upper left half-plane at exactly one point, see Figure 13  is also included here.

Figure 1 :
Figure 1: RLCs for the k-step BDF methods for 1 ≤ k ≤ 6.The stability region of the method in each case is the unbounded component of C.

Figure 2 :
Figure2: RLC for the unstable 7-step BDF method in red (left), and a close-up near the origin (right).For comparison, the curves from Figure1are also superimposed as dashed gray curves.

Figure 3 :Figure 4 :
Figure3: The black curve in the left figure shows the boundary ∂S of the stability region of the (unstable) 7-step BDF method; ∂S is non-differentiable at one point.The stability region is the unbounded outer component.The red curve segment near the origin is not part of ∂S, it is a subset only of the RLC as displayed in Figure2.The small brown rectangle in the center is shown in detail in the right figure.The red curve in the right figure is again the RLC.The 6 black dots depict the set of µ values such that P 1 (•, µ) in (9) has multiple roots (there are no other µ ∈ C parameters with this property for k = 7).The polynomial P 1 (•, µ) has 1, 2 and 3 roots outside the unit disk for µ values in the dark brown, light brown and orange regions, respectively; P 1 (•, µ) cannot have 4 or more roots outside the unit disk.Each of the three self-intersections of the RLC in this figure (as well as the self-intersection of the RLC seen only in the left figure) corresponds to a µ value for which P 1 (•, µ) has two distinct roots with modulus 1. Exactly computing, for example, the unique value of µ † ≈ −2.68886 • 10 −6 + 0.275988i in the open upper half-plane where the RLC crosses itself was a non-trivial task: it took Mathematica 86 minutes to explicitly determine the coefficients of the integer polynomial defining µ † and having degree 30.The RLCs for the k-step BDF methods with 1 ≤ k ≤ 6 do not have any self-intersections; other singularities may occur, see Figure6.

Figure 5 :
Figure 5: RLCs for the unstable 8-step Enright method in red.The stability region S is not connected, C \ S is the annulus-like region.For comparison, the curves from Figure 4 are displayed as dashed gray curves.
holds for some t ∈ R, then F (a, b) = 0 should hold with some a, b ∈ R. Hence, after identifying C with R 2 , we see that the RLC can be represented as the implicit algebraic curve C ∪ M −1 withC := {(a, b) ∈ R 2 : F (a, b) = 0}.Assuming that the set M −1 is finite (it has at most two elements in the case of the BDF and Enright methods we are interested in), we ignore this component and focus only on C. Suppose now that a line L passes through the origin and touches C in the open upper left half-plane at some (a 0 , b 0 ) with a 0 < 0 < b 0 .By assuming that C can be represented locally as the graph of an implicit function near (a 0 , b 0 ) ∈ C, we easily get, by differentiating a → F (a, b(a)), that (a 0 , b 0 ) satisfies (a, b) := 432 108a 6 − 1188a 5 + 9a 4 36b 2 + 439 − 2a 3 1188b 2 + 3121 + 9a 2 36b 4 + 394b 2 + 547 − 54a 22b 4 + 17b 2 + 30 + 27b 4 4b 2 − 15 .Now b is eliminated from the first two equations of (17), and we get that the possible choices for a 0 are the negative real roots of a 4 (24a − 25) 4 (5324a + 405) 2 6a 2 − 13a + 9 2 = 0, yielding the unique value a 0 = −405/5324.Substituting this a 0 into (17) we get the unique value b 0 = 987 √ 35/5324, hence tan(α) = −b 0 /a 0 = (329 7/5)/27 is the only possible value for the tangent of the stability angle.Finally, we verify that the corresponding tangent line L passing through the origin has no other intersection point with ∂S in the open upper left quadrant, and ∂S lies on the right side of L.

Enr 3 is
an algebraic number of degree 22, being the unique positive root of the following even polynomial with coefficients {6621625501626720011970719022734459520000000000000000,

7 by≈Remark 3 . 6 .
175000 characters in Mathematica.Let us consider the Enright stability region corresponding to k = 7.As we already remarked earlier, there are exactly two lines that pass through the origin and are locally tangent to the boundary curve at some point in the open upper left half-plane, see Figure7.Within the BDF family for 1 ≤ k ≤ 6 or in the Enright family for 1 ≤ k ≤ 7, this phenomenon occurs only in the present case.

Figure 7 :
Figure 7: Part of the boundary of the stability region of the 7-step Enright method near the origin (solid black curve) together with the two (dashed red and black) lines that pass through the origin and are locally tangent to the boundary curve at some point in the open upper left half-plane.Due to the scaling, the dashed black line is seen only in the larger plot window on the right.The stability angle α E 7 ≈ 37.6 • of the method is determined by the dashed black line; the red line has additional intersection points with the boundary curve.The angle between the dashed red line and the negative half of the real axis has also been computed exactly; its approximate value is ≈ 89.9999527 • .

Figure 8 :
Figure 8: The largest inscribed disk |z + r| ≤ r (with red boundary) in the stability region of the k-step BDF method for k = 4 (left) and k = 6 (right), see Table3.

Figure 12 Figure 12 :
Figure 12: The implicit curve {(ξ, η) ∈ R 2 : F opt (ξ, η) = 0} (see Remark A.1), being the boundary of the optimal stability region S 3/8,3/4 in the left half-plane, is shown in the left figure and a close-up in the right figure in black.The dashed red lines represent the boundary of the largest infinite sector A 1/2 that can be included in the stability region.
, Chapter V.3, p. 276, Exercise 2] the stiff instability of the Enright formulas for k ≥ 8 is still mentioned as an open problem.B The stability radius of a multistep method is the largest number r > 0 such that the inclusion {z ∈ C : |z + r| ≤ r} ⊂ S holds.The stability radius plays an important role when analyzing boundedness properties of multistep methods.For example, it has been proved [42, Theorem 3.1] that this radius is the largest step-size coefficient for linear boundedness of a LMM satisfying some natural assumptions.
The class Sch is referred to as strongly stable polynomials in[4, p. 345].The property Q ∈ svN is often expressed by saying that Q satisfies the root condition.

Table 2 :
Stability angles α Enr Remark 3.4.By rounding the values of α Enr k given in Table 2 to two decimal places, we recover the approximate values of these stability angles in [17, Chapter V.3, Table 3.1].