{'fraction_non_alphanumeric': 0.05227210323997803, 'fraction_numerical': 0.016474464579901153, 'mean_word_length': 3.9217905405405404, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On computing the instability index of a non-selfadjoint differential operator associated with coating and rimming flows

We study the problem of finding the instability index of certain non-selfadjoint fourth order differential operators that appear as linearizations of coating and rimming flows, where a thin layer of fluid coats a horizontal rotating cylinder. The main result reduces the computation of the instability index to a finite-dimensional space of trigonometric polynomials. The proof uses Lyapunov's method to associate the differential operator with a quadratic form, whose maximal positive subspace has dimension equal to the instability index. The quadratic form is given by a solution of Lyapunov's equation, which here takes the form of a fourth order linear PDE in two variables. Elliptic estimates for the solution of this PDE play a key role. We include some numerical examples.


Introduction
The stability of steady states is a basic question about the dynamics of any partial differential equation that models the evolution of a physical system. Frequently, the first step is to linearize the system about a given equilibrium. Linearized stability is determined by the spectrum of the resulting differential operator A. If A has discrete spectrum, an important quantity is the instability index, κ(A), which counts the number of eigenvalues in the right half plane (with multiplicity).
In order to numerically evaluate the instability index of a given differential operator, its computation should be reduced to a problem of linear algebra. Particularly for problems with periodic boundary conditions, it seems natural to restrict A to a finite-dimensional space of trigonometric polynomials. Under what conditions can κ(A) be computed from the resulting finite matrix? One difficulty is that the entries of the infinite matrix corresponding to the differential operator A grow with the row and column index, so that any truncation is not a small perturbation.
If A is a self-adjoint semi-bounded differential operator of even order, then the computation of its instability index is well-understood through the classical work of Morse [18] who solved this problem completely in the space of vector functions in one independent variable. The instability index of A agrees with the dimension of the positive cone of the corresponding quadratic form. It is invariant under congruence transformations that replace A with T * AT . The instability index can be estimated by variational methods, or computed directly from the zeroes of the corresponding Evans function.
Understanding the spectrum of a non-selfadjoint operator is a much harder problem. It is not at all obvious how to restrict the computation of its instability index to a finite-dimensional subspace, or how to even estimate its dimension. Furthermore, the numerical calculation of eigenvalues can be extremely ill-conditioned even in finite dimensions. One impressive example is the matrix Note that B differs from A only in the signs of two off-diagonal entries. The chance of encountering a matrix with moderately-sized entries and a badly conditioned eigenvalue problem increases rapidly with the dimension of the matrix (see [14,24]). Such examples demonstrate that the stability problem for a non-selfadjoint operator cannot be easily solved by direct computations of the spectrum.
In this paper, we examine the computation of the instability index for differential operators of the form acting on 2π-periodic functions. Such operators appear as linearizations of models for thin liquid films moving on the surface of a horizontal rotating cylinder. The resulting flows are called coating, if the fluid is on the outside of the cylinder, and rimming, if the fluid is on the inside of a hollow cylinder. They appear in many applications, including coating of fluorescent light bulbs when a coating solvent is placed inside a spinning glass tube, different type of moulding processes and paper productions. One would expect the flow to become unstable, if the fluid film is thick enough so that drops of fluid can form on the bottom of the cylinder (in case of a coating flow) or on its ceiling (in case of a rimming flow). In both cases, surface tension and higher rotation speeds should help to stabilize the fluid, but may also allow for more complicated steady states.
The operators in Eq.(1.1) appear as linearizations of the flows about steady states, when the dependence on the longitudinal variable in the cylinder is neglected. Benilov, O'Brien and Sazonov [5] studied the convection-diffusion equation with periodic boundary conditions on [0, 2π], which corresponds to a singular limit of a rimming flow where surface tension is neglected. This operator has remarkable properties: For |ε| < 2, all its eigenmodes are neutrally stable, but the Cauchy problem is ill-posed in any Sobolev and Hölder space of 2π-periodic functions. The underlying cause is the sign change of the diffusion coefficient as x → x + π. This phenomenon of explosive instability of a system with purely imaginary spectrum was studied analytically by Chugunova, Karabash and Pyatkov [9], who explained it in terms of the absence of the Riesz basis property of the set of eigenfunctions. The spectral and asymptotic properties of A are of interest in operator theory and were analyzed in [12,26,8,10]. One should expect the explosive instability to disappear in complete models that includes the smoothing effect of surface tension. Such models have been proposed, for example, by [20,21]. In [6,7], the authors linearized this model about some approximation of a positive steady state solution to obtain with periodic boundary conditions. Here, the parameter α 1 is related to the gravitational drainage, α 2 is related to the hydrostatic pressure (in lubrication approximation model this coefficient is very small), and the parameter α 3 describes surface tension effect. They showed numerically that a sufficiently strong surface tension can stabilize the film provided that the other coefficients are not too small. For smaller values of α 1 and α 2 , capillary effects destabilize the film. The number of unstable eigenvalues of A grows if α 3 is decreased.
We will consider operators given by Eq. (1.1) acting on L 2 [0, 2π] with periodic boundary conditions. We assume that the coefficients a(x), b(x) and c(x) are bounded smooth periodic functions. We will show that the instability index of A is determined by its projection to a sufficiently large finite-dimensional subspace of L 2 . The dimension of the space depends on a suitable norm of the distributional solution U of the partial differential equation with periodic boundary conditions on [0, 2π] × [0, 2π]. Here, the differential operator A is defined by applying the single-variable differential operator A to the functions F (·, y) and F (x, ·) and adding the results; symbolically We note that Eq. (1.3) has a unique solution if the spectra of A and −A * are disjoint [2]. Let U 0 (x, y) be the solution of Eq. (1.3) with a(x) = b(x) = 0 and c(x) = 1. We will see below that U 0 is piecewise smooth, with a jump in the third derivative across the line x = y, and that U(x, y) − U 0 (x, y) ∈ H 4 . To describe our results, denote by P N the standard projection onto the space of trigonometric polynomials of order N, 5) provided that The constant is given by The significance of Eq. (1.5) is that it allows to compute the instability index of A from the finite matrix that describes the restriction of U to the finite-dimensional subspace The weakness of this result is that both the condition on N and the computation of the subspace involve the unknown function U, which is defined as the solution of a partial differential equation.
The existence of such a solution, and its norm, depend sensitively on the spectrum of A, which is exactly the unknown quantity we are concerned with. It is tempting to consider instead the matrix obtained by truncating the Fourier representation of A at a suitable high order N. Our main result, Proposition 7.4, guarantees that ( 1.8) provided that Note that only the norm of the unknown function enters into the condition on N, and that the identity in Eq. (1.8) does not involve U at all. The selection of N and the problem of estimating this norm will be discussed at the end. Let us add a few words about the proofs. Our analysis relies on the indefinite quadratic form defined by the self-adjoint operator U. Classical results, which will be discussed in the next section, state that κ(A) = κ(U) , and that the positive and negative cones of U contain the invariant subspaces associated with the spectrum of A in the right and left half planes, respectively. The key to Eq. (1.5) is that the quadratic form is negative on high Fourier modes, because the fourth order term in A dominates the lower order derivatives. As part of the argument, we derive an addition formula for the instability index of a self-adjoint operator in terms of its restriction to suitable subspace. The proof of Eq. (1.8) combines Eq. (1.5) with estimates for the off-diagonal terms in the Fourier representation for U.
One of the possible extensions of our results could be an application of a similar method to obtain the estimations on the size of the finite dimensional truncation in the case of a more general forth order differential operator with the third order derivative term which is absent in 1.2.

Lyapunov's equation
The partial differential equation (1. The problem of obtaining information about the sign of eigenvalues of A in situations where both V and U may be indefinite and have non-trivial kernels remains an area of active research.

3) is an instance of Lyapunov's equation
Lyapunov's equation has many applications in stability theory and optimal control. In typical applications, κ(A) = 0, so that the system is asymptotically stable, and U is used to estimate the rate of convergence. Eq. (2.1) is a special case of Sylvester's equation which has been studied extensively in Linear Algebra, Operator Theory, and Numerical Analysis. It is known to be uniquely solvable, if and only if the matrices A and B have no eigenvalues in common. In particular, Eq. (2.1) has a unique solution if the spectra of A and −A * are disjoint. Since V is self-adjoint, a unique solution U is automatically self-adjoint as well. These results were extended to bounded operators on infinite-dimensional Hilbert spaces by Daleckii and Krein [11] and to unbounded operators by Belonosov [2,3].
Before stating Belonosov's result, we recall that a closed densely defined operator A on a Banach space is sectorial, if the spectrum of A is contained in an open sector with vertex at λ 0 ∈ R and opening angle θ < π/2, and the resolvent R λ (A) = (A − λI) −1 is uniformly bounded for λ outside S. Sectorial operators are precisely the generators of analytic semigroups. The sector S is invariant under similarity transformations, and does not change if the norm on the space is replaced by a equivalent norm.

Then for any bounded operator V on H, the Lyapunov equation (2.1) has a unique solution U in the class of bounded operators on H. Then U is invertible in the general sense, i.e. its inverse is densely defined but can be unbounded operator
Belonosov actually proved more general existence and uniqueness results for the Sylvester's equation in Banach spaces.
To explain the geometric meaning of Lyapunov's equation, we introduce on H the indefinite inner product If U has trivial nullspace and κ(U) < ∞, then H equipped with [·, ·] is called a Pontryagin space, and will be denoted by Π. The concepts of orthogonality and adjointness are defined on Π in the natural way with respect to the indefinite inner product Lyapunov's equation guarantees that the value of the quadratic form Denote by M + (A) the invariant subspace associated with the part of the spectrum of A located in the right half plane. If φ ∈ M + (A), then , which corresponds to the spectrum of A in the left half plane, is a negative subspace of Π. Since Eq. (2.1) excludes purely imaginary eigenvalues, these subspaces are maximal, and consequently κ(A) = κ(U). One can also interpret Lyapunov's equation as a dissipativity condition on A with respect to the Pontryagin space Π. In general, a densely defined linear operator i.e., −A is dissipative. In this framework, the analogue of Belonosov's theorem was proven by Azizov [1] (but note that Azizov formulates the result in terms of Im rather than Re ): (Azizov). Let A be an operator on Π such that −A is maximally dissipative. Then there exist a maximal nonnegative subspace Π + and a maximal nonpositive subspace Moreover, we can choose Π + and Π − to be invariant subspaces for A, and The second part of Azizov's theorem implies that κ(A) = κ(U) provided that V in Eq. (2.1) is positive definite. This agrees with the conclusion of Theorem 2.2, but note the difference in the hypotheses: Belonosov's assumption that A is sectorial provides resolvent estimates that allow to represent U as a contour integral (thereby proving existence), and the analytic semigroup e tA appears in the proof that κ(A) = κ(U), as sketched above. In contrast, Azizov's theorem does not require A to be sectorial, but starts instead from a given solution to Eq. (2.1). In the special case where κ(U) = 0, Theorem 2.3 reduces to a theorem of Phillips that characterizes maximal dissipative operators as generators of strongly continuous contraction semigroups. In particular, the spectrum of A lies in the closed left half plane (see [25], Corollary 1 in Section IX.4).
In the case where A is a sectorial differential operator of even order on an interval [a, b] Belonosov proved that the solution of Lyapunov's equation with V = I is given by a self-adjoint bounded operator [4]. His results are formulated for "split" boundary conditions that do not couple the values at the two ends of the interval. Belonosov's results were extended to second-order sectorial differential operators with non-split boundary conditions by Tersenov [23]. The operators we consider here are of fourth order with periodic boundary conditions.
It is an interesting open question how to take advantage of the freedom to choose an arbitrary positive definite self-adjoint bounded operator V for the right hand side of Eq. (2.1). For instance, if A is a sectorial non-selfadjoint differential operator, can V be chosen in such a way that the solution U is the inverse of a differential operator?

Spaces and norms
We start with some estimates for the differential operator in Eq. (1.1). We will work in L 2 = L 2 [0, 2π] , and will use periodic boundary conditions throughout. The inner product and norm are denoted by For the Fourier coefficients we use the conventionŝ In the hope of minimizing confusion, we will denote functions on [0, 2π] by lowercase letters such as (f , φ, . . . ), and functions on the square [0, 2π] × [0, 2π] by uppercase letters (F , Φ, ....). Abusing notation, we will identify a function F (x, y) ∈ L 2 with the corresponding integral operator F on L 2 . By Schwarz' inequality, and consequently Operators on functions of two variables will be denoted by calligraphic letters (F , G, . . . ). Given a single-variable operator F , we denote by F x or F y the operators that acts on the xor y-variable of a function Φ(x, y) while keeping the other one fixed. On the Sobolev spaces H s = H s [0, 2π] with periodic boundary conditions, we use the norms The corresponding Sobolev spaces of doubly periodic functions on [0, 2π] × [0, 2π] will be denoted by H s , and their norms are defined by Note that for s = 0, this agrees with the definition of the L 2 -norm as the square integral. The choice of the Fourier multipliers (1 + p 4 ) s/2 and (2 + p 4 + q 4 ) s/2 in place of the standard (1 + p 2 ) s and (1 + p 2 + q 2 ) s allows for an easier comparison between functions of one and two variables. Finally, we denote by D the unique positive definite self-adjoint operator on L 2 such that This is a first-order pseudodifferential operator that provides an isometry from H s+1 onto H s for every value of s. The domain of the operator A in Eq. (1.2) consists of periodic functions in H 4 [0, 2π], and its adjoint is given by Lemma 3.1. For any a ∈ H 1 and every φ ∈ L 2 , we have In particular, where M is the constant from Eq. (1.6).
For the second claim, we use that for φ ∈ H 2

Lemma 3.2. A is sectorial.
Proof. It suffices to show that the Hausdorff set f, Af | f ∈ Dom(A), ||f || = 1 is contained in a closed sector with some vertex λ 0 and opening angle θ < π, and that A − λ 0 I is invertible (see p. 280 of [17]). Choose We estimate, for f ∈ H 2 This shows that the spectrum of A lies in the half plane Re z < λ 0 − 1 2 . Similarly, For ||f || = 1 it follows that which yields the claim.
The lemma implies that the Cauchy problem for A has a unique solution for every initial value h 0 ∈ L 2 . This solution is analytic in t for t > 0, and for any fixed t > 0, the function h(t, ·) ∈ Dom(A). If the coefficients of A are analytic, then h is analytic in both variables for t > 0. An application of the Lax-Milgram theorem similar to Lemma 4.1 below shows that (λ 0 − A) −1 maps L 2 into H 2 . It follows that the resolvent is a compact operator of the Hilbert-Schmidt type, and that the spectrum of A is discrete.

The integral kernel U (x, y)
Let A be the differential operator from Eq. (1.1). Theorem 2.2 implies that Lyapunov's equation has a unique solution U, provided that the spectra of A and A * are disjoint. Our goal is to show that U admits an integral representation and to derive bounds on U(x, y). Equation (2.1) requires that (A * U + UA * )φ, ψ = φ, ψ for all smooth periodic test functions functions φ, ψ. This means that U(x, y) is a distributional solution of the partial differential equation (1.3).
Let us solve Eq. (1.3) in the special case A 0 = −D 4 , given by By our choice of norms, −A 0 defines an isometry from H 4 onto L 2 . Since A 0 has constant coefficients, the unique solution can be written as U 0 (x, y) = u 0 (x − y), where in other words, 2U 0 (x, y) is the Green's function of A 0 on [0, 2π] with periodic boundary conditions. One can compute u 0 (x) explicitly as a linear combination where the coefficients are adjusted so that u 0 is periodic and twice differentiable, and its third derivative jumps by −1/2 at x = 0. From this representation, it is clear that U 0 is smooth away from the line x = y, and that U 0 ∈ C 2,1 ⊂ H 3 . Alternately, we easily obtain from the Fourier representation of A 0 thatû 0 (p) = − 1 4π(1+p 4 ) , and (4.1) In particular, U 0 (x, y) ∈ H s for all s < 7 2 , and ||U 0 || H 3 ≤ 1. It remains to analyze the difference K(x, y) := U(x, y) − U 0 (x, y) .
By definition, K solves the partial differential equation The second order differential operator A * − A 0 maps maps H 2 into L 2 , see Lemma 3.1. A weak solution of this equation is provided by the next lemma.

Lemma 4.1 (Construction of K). The resolvent of A * is compact and maps
Proof. Let λ 0 be the vertex of the sector computed in Lemma 3.2, and assume that F (x, y) ∈ L 2 . We verify that the equation satisfies the assumptions of the Lax-Milgram theorem, as stated in [Evans, PDE, p. 297] [13]. Define a bilinear form on on smooth doubly periodic functions Φ, Ψ Then B is extended continuously to H 2 by On the other hand, it follows from Eq. (3.3) that defines a continuous linear form on H 2 . The Lax-Milgram theorem asserts that there exists a unique function K(x, y) ∈ H 2 such that for all Ψ ∈ H 2 . By the resolvent identity, the equation has a unique weak solution in H 2 for every value of λ that is not an eigenvalue of A * and every F (x, y) ∈ L 2 .
Proof. If K(x, y) solves Eq. (4.2), then and we conclude that The proof is completed with Lemma 3.1.

Estimates for the operator U
In this section, we derive bounds for U = U 0 + K as an operator on L 2 . Since K(x, y) ∈ H 4 , while U 0 (x, y) ∈ H s only for s < 7/2, the Fourier coefficients of K(x, y) decay more quickly than the Fourier coefficients of U 0 (x, y). This in turn implies that the restriction of U to high Fourier modes is dominated by U 0 . In this section, we provide the relevant estimates. As a consequence of the regularity result in Lemma 4.2 we see that U defines a bounded linear operator from L 2 to H 4 , with We have used that D 4 x U 0 = 1 2 δ and applied Eq. (3.1) to D 4 K(x, y). One attractive property of the H 4 -norm is that it depends only on the magnitude of the Fourier coefficients, not on the phases. In contrast, the operator norm can change drastically if we replaceF (p, q) by |F (p, q)|. This dependence on cancelations can cause difficulties in estimates: Multiplying the Fourier coefficients of F with factors α(p, q) ∈ [0, 1] will not necessarily decrease the operator norm. On the other hand, the H 4 -norm provides only a rather loose bound on the norm of the corresponding integral operator. For instance, the kernles U 0 (x, y) (and consequently U(x, y) does not lie in H 4 , even though ||U 0 || L 2 →H 4 = 1 2 . We find it useful to introduce another norm on integral kernels that lies between the H 4 -norm (as a function of two variables), and the operator norm (as a linear transformation from L 2 to H 4 ). By construction, this norm depends only on the modulus of the Fourier coefficients.

Then
|||F ||| ≤ ||F (x, y)|| H 4 , Proof. From the Fourier representation, we see that On the other hand, and similarly We note that if F has positive Fourier coefficients, then |||F ||| agrees with the operator norm of A 0 F as a linear transformation from L 2 into itself. In particular, |||U 0 ||| = 1.
Proof. Using Eq. (4.2) together with the definition of the norm, we obtain The first summand is bounded below by 1 because A 0 U 0 = I, and the second summand is bounded by MN −2 |||U||| according to Lemma 5.2. We conclude that (1 − MN −2 )|||U||| ≥ 1 for each N, and the claim follows upon taking N → ∞.

Addition rule for the instability index
We return to the Pontryagin space Π introduced in Section 2, with the indefinite inner product given by Eq. (2.2). Let Π 1 be a finite-dimensional subspace of Π, and let be its U-orthogonal complement. By construction, dim Π 1 = codim Π 2 . The natural question is can we compute κ(U) from the restrictions κ(U| Π 1 ) and κ(U| Π 2 )? The difficulty is that Π need not be a direct sum of Π 1 and Π 2 , because the two subspaces may intersect non-trivially in a subspace where the quadratic form vanishes.
Two finite-dimensional neutral subspaces X and Y of H are Π-skewly linked, if dim X = dim Y and the inner product [., .] does not degenerate on the direct sum X+Y . In particular, no vector of X different from 0 is orthogonal to the skewly linked subspace Y , and vice versa. Theorem 6.1 (Theorem 3.4 [16]). Let Π 1 be an arbitrary subspace of H, let Π 2 be its U-orthogonal complement, and let X = Π 1 ∩ Π 2 be their intersection. There exists a neutral subspace Y ⊂ Π that is skewly linked to X and provides a U-orthogonal decomposition The theorem was originally formulated for the case of regular Pontryagin spaces, where the quadratic form U is a bounded operator with bounded inverse. Under the assumption that Π 1 is finite-dimensional, the result easily extends to the situation where the inverse of U is unbounded but densely defined. Although the above decomposition is not unique in general, it yields the following addition formula for instability indices: Proposition 6.2. Let Π 1 be a finite-dimensional subspace of Π, and let Π 2 be its U-orthogonal complement, Then its instability index is given by Proof. Theorem 6.1 provides subspaces Π ′ 1 and Π ′ 2 such that κ(U) = κ(U| Π ′ 1 ) + κ(U| Π ′ 2 ) + κ(U| X+Y ) ) . By construction, we have κ(U| Π ′ 1 ) = κ(U| Π 1 ) and κ(U| Π ′ 2 ) = κ(U| Π 2 ). It remains to compute κ(U| X+Y ) ).
Since X and Y are skewly linked and finite-dimensional, there exists for each basis φ 1 , φ 2 , ...φ m of X a basis ψ 1 , ψ 2 , ...ψ m of Y such that [φ i , ψ j ] = δ ij (i, j = 1, ..., m). By expanding an arbitrary element h ∈ X+Y as the indefinite inner product can be expressed as This is an explicit representation of the indefinite inner product in terms of positive and negative squares, which shows that κ(U| X+Y ) = dim(X).

Restriction to finite dimensions
We first prove the claim in Eq. (1.5).
For our final result, we want to replace Π 1 by the range of the projection P N from Eq. (1.4). The next two lemmas concern the restriction of U to the range of P N . Lemma 7.2 (Lyapunov equation for P N AP N ). Let A be given by Eq. (1.1). Assume that the spectrum of A and −A * are disjoint, and let U(x, y) be the kernel of the unique solution of Lyapunov's equation that we constructed in Section 4. If 4) where c N = 1 − M 2 N −4 |||U|||. In particular, Proof. For φ ∈ L 2 , we write φ 1 = P N φ, φ 2 = (I − P N )φ and decompose From Eq. (2.1), we see that U 11 solves Lyapunov's equation We claim that the right hand side is positive definite on the range of P N .
To prove this claim, first observe that we can replace A by A − A 0 and U by K = U − U 0 in the definition of V , because A 0 and U 0 are diagonal in the Fourier representation. We estimate by Lemma 5.2. It follows that V ≥ c N P N > 0 as quadratic forms on the range of P N . Since P N UP N is a finite matrix, the conclusion of the lemma follows with Taussky's theorem.

Lemma 7.3.
Under the assumptions of the previous lemma, P N UP N is invertible on the range of P N , and Proof. Let us write A N = D −2 P N AP N D 2 and U N = D 2 P N UP N D 2 . By Eq. (7.4), we have as quadratic forms on the range of P N . Here, c N = 1 − M 2 N −4 |||U|||, as in Lemma 7.2. We apply this inequality to an eigenfunction φ 0 of U N where λ 0 is the corresponding eigenvalue. Writing φ 0 = D −2 ψ 0 , and using once more Lemma 5.2, we conclude that We are finally ready for our main result.
where M is given by Eq. (1.6), then Proof. Since U solves Lyapunov's equation, Theorem 2.2 implies that κ(A) = κ(U), and we already know from Lemma 7.2 that κ(P N AP N ) = κ(P N UP N ). We want to apply Proposition 6.2 in the case where Π 1 is the range of P N . Since |||U||| ≥ 1 by Lemma 5.3, our assumption implies that δ N = MN −2 < 1. On we compute for the indefinite quadratic form By Eq. (7.3) of Proposition 7.1, the last term is negative on the nullspace of P N , and satisfies the bound as quadratic forms. To estimate the other summand, Lemma 5.2 yields and analogously The middle factor is controlled with Lemma 7.3 by We arrive at as quadratic forms. Since δ 2 N |||U||| 2 (1 + M) < 1 − δ N |||U||| 2 by assumption, U is negative definite on Π 2 . It follows from Proposition 6.2 that κ(U) = κ(P N UP N ), completing the proof.

Numerical examples
Before we look at examples, a few words about how to verify the hypothesis on N in Eq. (7.1) or Eq. (7.5). The conditions involve the solution of the partial differential equation 1.3. A useful consequence of Lemmas 5.2 and 5.3 is that for δ N := MN 2 < 1, This follows by using the triangle inequality |||U 0 + K||| ≤ 1 + |||P N KP N ||| + |||K − P N KP N ||| and solving for P N KP N in Lemma 5.2. We propose two ways to estimate the size of |||P N UP N |||.
• Solve the partial differential equation (1.3) by a Galerkin approximation, and use this solution to compute, approximately, the value of |||U|||. If Eq. (7.1) is satisfied for some value of N much below the dimension of the Galerkin approximation, we can apply Proposition 7.1, and restrict U to the subspace in Eq. (1.7). A basis for this subspace can be computed by using the Gram-Schmidt algorithm, with the inner product replaced by the indefinite inner product associated with U. If even Eq. (7.5) can be satisfied, then we can just restrict U to the range of P N .
• Start with a value of N such that δ N = MN −2 < 1. Write the matrix P N AP N in the Fourier representation, find its eigenvalues, and bring it into triangular form. Solve Lyapunov's equation P N A * P NŨN +Ũ N P N AP N = I for U N . In the Fourier representation, U N is a finite matrix. Compute λ max , the largest eigenvalue of the matrix (2 + p 4 + q 4 )|U N (p, q) − 1| |p|,|q|<N .
Then (1 + λ max )/(1 − δ N ) is our best estimate for |||U|||. If the condition in Eq. 7.5 holds with the current value of N, we are satisfied and accept the value of κ(A N ) as the instability index for A. Else, we increase N accordingly, and repeat the above steps. Proposition 7.4 reduces the computation of the stability index of A to a finite-dimensional linear algebra problem. This is illustrated in Fig. 1 for the particular example the operator A from [7], see Eq. (1.2). In this example, we have a(x) = 1 + α 2 α 3 sin x , b(x) = 1 − (α 1 + α 2 ) cos x α 3 , c(x) = 0 .
The results of our computations are shown in Figure 1. We see that if the parameter α 3 is small, then the surface tension is not strong enough to overcome the gravity and the model is unstable with the number of the unstable eigenvalues growing as the parameter α 3 decreases.