Optimal subsets in the stability regions of multistep methods

TL;DR

This paper develops computational techniques to precisely analyze and optimize the stability regions of multistep numerical methods for differential equations, including exact calculations of stability angles and radii.

Contribution

It introduces straightforward methods for exact stability region analysis and optimization within parametric families of multistep methods using algebraic and recursive techniques.

Findings

Exact stability angles for BDF and second-derivative methods computed.

Largest stability radius determined for BDF methods.

Optimal methods identified within parametric families for stability angle and region.

Abstract

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.

Equations

$$\{z\in \mathbb{C}\setminus \{0\}:|\arg (-z)|\le \alpha \}\subset \mathcal{S},$$

$$\{z\in \mathbb{C}:|z+r|\le r\}\subset \mathcal{S}$$

$$Q^{\mathbf{r}}(z):=\frac{\overline{a_n}\cdot \left({\sum }_{j=0}^n{a}_j{z}^j\right)-a_0\cdot \left({\sum }_{j=0}^n\overline{a_{n-j}}{z}^j\right)}{z}=$$

$$\sum _{j=1}^n\left(\overline{a_n}\cdot a_j-a_0\cdot \overline{a_{n-j}}\right)z^{j-1},$$

$$y^{\prime }(t)=f(t,y(t)),y(t_0)=y_0$$

$$\{\begin{array}{rl}& \sum _{j=0}^s\sum _{\ell =0}^k{a}_{j,\ell }{\mu }^j{y}_{n+\ell }=0,n\in \mathbb{N},\\ & a_{j,\ell }\in \mathbb{R},\sum _{j=0}^s|a_{j,k}|>0,\mu :=h\lambda .\end{array}$$

$$\Phi (\zeta ,\mu ):=\sum _{j=0}^s\sum _{\ell =0}^k{a}_{j,\ell }{\mu }^j{\zeta }^{\ell }(\zeta \in \mathbb{C}).$$

$$\mathcal{S}:=\{\mu \in \mathbb{C}:\text{the degree of }\Phi (\cdot ,\mu )\text{ is exactly }k,\text{ and }\Phi (\cdot ,\mu )\in \mathbf{svN}\}.$$

$$\{\mu \in \mathbb{C}:\Phi (\cdot ,\mu )\in \mathbf{Sch}\},$$

$$\{\mu \in \mathbb{C}:\text{all roots }{\zeta }_j(\mu )\text{ of }\zeta \mapsto \Phi (\zeta ,\mu )\text{ satisfy }|{\zeta }_j(\mu )|\le 1,$$

$$\text{ and multiple roots satisfy }|{\zeta }_j(\mu )|<1\},$$

$$\{\mu \in \overline{\mathbb{C}}:\text{roots }{\zeta }_j\text{ of }\Phi (\zeta ,\mu )=0\text{ satisfy }|{\zeta }_j(\mu )|\le 1,$$

$$\text{ and if }|{\zeta }_j|=1,\text{then it is a simple root}\},$$

$$\mathcal{E}:=\{\mu \in \mathbb{C}:C_k(\mu )=0\}.$$

$$\Gamma :=\{\mu \in \overline{\mathbb{C}}:\exists \zeta \in \mathbb{C}\text{ with }|\zeta |=1\text{ and }\Phi (\zeta ,\mu )=0\}.$$

$$\sum _{j=0}^k({\alpha }_j{y}_{n+j}-h{\beta }_j{f}_{n+j})=0,$$

$$\varrho (\zeta ):=\sum _{j=0}^k{\alpha }_j{\zeta }^j\text{and}\sigma (\zeta ):=\sum _{j=0}^k{\beta }_j{\zeta }^j,$$

$$\Phi (\zeta ,\mu )\equiv P_1(\zeta ,\mu ):=\varrho (\zeta )-\mu \sigma (\zeta ).$$

$$[0,2\pi ]\ni \vartheta \mapsto \mu (\vartheta ):=\frac{\varrho \left(e^{i\vartheta }\right)}{\sigma \left(e^{i\vartheta }\right)}.$$

$$\sum _{j=1}^k\frac{1}{j}{\nabla }^j{y}_{n+1}=h{f}_{n+1},$$

$$\mu (\vartheta )\equiv \sum _{j=1}^k\frac{1}{j}(1-e^{-i\vartheta }{)}^j.$$

$$\sum _{j=0}^k({\alpha }_j{y}_{n+j}-h{\beta }_j{f}_{n+j}-h^2{\gamma }_j{g}_{n+j})=0,$$

$$\Phi (\zeta ,\mu )\equiv P_2(\zeta ,\mu ):=\sum _{j=0}^k({\alpha }_j-\mu {\beta }_j-{\mu }^2{\gamma }_j){\zeta }^j.$$

$$[0,2\pi ]\ni \vartheta \mapsto {\mu }_{1,2}(\vartheta ),$$

$$y_{n+1}=y_n+h{f}_{n+1}-h\sum _{j=1}^k\left(\frac{1}{j}\left(\sum _{\ell =j}^k{\nu }_{\ell }\right){\nabla }^j{f}_{n+1}\right)+h^2\left(\sum _{\ell =0}^k{\nu }_{\ell }\right)g_{n+1},$$

$${\nu }_{\ell }:=(-1{)}^{\ell }{\int }_0^1(\tau -1)(\genfrac{}{}{0ex}{}{1-\tau }{\ell })d\tau (0\le \ell \le k)$$

$$\vartheta =2\arctan (t)(t\in \mathbb{R}),$$

$$\Phi \left(\frac{i-t}{i+t},\mu \right)=0$$

$$M_{-1}:=\{\mu \in \mathbb{C}:\Phi \left(e^{i\pi },\mu \right)=0\}$$

$$\{\begin{array}{rl}{Q}_{\text{re}}(t,a,b)& =0\\ Q_{\text{im}}(t,a,b)& =0.\end{array}$$

$$\{\begin{array}{rl}F(a_0,b_0)& =0\\ a_0\cdot {\partial }_{1}F(a_0,b_0)+b_0\cdot {\partial }_{2}F(a_0,b_0)& =0\\ a_0& <0\\ b_0& >0.\end{array}$$

Full text

First we describe a simple procedure to exactly calculate the stability angle $\alpha$ in the definition of $A(\alpha)$-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\leq k\leq 6$) and discover that the values of $\tan(\alpha)$ are surprisingly simple algebraic numbers of degree 2, 2, 4 and 2, respectively. In contrast, the corresponding values of $\tan(\alpha)$ for the $k$-step second-derivative multistep methods of Enright ($3\leq k\leq 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\leq k\leq 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.

1 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 ${\mathcal{S}}\subset{\mathbb{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 ${\mathcal{S}}$, or

• b)

a family of stability regions ${\mathcal{S}}_{\beta}$ parametrized by some $\beta\in\mathbb{R}^{d}$,

and a family of subsets of ${\mathbb{C}}$, denoted by ${\mathfrak{F}}$. Due to their relevance in applications, we will consider the following three classes:

• •

${\mathfrak{F}}={\mathfrak{F}}^{\text{\,sect}}_{\alpha}$ is the family of infinite sectors in the left half of ${\mathbb{C}}$, with vertex at the origin, symmetric about the negative real axis, and parametrized by the sector angle $\alpha\in(0,\pi/2)$;

• •

${\mathfrak{F}}={\mathfrak{F}}^{\text{\,disk}}_{r}$ is the family of disks in the left half of ${\mathbb{C}}$, symmetric with respect to the real axis, touching the imaginary axis, and parametrized by the disk radius $r>0$;

• •

${\mathfrak{F}}={\mathfrak{F}}^{\text{\,para}}_{m}$ is the family of parabolas in the left half of ${\mathbb{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\in{\mathfrak{F}}$ with the largest parameter ($\alpha$, $r$, or $m$) such that

• •

$H\subset{\mathcal{S}}$ in case a);

• •

$H\subset{\mathcal{S}}_{\beta_{\text{opt}}}$ for some stability region in the family in case b), but $H\not\subset{\mathcal{S}}_{\beta}$ for $\beta\neq\beta_{\text{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.

1.1 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\in{\mathbb{C}}:\mathrm{Re}(z)\leq 0\}$ belongs to ${\mathcal{S}}$. Many useful methods are not $A$-stable, still, ${\mathcal{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(\alpha)$-stability: a method is $A(\alpha)$-stable with some $0<\alpha<\pi/2$ if

[TABLE]

where the argument of a non-zero complex number satisfies $-\pi<\mathrm{arg}\leq\pi$. The largest $0<\alpha<\pi/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, $\overset{\circ}{\text{A}}$-stability, stiff stability, or asymptotic $A(\alpha)$-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(\alpha)$-stability for a given $\alpha$ 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(\alpha)$-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(\alpha)$-stable for a given $\alpha$; 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 $\alpha_{k}^{\text{BDF}}$ denoting the stability angle of the $k$-step BDF method for $3\leq k\leq 6$, we show that $\tan\left(\alpha_{k}^{\text{BDF}}\right)$ is an unexpectedly simple algebraic number, having degree 2 for $k\in\{3,4,6\}$, and degree 4 for $k=5$; see Table 1. For the $k$-step Enright methods with $3\leq k\leq 7$, the corresponding constants $\tan\left(\alpha_{k}^{\text{Enr}}\right)$ (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 $\alpha\in(0,\pi/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\in\{1,2\}$ are $A$-stable. For $k\geq 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\leq k\leq 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\in\{1,2\}$, see [17], and unstable for $k\geq 8$. More precisely, [11] proves that these methods are not $A_{0}$-stable for $k\geq 8$, hence they cannot be stiffly stable either, see [23, Theorem 3] (cf. [22, 26]). However, in [17, Chapter V.3, p. 276, Exercise 2] the stiff instability of the Enright formulas for $k\geq 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

[TABLE]

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.

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\leq k\leq 6$.

** C ** The RLC, as the graph of a $[0,2\pi]\to\mathbb{C}$ function (or a union of such functions for generalized LMMs), yields information about the boundary of the stability region, $\partial{\mathcal{S}}$. It is known, however, that in general the RLC does not coincide with $\partial{\mathcal{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 ${\mathbb{C}}$ determined by the RLC to see which component belongs to ${\mathcal{S}}$ and which one to ${\mathbb{C}}\setminus{\mathcal{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 ${\mathcal{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 non-trivial 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].

1.2 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.

2 Preliminaries

2.1 Notation

The set of natural numbers $\{0,1,\ldots\}$ is denoted by $\mathbb{N}$. For $z\in\mathbb{C}$, $\mathrm{Re}(z)$, $\mathrm{Im}(z)$, and $\overline{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\subset\mathbb{C}$ is $\partial H\subset\mathbb{C}$. When describing certain algebraic numbers of higher degree, a polynomial $\sum_{j=0}^{n}a_{j}x^{j}$ with $a_{j}\in\mathbb{Z}$, $a_{n}\neq 0$ and $n\geq 3$ will be represented simply by its coefficient list $\{a_{n},a_{n-1},\ldots,a_{0}\}$. For a polynomial $Q(z)=\sum_{j=0}^{n}a_{j}z^{j}$ with $0\leq n\in\mathbb{N}$, $a_{j}\in\mathbb{C}$ ($0\leq j\leq n$) and $a_{n}\neq 0$, we denote its degree, leading coefficient and constant coefficient by $\deg Q=n$, $\mathfrak{lc}\,Q=a_{n}$ and $\mathfrak{cc}\,Q=a_{0}$. The acronyms RLC and LMM stand for root locus curve and linear multistep method, respectively.

2.2 The Schur–Cohn reduction

In the rest of this section we assume that $Q$ is a univariate polynomial with $\deg Q\geq 1$, and follow the terminology of [37]—we have explicitly added the $\deg Q\geq 1$ condition, being implicit in [37]. We say that

• •

$Q$ is a Schur polynomial, $Q\in\mathbf{Sch}$, if its roots lie in the open unit disk;

• •

$Q$ is a von Neumann polynomial, $Q\in\mathbf{vN}$, if its roots lie in the closed unit disk;

• •

$Q$ is a simple von Neumann polynomial, $Q\in\mathbf{svN}$, if $Q\in\mathbf{vN}$ and roots with modulus 1 are simple.

Remark 2.1**.**

The class $\mathbf{Sch}$ is referred to as strongly stable polynomials in [4, p. 345].

Remark 2.2**.**

The property $Q\in\mathbf{svN}$ is often expressed by saying that $Q$ satisfies the root condition.

The reduced polynomial of $Q(z)=\sum_{j=0}^{n}a_{j}z^{j}$ is defined as

[TABLE]

[TABLE]

so we have $\deg Q^{\mathbf{r}}\leq(\deg Q)-1$. When this reduction process is iterated, we write $Q^{\mathbf{r}\mathbf{r}}$ for $\left(Q^{\mathbf{r}}\right)^{\mathbf{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 $\mathfrak{lc}\,Q\neq 0\neq\mathfrak{cc}\,Q$ and $\deg Q\geq 2$.

Theorem 2.3**.**

$Q\in\mathbf{Sch}\Leftrightarrow(|\mathfrak{lc}\,Q|>|\mathfrak{cc}\,Q|\text{ and }Q^{\mathbf{r}}\in\mathbf{Sch})$.

Theorem 2.4**.**

$Q\in\mathbf{vN}\Leftrightarrow\text{ either }(|\mathfrak{lc}\,Q|>|\mathfrak{cc}\,Q|\text{ and }Q^{\mathbf{r}}\in\mathbf{vN})\text{ or }(Q^{\mathbf{r}}\equiv 0\text{ and }Q^{\prime}\in\mathbf{vN})$.

Theorem 2.5**.**

$Q\in\mathbf{svN}\Leftrightarrow\text{ either }(|\mathfrak{lc}\,Q|>|\mathfrak{cc}\,Q|\text{ and }Q^{\mathbf{r}}\in\mathbf{svN})\text{ or }(Q^{\mathbf{r}}\equiv 0\text{ and }Q^{\prime}\in\mathbf{Sch})$.

Remark 2.6**.**

Let us consider the following example when applying the theorems above, e.g. Theorem 2.4. For any $\lambda>0$, we set $Q_{\lambda}(z):=z^{2}+\lambda iz+1$. Then the roots of $Q_{\lambda}$ satisfy $|z_{1}(\lambda)|<1<|z_{2}(\lambda)|$, so $Q_{\lambda}\notin\mathbf{vN}$, and $Q_{\lambda}^{\mathbf{r}}=2\lambda i$. 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^{\mathbf{r}}\in\mathbf{vN}$ is undefined. In these cases, when $Q^{\mathbf{r}}$ is a non-zero constant, notice that neither $|Q^{\mathbf{r}}|<1$, nor $|Q^{\mathbf{r}}|=1$, nor $|Q^{\mathbf{r}}|>1$ can help us in general to determine whether $Q\in\mathbf{vN}$ or not (of course, the other condition $|\mathfrak{lc}\,Q|>|\mathfrak{cc}\,Q|$ is violated now); cf. the sentence above [37, Theorem 5.1].

2.3 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

[TABLE]

can be analyzed by studying the stability region of the method. When an $s$-stage $k$-step method ($s\geq 1$, $k\geq 1$ fixed positive integers; for $k=1$ we have a one-step method, while for $k\geq 2$ a multistep method) with constant step size $h>0$ is applied to the linear test equation $y^{\prime}=\lambda y$ ($\lambda\in{\mathbb{C}}$ fixed, $y(0)=y_{0}$ given), the method yields a numerical solution $(y_{n})_{n\in\mathbb{N}}$ that approximates the exact solution $y$ at time $t_{n}:=t_{0}+nh$ and satisfies a recurrence relation of the form [27]

[TABLE]

The characteristic polynomial associated with the method takes the form

[TABLE]

With $\Phi(\cdot,\mu)$ abbreviating the polynomial $\zeta\mapsto\Phi(\zeta,\mu)$, the stability region of the method is defined as

[TABLE]

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

[TABLE]

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

[TABLE]

that is, essentially, $\Phi(\cdot,\mu)\in\mathbf{svN}$. In [27, Formula (2.5)] the stability region is given by

[TABLE]

with $\overline{\mathbb{C}}$ denoting the extended complex plane.

We can regroup the terms in (4) as $\Phi(\zeta,\mu)=\sum_{\ell=0}^{k}C_{\ell}(\mu)\zeta^{\ell}$ with some suitable polynomials $C_{\ell}$. The inequality condition in (3) implies that the leading coefficient $C_{k}$ does not vanish identically; it may happen that for some exceptional $\mu$ values the leading coefficient is zero:

[TABLE]

For example, for the implicit Euler (IE) method $\Phi(\zeta,\mu)=\varrho(\zeta)-\mu\sigma(\zeta)=(1-\mu)\zeta-1$ with $\varrho(\zeta):=\zeta-1$ and $\sigma(\zeta):=\zeta$, so ${\mathcal{E}}=\{1\}$. For the $2$-step BDF method (BDF2), $\Phi(\zeta,\mu)=(3-2\mu)\zeta^{2}-4\zeta+1$, hence ${\mathcal{E}}=\{3/2\}$. If definition (6) (or (7)) is interpreted formally, we have for the IE method that ${\mathcal{E}}=\{1\}\subset{\mathcal{S}}$ (because (6) is satisfied vacuously). Similarly, for the BDF2 method, ${\mathcal{E}}=\{3/2\}\subset{\mathcal{S}}$ (because then the unique root of $\Phi(\zeta,3/2)=0$ is $\zeta=1/4$).

*However, elements of ${\mathcal{E}}$ or ${\mathcal{E}}\cap{\mathcal{S}}$ can be problematic.

(i) For $\mu\in{\mathcal{E}}$, the order of the recursion (3) decreases, thus, in general, the starting values $y_{0},y_{1},\ldots,y_{k-1}$ of the numerical method cannot be chosen arbitrarily.

(ii) Some exceptional values $\mu\in{\mathcal{E}}\cap{\mathcal{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 $\mu\in{\mathcal{E}}\cap{\mathcal{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 $\partial{\mathcal{S}}$ of the stability region (see Sections 2.4–2.5 below). In [27, Definition (2.21)], the RLC is given by*

[TABLE]

It can happen that $\partial{\mathcal{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 $\partial{\mathcal{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 $\varrho$ and $\sigma$ have no common root and $\varrho/\sigma$ is univalent on the set $\{z\in\overline{\mathbb{C}}:|z-1|>1\}$, so $Q(\mu)=1/(1-\mu)$ has Property C. Thus for the IE method $\partial{\mathcal{S}}=\Gamma$. As we have seen above, $1\in{\mathcal{E}}\cap{\mathcal{S}}$, so $1\in\partial{\mathcal{S}}$. On the other hand, $\Phi(\zeta,1)=\varrho(\zeta)-\sigma(\zeta)=-1$, so $1\notin\Gamma=\partial{\mathcal{S}}$. This apparent contradiction seems to indicate that the authors of [27] interpreted definition (7) intuitively: a root $\zeta=\infty$ is tacitly introduced as soon as the leading coefficient $C_{k}(\mu)$ becomes zero. So [27, Corollary 2.6], for example, actually relies on definition (5) rather than on definition (7) (or (6)).

The problem of vanishing leading coefficient is implicitly avoided in [33, p. 66], or in [40], because they impose a requirement on “all the roots $r_{s}$ ($s=1,\ldots,k$)”. Definition (5) above with the non-vanishing leading coefficient essentially appears, for example, in [41, Section 2.1] (where it is formulated for LMMs, that is, for $s=1$ in (3)), or in [42, Section 2].

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 $\partial{\mathcal{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.

2.4 The RLC of a LMM

A linear multistep method for (2) has the form

[TABLE]

where $f_{m}:=f(t_{m},y_{m})$, and the numbers $\alpha_{j}\in\mathbb{R}$ and $\beta_{j}\in\mathbb{R}$ ($j=0,\ldots,k$) are the suitably chosen method coefficients with $\alpha_{k}\neq 0$. The method is implicit, if $\beta_{k}\neq 0$. By setting

[TABLE]

the associated characteristic polynomial (4) becomes

[TABLE]

One way to study the stability region (5), or its boundary $\partial{\mathcal{S}}$ in the complex plane is to depict the RLC corresponding to the method [17]: observe that $P_{1}$ is linear in $\mu$, so $P_{1}(\zeta,\mu)=0$ implies $\mu=\varrho(\zeta)/\sigma(\zeta)$ (for $\sigma(\zeta)\neq 0)$. The RLC is then the image of the parametric curve

[TABLE]

2.4.1 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

[TABLE]

where $\nabla$ denotes the backward difference operator $\nabla y_{n+1}:=y_{n+1}-y_{n}$, and $\nabla^{j}y_{n+1}:=$ $\nabla^{j-1}y_{n+1}-$ $\nabla^{j-1}y_{n}$ (for $j>1$). It is known [17] that the corresponding RLC is

[TABLE]

Figures 1–3 show the RLCs for some BDF methods.

2.5 The RLC of a multiderivative multistep method

A second-derivative multistep method is more general than (8) and can be written as

[TABLE]

where $g_{n}:=g(t_{n},y_{n})$ with $g(t,y):=\partial_{1}f(t,y)+\partial_{2}f(t,y)\cdot f(t,y)$, and the method is determined by the coefficients $\alpha_{j}$, $\beta_{j}$ and $\gamma_{j}$, see [17]. Now the associated characteristic polynomial (4) becomes

[TABLE]

This time we have two RLCs:

[TABLE]

where $\mu_{1,2}$ are the two solutions of $P_{2}\left(e^{i\vartheta},\mu\right)=0$. For any choice of the method coefficients $\alpha_{j}$, $\beta_{j}$ and $\gamma_{j}$, one can construct $\mu_{1,2}$ explicitly, since $P_{2}$ is only quadratic in $\mu$.

2.5.1 RLCs for the Enright methods

The Enright methods are special cases of (12), and for $k\geq 1$ they are defined [17] as

[TABLE]

where

[TABLE]

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

3.1 The RLC in implicit algebraic form

Computing the stability angle of a method with stability region ${\mathcal{S}}$ is equivalent to finding the slope of the unique line $L$ that passes through the origin, touches $\partial{\mathcal{S}}$ at some point in the open upper left half-plane such that $\partial{\mathcal{S}}$ lies on the right-hand side of $L$ (viewed from the origin) in this quadrant. This last requirement is necessary since $\partial{\mathcal{S}}\cap L$ can consist of more points, even in the open upper left half-plane, see Figure 7.

Assume now that $\partial{\mathcal{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 $\mu(\cdot)$ in (10) for LMMs, or the union of the images of the functions $\mu_{1,2}(\cdot)$ in (13) for second-derivative multistep methods. The function $\mu$ is given as a simple ratio, but to get the explicit forms of $\mu_{1,2}$, one should solve a quadratic equation. As the value of $k$ gets larger, these explicit formulae for $\mu_{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 polynomial equations. Let us consider the equation $\Phi(e^{i\vartheta},\mu)=0$ (see (4)). By using the well-known Weierstrass substitution [43, pp. 382-383]

[TABLE]

we have $e^{i\vartheta}={(i-t)}/{(i+t)}$; so instead of solving $\Phi(e^{i\vartheta},\mu)=0$ for $\mu$, we can solve

[TABLE]

without trigonometric functions. Notice that originally we have $\vartheta\in[0,2\pi]$ in $e^{i\vartheta}$, or equivalently, $\vartheta\in(-\pi,\pi]$, but $\pi$ is not in the range of the function $2\arctan$, therefore we define

[TABLE]

to restore the missing $\mu$ value(s) due to the reparametrization. Then, clearly, (15) can be brought to the form $Q(t,\mu)/R(t)=0$ with some (complex) polynomials $Q$ and $R$. By writing $\mu=a+bi$ ($a,b\in\mathbb{R}$) we get that there exist two real polynomials $Q_{\text{re}}:\mathbb{R}^{3}\to\mathbb{R}$ and $Q_{\text{im}}:\mathbb{R}^{3}\to\mathbb{R}$ such that the solutions of $Q(t,\mu)=0$ for any fixed $t\in\mathbb{R}$ are obtained as the solutions of the system

[TABLE]

Now we eliminate $t$ by taking the resultant [14] of $Q_{\text{re}}$ and $Q_{\text{im}}$ with respect to this parameter, and get that there exists a real polynomial $F:\mathbb{R}^{2}\to\mathbb{R}$ such that if (16) holds for some $t\in\mathbb{R}$, then $F(a,b)=0$ should hold with some $a,b\in\mathbb{R}$. Hence, after identifying $\mathbb{C}$ with $\mathbb{R}^{2}$, we see that the RLC can be represented as the implicit algebraic curve $C\cup M_{-1}$ with $C:=\{(a,b)\in\mathbb{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})\in C$, we easily get, by differentiating $a\mapsto F(a,b(a))$, that $(a_{0},b_{0})$ satisfies

[TABLE]

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 $\alpha$ denoting the angle (in radians) between $L$ and the negative half of the real axis, we get that $\tan(\alpha)=-{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})\in\partial{\mathcal{S}}\subset\mathbb{C}=\mathbb{R}^{2}$, and determine whether $\partial{\mathcal{S}}$ lies on the right-hand side of $L$. The appropriately chosen $\alpha$ angle then yields the desired stability angle.

3.2 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 $\mathbb{R}^{2}$ after a simplification:

[TABLE]

[TABLE]

After eliminating the trigonometric functions, (15) can be written as

[TABLE]

Then $Q_{\text{re}}$ and $Q_{\text{im}}$ in (16) become

[TABLE]

We eliminate $t$ from this system and obtain

[TABLE]

[TABLE]

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

[TABLE]

yielding the unique value $a_{0}=-{405}/{5324}$. Substituting this $a_{0}$ into (17) we get the unique value $b_{0}={987\sqrt{35}}/{5324}$, hence $\tan(\alpha)=-b_{0}/a_{0}=({329\sqrt{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 $\partial{\mathcal{S}}$ in the open upper left quadrant, and $\partial{\mathcal{S}}$ lies on the right side of $L$.

Remark 3.1**.**

The above RLC for the $3$-step BDF method can also be parametrized as

[TABLE]

Here $M_{-1}=\{(20/3,0)\}\subset\mathbb{R}^{2}$, corresponding to the $t\to\pm\infty$ limiting value of the parametrization.

The remaining stability angle values for $4\leq k\leq 6$ can be computed analogously, so Table 1 shows only the final exact results.

Remark 3.2**.**

For $3\leq k\leq 6$, the BDF stability region includes an interval along the imaginary axis and containing the origin if and only if $k=5$ or $k=6$. For $k=5$ and $k=6$ the two intervals are

[TABLE]

and

[TABLE]

respectively.

Remark 3.3**.**

The boundary curve of the stability region of the $6$-step BDF method contains two cusp singularities, see Figure 6 (and compare with Figure 3). No other $\partial{\mathcal{S}}$ curve has this type of degeneracy in the BDF family for $1\leq k\leq 5$ or $k=7$. Since the cusp points for $k=6$ are not part of ${\mathcal{S}}$, the stability region in this case is not closed (nor open).

3.3 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_{k}^{E}$ values are much more complicated algebraic numbers than the corresponding $c_{k}^{B}$ constants in Table 1, Table 2 contains only a numerical approximation to the exact stability angles.

Remark 3.4**.**

By rounding the values of $\alpha_{k}^{\text{Enr}}$ 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].

It turns out that $c_{3}^{\text{Enr}}$ is an algebraic number of degree 22, being the unique positive root of the following even polynomial with coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Remark 3.5**.**

Besides the stability angle, there are other measures of stability for $A(\alpha)$-stable methods. One of these characteristics is the stiff stability abscissa, being the smallest constant $D>0$ such that $\{z\in\mathbb{C}:\mathrm{Re}(z)\leq-D\}\subset{\mathcal{S}}$. For example, for the $3$-step Enright method, Table 3.1 in [17, Chapter V.3] contains the approximate value $D\approx 0.103$. By using our implicit representation of $\partial{\mathcal{S}}$, it is straightforward to determine the exact value of $D\approx 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.

As for the $k=4$ case, the algebraic degree of $c_{4}^{\text{Enr}}$ is 28. The constants $c_{5}^{\text{Enr}}$, $c_{6}^{\text{Enr}}$ and $c_{7}^{\text{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_{7}^{\text{Enr}}$ by $\approx 175000$ characters in Mathematica.

Remark 3.6**.**

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 Figure 7. Within the BDF family for $1\leq k\leq 6$ or in the Enright family for $1\leq k\leq 7$, this phenomenon occurs only in the present case.

4 Optimal disk inclusions

As for the largest inscribed disk in the stability region ${\mathcal{S}}$, we again expect—similarly to Section 3.1—that $\partial{\mathcal{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

[TABLE]

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_{k}^{\text{\,BDF}}$ for the $k$-step BDF methods ($3\leq k\leq 6$) are found in Table 3; see Figure 8 also. The degree of the algebraic number $r_{k}^{\text{\,BDF}}$ is $2,3,5,5$ for $3\leq k\leq 6$, respectively.

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_{3}^{\text{\,BDF}}\approx 7.0497$. For the 4-step BDF method, the corresponding $r$-polynomial has degree 52 and a unique positive root $\approx 2.7272$. The $r$-polynomial for the 5-step BDF method has degree 88 and a unique positive root $\approx 1.3579$. Finally, the $r$-polynomial for the 6-step BDF method has degree 128, and a unique positive root $\approx 0.5599$.

5 Optimal stability angle in a family of multistep methods

In [18], ODEs of the form $u^{\prime}(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 $\beta_{j}$):

[TABLE]

The values of $\beta_{2}:=-3\beta_{0}-2\beta_{1}+3$ and $\beta_{3}:=2\beta_{0}+\beta_{1}-\frac{3}{2}$ are determined from the order conditions, so (18) becomes a 2-parameter family of methods, with real parameters $\beta_{1}$ and $\beta_{0}$. The three figures in [18, Figure 1] then depict the $A(\alpha)$-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(\alpha)$-stability investigation, the authors of [18] define the usual linear test functions $F(u):=\hat{\lambda}u$ and $G(u):=\lambda u$. They then assume that $\Delta t\cdot\hat{\lambda}=i\eta$ and $\Delta t\cdot\lambda=\xi$ with $\eta\in\mathbb{R}$ and $\mathbb{R}\ni\xi\leq 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)]:

[TABLE]

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 $\xi\to-\infty$” to estimate the stability angles111When the stability angle $\alpha$ of a method is defined in [18, Sections 2.3 and 3.2.1] notice that we should require that the sector

$\xi\leq 0,\quad|\eta/\xi|\leq\tan(\alpha)\quad\text{ with angle }\alpha\leq\pi/2$

be included in the stability region in the $(\xi,\eta)$-plane (with the $\xi=0$ and $\alpha=\pi/2$ cases interpreted appropriately). In other words, $\arctan(\alpha)$ in [18] is to be replaced by $\tan(\alpha)$, otherwise the sector would not “open wide enough” and $A$-stability would not be recovered in the $\alpha\to\pi/2^{-}$ limit. See also Footnote 2..

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-intersections, and it is often not obvious a priori whether a particular segment of the RLC coincides with the stability region boundary or not.

5.1 Summary of the main steps and results

By rearranging (19) and inserting the values of $\beta_{2}$ and $\beta_{3}$ given below (18), we define

[TABLE]

[TABLE]

where $\zeta\in\mathbb{C}$, $(\beta_{1},\beta_{0})\in\mathbb{R}^{2}$, $\xi\leq 0$ and $\eta\in\mathbb{R}$. Our goal is to find the parameters $(\beta_{1},\beta_{0})$ such that the stability region

[TABLE]

contains the infinite sector

[TABLE]

with the largest $m>0$ in the definition of $A(\alpha)$-stability. In other words, we are to find $(\beta_{1},\beta_{0})$ such that

[TABLE]

holds with the largest possible $m>0$. Note that for convenience we have identified $\mathbb{C}$ with $\mathbb{R}^{2}$, hence stability regions in this section are subsets of $\mathbb{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 $\xi\to-\infty$, $\eta=0$ limiting values. At this point it is convenient to recall the notion of $\overset{\circ}{\text{A}}$-stability [17, Chapter V.2]: a method is $\overset{\circ}{\text{A}}$-stable, if its stability region includes the non-positive reals $\{\xi\in\mathbb{R}:\xi\leq 0\}$. Clearly,

[TABLE]

Lemma 5.1**.**

Let us define

[TABLE]

Then a method of the form (18) is not $\overset{\circ}{\text{A}}$-stable for $(\beta_{1},\beta_{0})\notin W$.

As a consequence, from now on we can assume $(\beta_{1},\beta_{0})\in 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 $\beta_{1}$-axis is horizontal, while the $\beta_{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$.

Remark 5.2**.**

The assumption $(\beta_{1},\beta_{0})\in W$ implies $\beta_{0}>0$, so due to $\xi\leq 0$, the leading coefficient of (20), $1-\beta_{0}\xi$, cannot vanish (cf. Remark 2.7).

Then in Appendix A.2 we prove the main result of Section 5.

Theorem 5.3**.**

Suppose that $(\beta_{1},\beta_{0})\in W$. Then the largest $m>0$ such that (22) holds is $m\equiv m_{\text{opt}}:=1/2$.

In the proof we show that finding the optimal $(\beta_{1},\beta_{0})\in W$ is equivalent to finding the largest positive real root of a suitable polynomial in $m$ with coefficients depending on $\beta_{1}$ and $\beta_{0}$. We verify that this optimal root is located at $m_{\text{opt}}$, corresponding to the unique method with $(\beta_{1},\beta_{0})=W_{\text{opt}}:=(3/8,3/4)\in W$ and represented as a red dot in the parameter space in Figure 9. The black curve in the left half-plane in Figure 12 is the boundary of the optimal stability region, and the dashed red lines bound the largest inscribed infinite sector ${\mathcal{A}}_{1/2}$: the optimal stability angle satisfies $\tan(\alpha)=m_{\text{opt}}$. As a conclusion, the highest value in the scale adjacent to the leftmost figure in [18, Figure 1] should be exactly $\alpha=\arctan(1/2)\approx 0.463648$, that is, $\alpha\approx 26.5651^{\circ}$.

Remark 5.4**.**

Unlike in Section 6 (see Remark B.2), the boundary of the optimal sector ${\mathcal{A}}_{1/2}$ does not touch (or intersect) the boundary of the optimal stability region ${\mathcal{S}}_{3/8,3/4}$ in the open left half-plane.

Remark 5.5**.**

In [18, Section 3.2.1, (3.4)–(3.5)], the stability angles for two particular schemes from the family (18) are also approximated. For the IMEX-Shu(3,2) scheme

[TABLE]

[TABLE]

they obtain $\alpha_{\emph{Shu}}\approx 0.06$, and for the IMEX-SG(3,2) scheme

[TABLE]

they get $\alpha_{\emph{SG}}\approx 0.38$. Our technique easily yields the exact values

[TABLE]

and

[TABLE]

6 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 advection-diffusion equations, stability within a parabola222 Similarly to Footnote 1, an analogous typo is present in [18, Section 2.3] when the notion of “stability within a parabola” is defined. There we should have again $\tan$ instead of $\arctan$, that is,

$\xi\leq 0,\quad|\eta^{2}/\xi|\leq\tan(\beta)\quad\text{ with some angle }0<\beta\leq\pi/2.$

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_{\beta_{1},\beta_{0}}$ and ${\mathcal{S}}_{\beta_{1},\beta_{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

[TABLE]

that is, the inclusion

[TABLE]

holds. Clearly, we need $\overset{\circ}{\text{A}}$-stability again to have (25) with some $m>0$, so from now on, by Lemma 5.1, we can assume that $(\beta_{1},\beta_{0})\in 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 $\{(\xi,\eta)\in\mathbb{R}^{2}:F_{\beta_{1},\beta_{0}}(\xi,\eta)=0\}$, then invoke the notion of discriminant [14] to construct a polynomial in $m$ (and depending on the parameters $\beta_{1}$ and $\beta_{0}$) whose suitable root can yield the optimal value $\widetilde{m}_{\text{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 $\beta_{1}=1/5$ and $\beta_{0}=37/40$ give $\widetilde{m}_{\text{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 $\widetilde{m}_{\text{opt}}$ is a simple rational number. In summary, we have the following theorem.

Theorem 6.1**.**

Suppose that $(\beta_{1},\beta_{0})\in W$. Then the largest $m>0$ such that (25) holds is $m\equiv\widetilde{m}_{\text{opt}}:=6/5$.

Remark 6.2**.**

The authors of [18] observe that “for the methods considered in this paper, a large angle $\alpha$ will correspond to a large $\beta$” (with $\alpha$ and $\beta$ interpreted in our Footnotes 1 and 2). According to our results, the optimal $(\beta_{1},\beta_{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).

Appendix A Appendix

A.1 The proof of Lemma 5.1

Proof.

Let us fix some $(\beta_{1},\beta_{0})\in\mathbb{R}^{2}$. For $\xi<0$, $P_{\beta_{1},\beta_{0}}(\zeta,\xi,0)=0$ is equivalent to $\text{LHS}(\zeta)=\text{RHS}(\zeta)$ with

[TABLE]

and $\text{RHS}_{\xi}(\zeta):={(\zeta^{3}-{3\zeta^{2}}/{4}-{1}/{4})}/{\xi}$. Clearly, if $|\xi|$ is large enough, the coefficients of the RHS polynomial can be arbitrarily close to [math]. So by the fact that the roots of a polynomial are continuous functions of its coefficients, we get that “the $\zeta_{j}$ roots of $\text{LHS}(\zeta)=\text{RHS}(\zeta)$ can be made arbitrarily close to those of $\text{LHS}(\zeta)=0$ by choosing $|\xi|$ large”. To make the previous “statement” precise, we distinguish two cases according to whether the leading coefficient of LHS vanishes or not: for $\beta_{0}=0$, the LHS polynomial has at most two roots, whereas the difference $\text{LHS}-\text{RHS}$ has three.

Case I: $\beta_{0}\neq 0$. By the above statement we easily see that if $\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\notin\mathbf{vN}$, then $P_{\beta_{1},\beta_{0}}(\cdot,\xi,0)\notin\mathbf{svN}$ for $|\xi|$ large enough. We now show that

[TABLE]

So let us suppose in the rest of Case I that $(\beta_{1},\beta_{0})\notin W$ and $\beta_{0}\neq 0$.

Case I${}_{\text{a}}$. First we check the case when $\mathfrak{cc}\,\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)=0$. Then

[TABLE]

and, since now $2\beta_{1}-3\neq 0$, we can apply Theorem 2.4 to the above polynomial in $[\cdots]$: due to $[\cdots]^{\mathbf{r}}\equiv 0$ we have that $[\cdots]\in\mathbf{vN}$ if and only if $\zeta\mapsto[\cdots]^{\prime}=2\left(2\beta_{1}-3\right)\zeta-4\beta_{1}\in\mathbf{vN}$. But we directly see that this last linear polynomial $\notin\mathbf{vN}$, because $(\beta_{1},\beta_{0})\notin W$ and $\mathfrak{cc}\,\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)=0$ imply $\beta_{1}>3/4$.

Case I${}_{\text{b}}$. The conditions $\mathfrak{cc}\,\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\neq 0\neq\beta_{0}$ mean that we can apply Theorem 2.4 to $\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)$. It is easy to verify that $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}$ does not vanish identically, so $\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\in\mathbf{vN}$ if and only if

[TABLE]

We show in Cases I${}_{\text{b1}}$ and I${}_{\text{b2}}$ below that (27) never occurs. First we observe that the inequality constraint in (27) yields that $\mathfrak{lc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\neq 0$.

Case I${}_{\text{b1}}$. If $\mathfrak{cc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}=0$, then the polynomial $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\zeta)\right)^{\mathbf{r}}$ has exactly two roots: $\zeta_{1}=0$ and

[TABLE]

One directly checks that $\mathfrak{cc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}=0$ and $(\beta_{1},\beta_{0})\notin W$ imply $|\zeta_{2}|>1$.

Case I${}_{\text{b2}}$. If $\mathfrak{cc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\neq 0$, we apply Theorem 2.4 to get that the quadratic polynomial $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\in\mathbf{vN}$ if and only if either Case I${}_{\text{b2}\alpha}$ or I${}_{\text{b2}\beta}$ below occurs.

Case I${}_{\text{b2}\alpha}$: when $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{rr}}\equiv 0$ and $\left[\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\right]^{\prime}\in\mathbf{vN}$. In this case, however, the unique root of the polynomial $\left[\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\right]^{\prime}$,

[TABLE]

has absolute value $>1$.

Case I${}_{\text{b2}\beta}$: when $\left|\mathfrak{lc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\right|>\left|\mathfrak{cc}\,\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{r}}\right|$ and $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{rr}}$ $\in\mathbf{vN}$. But then the unique root of $\left(\text{LHS}_{\beta_{1},\beta_{0}}(\cdot)\right)^{\mathbf{rr}}$ is

[TABLE]

for which we again have $|\zeta_{1}|>1$, completing Case I.

Case II: $\beta_{0}=0$. Then

[TABLE]

and the leading coefficient of this cubic polynomial is 1. For each fixed $\beta_{1}\in\mathbb{R}$ we see that at least one of its coefficients is unbounded as $\xi\to-\infty$, so (by Vieta’s formulae) at least one of its roots $\zeta(\xi)$ is unbounded as $\xi\to-\infty$. Hence $(-\infty,0)\times\{0\}\subset{\mathcal{S}}_{\beta_{1},0}$ cannot hold. ∎

A.2 The proof of Theorem 5.3

Proof.

In the proof we suppose $m>0$ and, due to Lemma 5.1, that $(\beta_{1},\beta_{0})\in W$.

Step 1. Let us apply the same ideas as in Section A.1 but along the ray $\eta=-m\xi$. For $\xi<0$ we consider the roots of $P_{\beta_{1},\beta_{0}}(\cdot,\xi,-m\xi)$, and get that

[TABLE]

for some $|\xi|$ large enough, where the corresponding “modified left-hand side” is defined as

[TABLE]

and we have also taken into account that $\mathfrak{lc}\,\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)=\beta_{0}\neq 0$. (The corresponding “modified right-hand side” would be the same $\text{RHS}_{\xi}(\zeta)$ as in Section A.1.) Hence if the inclusion (22) holds with some $m>0$, then $\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\in\mathbf{vN}$.

Step 2. In this step we derive a necessary condition for $\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\in\mathbf{vN}$. First, one simply checks via Theorem 2.4 that $\mathfrak{cc}\,\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)=0$, $\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\in\mathbf{vN}$ and $m>0$ cannot be simultaneously true. So we can suppose

[TABLE]

We check that $\left(\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\right)^{\mathbf{r}}$ does not vanish identically, and that

[TABLE]

Then by Theorem 2.4 we have that

[TABLE]

Now we see that

[TABLE]

and $\left(\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\right)^{\mathbf{rr}}$ does not vanish identically. Thus Theorem 2.4 yields that

[TABLE]

if and only if

[TABLE]

and

[TABLE]

Clearly, $\deg\left(\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\right)^{\mathbf{rr}}\leq 1$, and we directly confirm that (28) implies that the degree is exactly 1. From this we obtain that (28) and (29) hold if and only if (28) and

[TABLE]

[TABLE]

hold. In particular, (28) guarantees that the denominators appearing in (30) are non-zero, hence from now on we can restrict the parameters $(\beta_{1},\beta_{0})\in W$ to the set $(\beta_{1},\beta_{0})\in W\setminus L$ with

[TABLE]

being the left edge of the wedge $W$; see Figure 10.

By defining

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

it is easily verified after some factorization and simplification that

[TABLE]

In particular, $\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\in\mathbf{vN}$ implies $Q_{\beta_{1},\beta_{0}}(m)\geq 0$.

Step 3. We see that $C_{4}<0$ and $C_{0}\geq 0$ for $(\beta_{1},\beta_{0})\in W\setminus L$, hence we can denote the largest real root of the polynomial $Q_{\beta_{1},\beta_{0}}(\cdot)$ by $m^{*}(\beta_{1},\beta_{0})\in[0,+\infty)$. Consequently, if $\text{MLHS}_{\beta_{1},\beta_{0},m}(\cdot)\in\mathbf{vN}$, then $m\leq m^{*}(\beta_{1},\beta_{0})$. We now conjecture (by using Mathematica’s Maximize command, for example) that

[TABLE]

and $m^{*}(\beta_{1},\beta_{0})=1/2$ occurs precisely for $(\beta_{1},\beta_{0})=({3}/{8},{3}/{4})$, see Figure 11. With this conjectured optimal $m^{*}$ value, we can prove (32) and the uniqueness property in an elementary way.

By introducing the shifted variable $M:=m-1/2$, we rewrite $Q_{\beta_{1},\beta_{0}}(m)$ as

[TABLE]

Then we check that

[TABLE]

Moreover, we have

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

On the one hand, these mean that (33) is negative for $M>0$ and $(\beta_{1},\beta_{0})\in W\setminus L$. On the other hand, for $M=0$ the polynomial (33) is zero if and only if $(\beta_{1},\beta_{0})=(3/8,3/4)$.

Therefore we have proved that if (22) holds with some $m>0$, then $m\leq 1/2$; and if $m=1/2$ is possible at all, then $(\beta_{1},\beta_{0})=({3}/{8},{3}/{4})$.

Step 4. In this final step we show that $m=1/2$ in (22) can be achieved, by showing that ${\mathcal{A}}_{1/2}\subset{\mathcal{S}}_{{3}/{8},{3}/{4}}$, that is,

[TABLE]

Let us fix such a pair $(\xi,\eta)$. One sees that

[TABLE]

and in the $\mathfrak{cc}\,P_{3/8,3/4}(\cdot,\xi,\eta)=0$ case (34) is easily verified to hold. Otherwise, if $\mathfrak{cc}\,\neq 0$, we check that $\left(P_{3/8,3/4}(\cdot,\xi,\eta)\right)^{\mathbf{r}}$ does not vanish identically, so by Theorem 2.5 we have that $P_{3/8,3/4}(\cdot,\xi,\eta)\in\mathbf{svN}$ if and only if

[TABLE]

We have that $\mathfrak{lc}\,\left(P_{3/8,3/4}(\cdot,\xi,\eta)\right)^{\mathbf{r}}\neq 0$. Moreover, $\mathfrak{cc}\,\left(P_{3/8,3/4}(\cdot,\xi,\eta)\right)^{\mathbf{r}}=0$ for $\xi=-2$ or $\xi=-2/3$, in which cases (35) holds. So we can suppose from now on that $\mathfrak{cc}\,\left(P_{3/8,3/4}(\cdot,\xi,\eta)\right)^{\mathbf{r}}\neq 0$. Then one proves that

[TABLE]

[TABLE]

so $\deg\left(P_{3/8,3/4}(\cdot,\xi,\eta)\right)^{\mathbf{rr}}=1$. Hence, by using Theorem 2.5 again, we get that (35) holds if and only if

[TABLE]

hold. Finally, we check that these last two conditions are satisfied for any $(\xi,\eta)\in{\mathcal{A}}_{1/2}$ pair not excluded earlier during the case separations. ∎

Remark A.1**.**

By defining

[TABLE]

[TABLE]

and applying Theorem 2.5, it is straightforward to show (cf. Step 4 in the above proof) that

[TABLE]

see Figure 12 and cf. Remark B.1.

Appendix B Appendix

B.1 Locating the candidate optimum for Theorem 6.1

For a given $(\beta_{1},\beta_{0})\in W$ pair, we can represent the RLC of the corresponding multistep method of the family (18) as an implicit curve of the form

[TABLE]

by using the transformations in Section 3.1 as follows. First we perform the substitution $\zeta\mapsto\frac{i-t}{i+t}$ in the polynomial (20), then eliminate $t\in\mathbb{R}$ by taking the resultant of the real and imaginary parts of $P_{\beta_{1},\beta_{0}}\left(\frac{i-t}{i+t},\xi,\eta\right)$. The resulting polynomial can be factored to get $2^{34}\cdot 9\cdot\left(1-\beta_{0}\xi\right)^{6}\cdot F_{\beta_{1},\beta_{0}}(\xi,\eta)$; the normalization with $F_{\beta_{1},\beta_{0}}(0,1)=9$ has been used to make this polynomial $F_{\beta_{1},\beta_{0}}$ unique. The term $\left(1-\beta_{0}\xi\right)^{6}$ (cf. the leading coefficient of $P_{\beta_{1},\beta_{0}}(\cdot,\xi,\eta)$) does not vanish now due to $\xi\leq 0$ and $(\beta_{1},\beta_{0})\in W$, hence (36) is obtained. We are not going to display the polynomial $F_{\beta_{1},\beta_{0}}(\xi,\eta)$ explicitly: it contains 82 terms in its expanded form and its degree in the variables/parameters $(\xi,\eta,\beta_{1},\beta_{0})$ is $(6,4,4,4)$.

Now supposing that the RLC (36) describes the boundary of the stability region of the multistep method determined by the given pair $(\beta_{1},\beta_{0})$, it is reasonable to expect that, say, the upper branch of the largest parabola inscribed in ${\mathcal{S}}_{\beta_{1},\beta_{0}}$, $\{(\xi,\eta)\in\mathbb{R}^{2}:\xi<0,\,\eta>0,\,\eta^{2}=-m\xi\}$, touches the RLC (36) at some finite point. In this case, the polynomial

[TABLE]

has a multiple root there—it is indeed a polynomial, because in our situation $F_{\beta_{1},\beta_{0}}(\xi,\eta)$ contains only even powers of $\eta$ (namely, $\eta^{2}$ and $\eta^{4}$). Moreover, we now have

[TABLE]

where $\widetilde{Q}_{\beta_{1},\beta_{0},m}(\cdot)$ is a quartic polynomial. The existence of a multiple root of $\widetilde{Q}_{\beta_{1},\beta_{0},m}(\cdot)$ implies that the discriminant of this polynomial (with respect to $\xi$), denoted by $\widetilde{\Delta}_{\beta_{1},\beta_{0}}(m)$, vanishes. Mathematica yields that

[TABLE]

[TABLE]

where the “$\ldots$” symbols contain 57 and 228 terms, respectively. We see that the factor $9\beta_{0}+4\beta_{1}-6$ above is always positive in $W$. In this way we can determine the parameter $m$ of the largest parabola within the region bounded by the RLC for any fixed $(\beta_{1},\beta_{0})\in W$.

Remark B.1**.**

By setting $(\beta_{1},\beta_{0})=(3/8,3/4)$ for example (corresponding to the “sector-optimal” method in Section A.2), we have that the RLC in (36) is identical to $3/16\cdot F_{\text{opt}}(\xi,\eta)$ in Remark A.1, implying that the RLC in the left half-plane $\xi\leq 0$ represents the boundary of the stability region ${\mathcal{S}}_{3/8,3/4}$. Now $\widetilde{\Delta}_{3/8,3/4}(m)$ can be written as

[TABLE]

from which we can prove that the largest parabola ${\mathcal{P}_{m}}$ contained in ${\mathcal{S}}_{3/8,3/4}$ has $m\approx 1.11226$ (being the unique positive root of the polynomial $\{36,1362,343,-2116\}$).

By studying the positive roots of $\widetilde{\Delta}_{\beta_{1},\beta_{0}}(\cdot)$ as $(\beta_{1},\beta_{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 $\beta_{1}=1/5$ and $\beta_{0}=37/40$, and in this case the RLC and the upper parabola branch touch each other at $(\xi,\eta)=(-10/7,2\sqrt{3/7})$. Since the polynomial $\widetilde{\Delta}_{\beta_{1},\beta_{0}}(m)$ is much more complicated than the corresponding polynomial $Q_{\beta_{1},\beta_{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.

Step 1. By assuming $(\beta_{1},\beta_{0})\in W$ throughout the step, we show that the point $(\xi_{0},\eta_{0}):=(-10/7,2\sqrt{3/7})$ belongs to precisely one stability region in the family, by verifying that

[TABLE]

To see this, first we check that $\mathfrak{lc}\,P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)\neq 0$. Moreover, it is easily seen that $\mathfrak{cc}\,P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)$ vanishes exactly for $\beta_{1}\leq 23/40$ and $\beta_{0}=\left(67-40\beta_{1}\right)/80$, and in this case the polynomial $P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)/\xi$ has $\deg=2$ but $\notin\mathbf{svN}$, as a recursive application of Theorem 2.5 shows. Then we can also prove that $\left(P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)\right)^{\mathbf{r}}$ does not vanish identically, and that

[TABLE]

Thus, according to Theorem 2.5,

[TABLE]

Now we repeat the above process with $\left(P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)\right)^{\mathbf{r}}$. We prove that

[TABLE]

and that $\left(P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)\right)^{\mathbf{rr}}$ does not vanish identically, so by Theorem 2.5 we have that

[TABLE]

But $\left(P_{\beta_{1},\beta_{0}}\left(\cdot,\xi_{0},\eta_{0}\right)\right)^{\mathbf{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 $\leq 1$ if and only if

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The product of the first two factors is strictly negative in $W$, and a standard constrained optimization computation shows that the third factor is $\geq 0$ in $W$ if and only if $(\beta_{1},\beta_{0})=\left(1/5,37/40\right)$, completing Step 1.

Step 2. Since ${\mathcal{P}_{m_{1}}}\subseteq{\mathcal{P}_{m_{2}}}$ is equivalent to $0<m_{1}\leq m_{2}$ (see (24)), and now $\left|\eta_{0}^{2}/\xi_{0}\right|=6/5$, the uniqueness property in the previous step implies that $m\geq 6/5$ in (25) can hold only for $(\beta_{1},\beta_{0})=\left(1/5,37/40\right)$. In this step we verify that (25) indeed holds with $m=6/5$ and $(\beta_{1},\beta_{0})=\left(1/5,37/40\right)$, that is, we show that $P_{1/5,37/40}\left(\cdot,\xi,\eta\right)\in\mathbf{svN}$ for any $(\xi,\eta)\in{\mathcal{P}_{6/5}}$.

Let us pick and fix an arbitrary point $(\xi,\eta)\in{\mathcal{P}_{6/5}}$. Then we easily see that

[TABLE]

and this $\mathfrak{cc}\,=0$ if and only if $\xi=-5/11$; in this case Theorem 2.5 tells us that $\zeta\mapsto P_{1/5,37/40}\left(\zeta,-5/11,\eta\right)=\zeta(125\zeta^{2}-132i\zeta\eta-58\zeta-7)/88\in\mathbf{svN}$. So for $\xi\neq-5/11$, again by Theorem 2.5 we get that

[TABLE]

provided that $\left(P_{1/5,37/40}\left(\cdot,\xi,\eta\right)\right)^{\mathbf{r}}$ does not vanish identically. But this non-vanishing condition is true because

[TABLE]

Moreover, since $\left|\mathfrak{lc}\,\left(P_{1/5,37/40}\left(\cdot,\xi,\eta\right)\right)^{\mathbf{rr}}\right|$ is also positive, the above with Theorem 2.5 imply that

[TABLE]

The positivity of $\left|\mathfrak{lc}\,\left(\ldots\right)^{\mathbf{rr}}\right|$ yields that $\left(P_{1/5,37/40}\left(\cdot,\xi,\eta\right)\right)^{\mathbf{rr}}$, a $\deg=1$ polynomial, has a unique root. The absolute value of this (real or complex) root is $\leq 1$ if and only if $(3\xi-10)(59\xi-30)\cdot\widetilde{F}_{\text{opt}}(\xi,\eta)\leq 0$, where

[TABLE]

[TABLE]

Now $(3\xi-10)(59\xi-30)>0$, and one checks that $\widetilde{F}_{\text{opt}}(\xi,\eta)\leq 0$ for $(\xi,\eta)\in{\mathcal{P}_{6/5}}$, completing Step 2.

Step 3. To complete the optimality proof, we finally show that

[TABLE]

that is, we cannot have $m>6/5$ in (25). We repeat the same two-step reduction process as above and get that $\varepsilon\in(0,1)$ guarantees that

[TABLE]

But this last $\left(\ldots\right)^{\mathbf{rr}}\in\mathbf{svN}$ condition is equivalent to

[TABLE]

so it cannot hold for any $\varepsilon\in(0,1)$.

Remark B.2**.**

In addition to the inequality $\widetilde{F}_{\text{opt}}(\xi,\eta)\leq 0$ in Step 2, we have that $\widetilde{F}_{\text{opt}}(\xi,\eta)=0$ for $(\xi,\eta)\in{\mathcal{P}_{6/5}}$ if and only if $(\xi,\eta)=(\xi_{0},\pm\eta_{0})$. Moreover, $\widetilde{F}_{\text{opt}}(\xi,\eta)\equiv 2000\cdot F_{1/5,37/40}(\xi,\eta)$ (see (36)). On the other hand, by using the reduction process one can actually prove that

[TABLE]

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 (and cf. Remarks 5.4 and B.1).