Orthonormal eigenfunction expansions for sixth-order boundary value problems

Sixth-order boundary value problems (BVPs) arise in thin-film flows with a surface that has elastic bending resistance. To solve such problems, we first derive a complete set of odd and even orthonormal eigenfunctions — resembling trigonometric sines and cosines, as well as the so-called “beam” functions. These functions intrinsically satisfy boundary conditions (BCs) of relevance to thin-film flows, since they are the solutions of a self-adjoint sixth-order Sturm–Liouville BVP with the same BCs. Next, we propose a Galerkin spectral approach for sixth-order problems; namely the sought function as well as all its derivatives and terms appearing in the differential equation are expanded into an infinite series with respect to the derived complete orthonormal (CON) set of eigenfunctions. The unknown coefficients in the series expansion are determined by solving the algebraic system derived by taking successive inner products with each member of the CON set of eigenfunctions. The proposed method and its convergence are demonstrated by solving two model sixth-order BVPs.


Introduction
Sixth-order parabolic equations arise in a number of physical contexts.For example, early work by King [1] on isolation oxidation of silicon led to a model featuring a sixth-order parabolic equation due to the effect of diffusion and reaction on substrate deformation during the oxidation process.Importantly, King [1] observed that the resulting partial differential equation is a higher-order version of the well-known nonlinear, degenerate parabolic equations of second and fourth order that describe the height of a thin liquid film during gravity-and surface-tension-driven spreading on a rigid surface.At vastly different scales, the flow of magma under earthen layers leads to the formation of domes called laccoliths.The model used the describe the formation and spread of these geological features [2,3,4,5] also led to sixth-order parabolic equations, similar to the one derived by King [1], but usually in axisymmetric polar coordinates.
In the last two decades, there has been significant interest in fundamental questions about the behavior of solutions of nonlinear sixth-order parabolic equations, motivated by problems of thin liquid films spreading underneath an elastic membrane [6,7,8,9,10,11].With the advent of cheap and rapid fabrication methods and improved imaging modalities, a growing area of application for these models and studies has been actuation in microfluidics [12,13], including impact mitigation [14].In these contexts, it is often of interest to understand infinitesimal perturbations of the film and their stability, leading to sixth-order eigenvalue problems.
In the mathematics literature, general properties of sixth-order eigenvalue problems have been discussed by Greenberg and Marletta [15,16].However, the eigenfunctions are seldom explicitly constructed, the studies focusing instead on calculation of eigenvalues for problems arising in hydrodynamic stability [17] or vibrations of sandwich beams [18,19].We surmise that constructing the eigenfunctions of sixth-order eigenvalue problems would be enlightening, and that these functions may resemble the well-known "beam" functions [17,20,21] of the fourth-order eigenvalue problem.
Thus, in the present work, motivated by recent interest in sixth-order parabolic equations, we study a model initial-boundary-value problem for a sixth-order equation (section 2), construct the eigenfunctions of its associated eigenvalue problem (section 3), use the latter to propose a Galerkin spectral method [22,23] for the solution of the former (section 4), and demonstrate the approach on two elementary model problems (section 5).

Problem formulation: elastic-plated thin film 2.1. Problem definition and flow geometry
Here, we study a thin fluid film of equilibrium height h 0 , underneath an elastic interface in a closed trough of width 2ℓ, as shown in figure 1.The fluid is considered incompressible and Newtonian with (constant) density ρ f and (constant) dynamic viscosity µ f .The elastic sheet has (constant) in-plane tension T and (constant) bending resistance B. The sheet is idealized as an interface with no mass (thus, no inertia during its motion).Gravity is the only body force considered to act on the fluid, and it acts in the −y direction.The elastic interface is free to move vertically along the lateral walls at x = ±ℓ while maintaining a 90 • contact angle (non-wetting).
Figure 1: Schematic of an elastic sheet (with in-plane elastic tension T , out-of-plane bending rigidity B, and negligible mass), overlaying a viscous fluid (density ρ f and dynamic viscosity µ f ) of equilibrium height h 0 in a closed trough of axial width 2ℓ.A non-uniform height h(x, t) = h 0 + u(x, t) drives a horizontal flow x (x, y, t) under the elastic interface, and vice versa.Gravity acts in the −y-direction, and the gravitational acceleration is denoted by g.

Derivation of the governing thin-film equation
We study the dynamics of the thin film in the long-wave limit [24,25], also known as the lubrication approximation [26,27], such that h 0 ≪ ℓ as well as max x,t h(x, t) ≪ ℓ during the film's evolution.Under this assumption, it can be shown [7,8,9,10] that the dimensional governing equations for the viscous fluid flow under the elastic sheet reduce to: conservation of mass: x-momentum balance: kinematic interface condition: See also [6,28] for more in-depth discussions of the dynamic boundary condition at a fluid-elastic solid interface and its linearization in the long-wave limit.
Under the lubrication approximation, the flow is approximately unidirectional and y ≪ x .The dominant velocity component can be obtained from (2.1b) as having enforced the no-slip conditions x (x, 0, t) = x (x, h, t) = 0 at the bottom of the trough (y = 0) and along the elastic interface (y = h).Based on (2.2), the volumetric flux per unit width is 3) The vertical velocity can be obtained from the conservation of mass equation (2.1a) using x from (2.2) as having enforce the no-penetration condition y (x, 0, t) = 0 at the bottom of the trough (y = 0).Using (2.3) and (2.4), the kinematic condition (2.1e) becomes Next, fixing the pressure at the interface y = h(x, t) to be the "ambient pressure" p a = const., equations (2.1c) and (2.1d) allow us to solve for the pressure distribution in the thin fluid film: Finally, using the pressure (2.6) to eliminate ∂p/∂x from the flow rate expression (2.3), and substituting the result into (2.5),we obtain the (long-wave) thin film equation: Observe that (2.7) is a nonlinear partial differential equation (PDE), nominally of parabolic type, but with possible degeneracy at a "contact line" where h → 0, which leads to a wealth of possible contact-line behaviors [7,29] beyond the scope of this work.

Nondimensionalization
We introduce the dimensionless ("hat") variables via the transformations: We focus on the case in which the dominant force on the interface is the bending resistance, thus the choice of time scale is t c = 12µ f ℓ 6 /Bh 3 0 , as in [14,30].Substituting the variables from (2.8) into (2.7) and dropping the hats, we arrive at: (2.9) The following dimensionless numbers have arisen: Elastic Bond number: Tension number: The elastic Bond number B quantifies the relative importance of gravity's role in deforming the elastic interface relative to its bending resistance (see also discussion in [31]).The tension number T quantifies the relative importance of in-plane tension's role in deforming the elastic interface relative to its bending resistance (see also the discussions in [28,32] and [12,13]).We are interested in small perturbations to a flat film, h = 1, so we introduce h(x, t) = 1 + εu(x, t) with ε ≪ 1 into (2.9) and keep only terms at O(ε), to obtain Equation (2.11) is the linear sixth-order PDE of interest in the remainder of this work.

Boundary conditions
Equation (2.11) requires six boundary conditions (BCs).Some are discussed in [8] for a thin film problem and in [28, section 4.4.2] in the context of beams.According to the problem posed in section 2.1 and shown in figure 1, we are interested in a nonwetting, non-pinned film in a closed trough.A non-wetting film will maintain a 90 • contact at the wall, thus Since the film is not pinned (i.e., it is free), u(±1, t) is free and not necessarily = 0. Similarly, since the film is confined in a trough with walls, the walls provide shear force/resistance to the motion, so that (∂ 3 u/∂x 3 )| x=±1 is not necessarily = 0, but this in turn requires that there be no moment at the wall: Finally, we observe from the derivation of (2.7), that if the fluid flux q is to vanish at x = ±1 (closed trough) then, taking into account (2.12a) and (2.12b), we must additionally require that If we restrict ourselves to cases in which tension in negligible compared to bending (T → 0), then (2.12c)

The sixth-order Sturm-Liouville boundary value problem
In the context of the PDE (2.11) subject to the BCs (2.12), we consider the following sixth-order initial boundary value problem (IBVP) for u = u(x, t): where u 0 (x) and f (x, t) are given smooth functions.Although, on physical grounds, we took T → 0 to derive BC (2.12c), we have kept T in (3.1a) for the purposes of the mathematical discussion that follow.
To solve IBVP (3.1), we will employ a Galerkin method (section 4).To this end, we must construct a suitable set of eigenfunctions, and determine their properties, which is the subject of the remainder of this section.

Self-adjointness, orthogonality and completeness
The Sturm-Liouville eigenvalue problem (EVP) associated with IBVP (3.1) reads: Before we proceed with the derivation of the solutions (eigenfunctions) of EVP (3.2), we first discuss why they possess certain properties required for the development of our Galerkin method in section 4. The desired results follow as direct applications of the comprehensive general theory on eigenvalue problems found in the seminal book of Coddington and Levinson [33].To facilitate this discussion for the reader, we quote or paraphrase the relevant definitions and theorems from [33] accordingly.We begin with the definitions.
Definition 3.1.[Adjoint of a linear operator (from [33], chapter 3, section 6)] Consider the linear differential operator where a 0 , a 1 , . . ., a n are continuous (possibly complex-valued) functions of Then, the adjoint of L n is another linear operator, denoted L + n , also of order n given by where the bar ( ¯) denotes the complex conjugate.
Throughout this work we will use the notation Also, ∥ • ∥ will denote the associated L 2 norm unless explicitly specified otherwise.
Definition 3.2.[Self-adjoint eigenvalue problem (from [33], chapter 7, section 2)] Let L n be the n-th order operator given by where the p i (x) are (possibly complex-valued) functions of class C n− j , j = 0, 1, . . ., n on the closed interval [a, b] and p 0 (x) 0 on [a, b].Let where the M jk and N jk are constants.Denote the boundary conditions U j ϕ = 0, j = 1, . . ., n by U ϕ = 0.
The eigenvalue problem The first important consequence of self-adjointness is reflected in theorem 3.1.
[From [33], theorem 2.1 of chapter 7, section 2] A self-adjoint EVP has a real, enumerable set of eigenvalues with corresponding distinct and mutually-orthogonal eigenfunctions.
We now show the following proposition.The second important result regarding self-adjoint EVPs was shown in the expansion and completeness theorems, namely, theorem 4.1 and theorem 4.2 of chapter 7, section 4, of [33].We paraphrase this result as follows.
where the equality implies convergence in the L 2 norm.

Constructing the orthonormal eigenfunctions
We now solve EVP (3.2) and derive our complete orthonormal (CON) set of solutions.The characteristic equation corresponding to (3.2a) has roots Therefore, the general solution of the homogeneous version of (3.2a) can be written in the form . (3.14) We have chosen to express the solution in this way, following the example of the so-called "beam" functions [17,20,21], in order to have odd and even sets of eigenfunctions resembling trigonometric sines and cosines.Imposing BCs (3.2b) on each of the ψ s (x) and ψ c (x) separately yields homogeneous linear systems for the vectors c odd = (c 1 , c 5 , c 6 ) ⊤ and c even = (c 2 , c 3 , c 4 ) ⊤ , respectively.To find nontrivial solutions we set the determinants of the coefficient matrices equal to zero and arrive at the eigenvalue relations: Each solution λ c m of (3.15a) is an eigenvalue corresponding to an even eigenfunction ψ c m of EVP (3.2), whereas each solution λ s m of (3.15b) corresponds to an odd eigenfunction ψ s m .Eigenvalue relations (3.15) were solved numerically using a highly accurate numerical solver-Mathematica's FindRoot [34] with 16 digits of working precision.We also derived large-λ asymptotic formulas for the eigenvalues.Using (3.15a), we found that λ c m ∼ (m + 1/6)π as λ → ∞, whereas (3.15b) yields λ s m ∼ (m − 1/3)π as λ → ∞.The results are presented in table 1.As we can see, the asymptotic formulas are extremely accurate even for m = 6, with the formulas for λ c m and λ s m agreeing with the results obtained with FindRoot for 12 and 11 digits, respectively.Thus, in practice when implementing a Galerkin method, there will be no need to use numerical root-finding to determine λ c m and λ s m for m ≥ 7. What remains is to determine the constants c 2 , c 3 and c 4 for the set of even functions {ψ c m } and c 1 , c 5 and c 6 for the set of odd functions {ψ s m }, which were introduced in (3.14).These constants are found by imposing BCs (3.2b), utilizing the eigenvalue relations (3.15), and normalizing the eigenfunctions Table 1: "Even" and "odd" eigenvalues, evaluated by solving the transcendental equations (3.15) with Mathematica's [34] FindRoot, as well as the proposed asymptotic formulas.even odd with respect to the L 2 norm.After a lengthy calculation, using Mathematica's algebraic manipulation capabilities, we arrive at the expressions: where for n, m ∈ N.
It is important to note here that λ c 0 = 0 is an eigenvalue of EVP (3.2) with corresponding eigenfunction ψ c 0 (x) = 1.The notation is due to the fact that ψ(x) = 1 is an even function.The profiles of the first members of the two sets of eigenfunctions are presented in figure 2.
We conclude the this section by stating the following proposition, which forms the foundation of our spectral method.
with the series (3.17 The proposition follows from proposition 3.1, theorem 3.1 and proposition 3.2.It is clear that the only nonzero coefficients in the spectral expansion (3.17) of any even (or odd) function f ∈ L 2 [−1, 1] are the u c m (or u s m ).In this section n indicated the order of the EVP, and m was the index used for the countable infinity of eigenvalues and eigenfunctions, we have now restricted ourselves to a specific sixth-order EVP.Thus, in the following sections, we reuse n as the index for the countable infinity of eigenvalues and eigenfunctions to simplify the notation.

The Galerkin spectral method 4.1. Method formulation
Having established the necessary theoretical foundation and derived our CON basis functions, we proceed with the development of the proposed Galerkin spectral method [22,23] for IBVP (3.1).To solve the sixth-order PDE (3.1a) using a Galerkin method based on the results above, we expand the sought function (i.e., the solution of the IBVP) u(x, t) as in (3.17), allowing the expansion coefficients to be functions of time, and truncate the series at M terms: Note that the spectral expansion (4.1) will intrinsically satisfy BCs (3.1b), which is an important advantage of the classical Galerkin approach.
Next, we must express all the x-derivatives that appear in the PDE (3.1a) as linear combinations of the same basis functions: Note that, while β s n0 = γ s n0 = 0 by definition, it can be shown that β c n0 = 0 but γ c n0 0. The expressions for β c,s nm and γ c,s nm will be given in section 4.2 below.It follows that Next, substituting (4.1) and ( 4.3) into (3.1a),taking successive inner products with ψ c 0 (x) = 1, ψ c ℓ (x) M ℓ=1 and ψ s ℓ (x) M ℓ=1 , and using the orthogonality of the eigenfunctions, we obtain a semi-discrete (dynamical) system: where { f c 0 (t), f c m (t), f s m (t)} are the expansion coefficients of the known right-hand side f (x, t) of (3.1a), according to proposition 3.3.System (4.4) can be solved using a Crank-Nicolson-type approach for the time discretization, but its solution is the object of future work.

Expansion formulas
We first present the formula for expanding the second derivative of an even eigenfunction into a series of the even basis functions.After a lengthy calculation, using Mathematica's algebraic manipulation The decay rate of the odd coefficients (not plotted here) is |γ s 3m | ∼ m −3 , matching that of its even counterparts.This was confirmed for all the examined values of n and is expected from the expression (4.8) for γ s nm .For the model problems considered in section 5, the following expansion formulas for even powers of x are also required: where p = 2, 4, . . ., 12.The expressions for the coefficients χ {p} m for m ≥ 1 are given in the Appendix.The case m = 0 is not needed, as will become apparent through (5.4) and (5.7) below.

Model problems: Results and discussion
In this first foray into Galerkin methods for sixth-order-in-space parabolic equations, and the construction of the associated higher-order beam eigenfunctions, we would like to restrict to the steady case, such that ∂u/∂t = 0.In this context, we would like to demonstrate the spectral accuracy of the eigenfunction expansion on two simple model boundary-value problems (BVPs), which have simple exact solutions.

Model problem I
Model problem I is the BVP: (5.1c) The problem (5.1) admits the exact solution: u exact (x) = (x − 1) 6 (x + 1) 6 . (5. 2) The exact solution (5.2) as well as f (x) in (5.1a) are even functions, and BCs (5.1b) are symmetric.Thus, we need only use the even eigenfunctions ψ c 0 (x) = 1 and ψ c n (x) ∞ n=1 for the spectral expansion.Applying the Galerkin spectral method introduced in section 4 leads us to introduce the truncated series Now, we substitute (5.3) into (5.1a),take inner products with ψ c 0 (x) = 1 and ψ c m (x) M m=1 , and employ the orthogonality relations, as in the derivation of (4.4).First, the coefficient u c 0 is found as where δ nm is Kronecker's delta.
The exact solution (5.2) and the spectral approximation (5.3)-(5.5)using the sixth-order eigenfunctions are compared in figure 4(a), showing they agree identically (at least visually).As seen in figure 4(b), the spectral expansion is highly accurate with the maximum absolute error being ≤ 5 × 10 −13 .Despite the Gibbs effect near x = ±1, figure 4(c) shows that the overall convergence rate of the spectral expansion is O(n −8 ).This extremely rapid convergence is due to the very definition of our eigenfunctions, namely that they satisfy a sixth-order BVP.Put differently: when solving the algebraic system (5.5),we are essentially dividing with quantities involving (λ c n ) 6 , and it is thanks to this fact that we have the impressive convergence rate, overcoming the more moderate decay rate of the coefficients appearing in the right-hide side of (5.5).

Model problem II
Model problem II is the BVP: (5.6c)

Proposition 3 . 3 .
[Complete orthonormal set] The two sets of solutions ψ c m and ψ s m of Sturm-Liouville problem (3.2) given by (3.16) and supplemented by ψ c 0 (x) = 1 form a complete orthonormal set of