Practical Error Estimates for Reynolds' Lubrication Approximation and its Higher Order Corrections
aa r X i v : . [ m a t h . A P ] J un PRACTICAL ERROR ESTIMATES FOR REYNOLDS’ LUBRICATIONAPPROXIMATION AND ITS HIGHER ORDER CORRECTIONS ∗ JON WILKENING † Abstract.
Reynolds’ lubrication approximation is used extensively to study flows betweenmoving machine parts, in narrow channels, and in thin films. The solution of Reynolds’ equationmay be thought of as the zeroth order term in an expansion of the solution of the Stokes equations inpowers of the aspect ratio ε of the domain. In this paper, we show how to compute the terms in thisexpansion to arbitrary order on a two-dimensional, x -periodic domain and derive rigorous, a priorierror bounds for the difference between the exact solution and the truncated expansion solution.Unlike previous studies of this sort, the constants in our error bounds either are independent of thefunction h ( x ) describing the geometry or depend on h and its derivatives in an explicit, intuitiveway. Specifically, if the expansion is truncated at order 2 k , the error is O ( ε k +2 ), and h enters intothe error bound only through its first and third inverse moments R h ( x ) − m dx , m = 1 ,
3, and viathe max norms k ℓ ! h ℓ − ∂ ℓx h k ∞ , 1 ≤ ℓ ≤ k + 2. We validate our estimates by comparing with finiteelement solutions and present numerical evidence that suggests that even when h is real analytic andperiodic, the expansion solution forms an asymptotic series rather than a convergent series. Key words. incompressible flow, lubrication theory, asymptotic expansion, Stokes equations,thin domain, a priori error estimates
AMS subject classifications.
DOI.
1. Introduction.
Reynolds’ lubrication equation [22, 20, 16, 12] is used exten-sively in engineering applications to study flows between moving machine parts, e.g.,in journal bearings or computer disk drives. It is also used in microfluid and bio-fluid mechanics to model creeping flows through narrow channels and in thin films.Although there is a vast literature (including several textbooks) on viscous flows inthin geometries, the equations are normally derived either directly from physical ar-guments [16] or using formal asymptotic arguments [12]. This is acceptable in mostcircumstances as the original equations (Stokes or Navier–Stokes) have also been de-rived from physical considerations, and by now the lubrication equations have beenused frequently enough that one can draw on experience and intuition to determinewhether they will work well for a given problem.On the other hand, as soon as the geometry of interest develops (or approaches)a singularity, or if we wish to compute several terms in the asymptotic expansion ofthe solution in powers of the aspect ratio ε , we rapidly leave the space of problems forwhich we can use experience as a guide; thus, it would be helpful to have a rigorousproof of convergence to serve as a guide to identify the features of the geometry thatcould potentially invalidate the approximation. For example, in [25], Wilkening andHosoi used lubrication theory to study the optimal wave shapes that an animal suchas a gastropod should use as it propagates ripples along its muscular foot to crawlover a thin layer of viscous fluid. In certain limits of this constrained optimizationproblem, the optimal wave shape develops a kink or cusp in the vicinity of the regionclosest to the substrate, and there is a competing mechanism controlling the size of ∗ This work was supported in part by the Director, Office of Science, Advanced Scientific Com-puting Research, U.S. Department of Energy contract DE-AC02-05CH11231. † Department of Mathematics and Lawrence Berkeley National Laboratory, University of Califor-nia, Berkeley, CA 94720 ([email protected]).1
JON WILKENING the modeling error (singularity formation versus nearness to the substrate). We foundthat shape optimization within (zeroth order) lubrication theory drives the geometryout of the realm of applicability of the lubrication model; however, by computinghigher order corrections and monitoring the errors (using the results of this paper),we learned that cusp-like singularities are appropriately penalized by the full Stokesequations, yielding nonsingular optimal solutions; see [25] for further details.
In most of the following papers, the Stokes or Navier–Stokes equations are solved in a domain Ω ε bounded below by a flat substrate andabove by a curved boundary y = εh ( x ) in two dimensions, or z = εh ( x, y ) in threedimensions, where ε is a small parameter and the function h is fixed. These solutionsare then compared to the solution of Reynolds’ equation (or to a truncated expansionsolution of the Stokes or Navier–Stokes equations), and the error is shown to convergeto zero in the limit as ε → ε → ε =1 is held fixed and the equations contain thesmall parameter). Cimatti assumes h has four weak derivatives (whereas, we requireonly h ∈ C , ) and shows that for any compact set K ⊂ Ω,(1.1) k εu − ¯ u k L (Ω) ≤ Cε, max (cid:0) k ε p x − ¯ p x k L ( K ) , k ε p y k L ( K ) (cid:1) ≤ Cε / , where u is the x -component of velocity, p is the pressure, a bar denotes the solution ofReynolds’ equation, and C is independent of ε but depends on h in the first inequalityand on h and K in the second. The scaling here in not standard: he imposes theboundary condition εu ( x,
0) = ¯ u ( x,
0) = const, which accounts for the extra factorof ε in each of the left-hand sides of (1.1). There are a few problems with Cimatti’sanalysis, notably the dependence of C on L (the “arbitrary cutoff” used to make theunbounded domain bounded) and the fact that some of his arguments seem to require ε to be small in comparison to C − ; however, his basic approach is interesting andinspired much of the work that followed in this subject.In 1986, Bayada and Chambat [3] generalized Cimatti’s work to three dimensions.They analyze the Stokes equations directly rather than using a stream function for-mulation, assume less regularity of h (apparently only h ∈ C ), and state their resultsin terms of limits (i.e., the quantities u εi , ε∂ x u εi , ∂ y u εi , and p ε in the solution of theStokes equations converge in L to the corresponding quantities in the solution ofReynolds equations as ε → ε . He proved that if h ( x, y ) is smooth, then thereis a constant C depending on h , N , and the boundary conditions such that(1.2) (cid:13)(cid:13) u − u N (cid:13)(cid:13) H + (cid:13)(cid:13) p − ε − p N (cid:13)(cid:13) L ≤ Cε N − / , where ( u , p ) is the solution of the Navier–Stokes equations, u N and p N are the termsof the asymptotic expansion truncated at the N th order (including a boundary layerexpansion near the lateral edges of the thin domain), and the norms are taken on the RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION ε (rather than the rescaled domain Ω). As a corollary, if the expansionis computed with “superfluous” terms that are afterwards treated as remainders, heobtains the optimal estimate(1.3) (cid:13)(cid:13) u − u N (cid:13)(cid:13) L + ε / (cid:13)(cid:13)(cid:13)(cid:16) π ε / ∇ (cid:17) (cid:0) u − u N (cid:1)(cid:13)(cid:13)(cid:13) L + (cid:13)(cid:13) p − ε − p N (cid:13)(cid:13) L ≤ Cε N +1 . Nazarov’s paper is concise to the point of being impenetrable at times. We interpret π ε / ∇ = ( ∂ x , ∂ y , ε / ∂ z ), but this symbol was not defined and may actually be a vari-able coefficient operator that incorporates the boundary conditions in its definition.We are also unsure of the definition of p and p N , as we would have expected p − ε − p N to appear together.In a later paper [19], Nazarov studies the asymptotics of the solution of the Stokesequations in a domain in which two smooth surfaces meet at a point. This problem isalso studied in a recent paper of Ciuperca, Hafidi, and Jai in [9]. This singular limit isinteresting in that deriving even the first correction to the zeroth order approximationin the asymptotic expansion remains an open problem.Assemien, Bayada, and Chambat [2] have studied the important question of theeffect of inertia on the asymptotic behavior of a thin film flow, which can in manycases be significant, requiring that the Navier–Stokes equations be used in place ofthe Stokes equations as the underlying model for the asymptotic expansion. We alsomention that there is a large body of literature on the long-time behavior of solutionsof the Navier–Stokes equations on thin domains; see, e.g., [21, 17].In 2000, Duvnjak and Maru˘si´c-Paloka [11] showed how to rigorously analyze thelubrication approximation of the Navier–Stokes equations for a slipper bearing in acircular geometry. The focus of their paper is on formulating the problem in cylindricalcoordinates and showing how to adapt the zeroth order case of Nazarov’s proof tohandle the change of variables. Elrod’s pioneering 1960 paper [12] is also concernedwith the (formal) relationship between the Navier–Stokes equations and Reynolds’equation for this geometry. None of the studies described above showshow the constant C bounding the error depends on the function h ( x ) describing thegeometry. This is because most theorems of analysis give constants that depend onthe domain Ω, which is usually fixed. But in our case, the data h ( x ) of the problemactually specifies the domain; therefore, to obtain bounds that are independent of h ,one must avoid or modify standard arguments for flattening the boundary, etc., so asnot to lose track of h ( x ) in the analysis. Moreover, arguments based on the closedgraph theorem or Rellich’s compactness theorem must be avoided entirely, as thesealso depend on the geometry. This forces us to look for new ways to analyze oldproblems using tools that furnish explicit constants.In this paper, we consider only the two-dimensional, periodic Stokes equationswith a specific choice of boundary conditions, but we derive error estimates thatdepend on h in an explicit, intuitive way. Our main result is summarized in Theo-rem 4.11, which may be stated as follows: Let T = [0 , p be the periodic unit interval.If k ≥ h ∈ C k +1 , ( T ), 0 < h ≤ h ( x ) ≤ x ∈ T , and ε ≤ r / ε -weighted Sobolev norms) at order 2 k (keeping in mindthat only even powers of ε appear in these expansions) is bounded by(1.4) p I ( | V | + | V | ) " θ k εr k r I I ερ k r k (cid:19) k +2 , JON WILKENING where V and V are prescribed tangential velocities on the lower and upper boundariesof the domain(1.5) r k = max ≤ ℓ ≤ k +2 ((cid:13)(cid:13)(cid:13)(cid:13) ℓ ! h ℓ − ∂ ℓx h (cid:13)(cid:13)(cid:13)(cid:13) /ℓ ∞ )! − , I m = Z h ( x ) − m dx and ρ k , θ k are constants independent of h . The bound on pressure has another terminvolving h ; see (4.110) below.The constants in (1.4) have been divided into two types: those that are (1) givenin the problem statement or easily computable from h ; or (2) difficult to compute butuniversal (independent of h ). We show how to compute the constants in the lattercategory ( ρ k and θ k ) in section 4; see Table 4.4. The constants in the former category( r k and I m ) help us understand the competing mechanism of singularity formationversus proximity to the substrate: the curvature and higher derivatives are allowed todiverge as long as the gap size simultaneously approaches zero in such a way that thehomogeneous products ℓ ! h ℓ − ∂ ℓx h remain uniformly bounded. Although the factors √ I and p I /I in (1.4) also diverge in this limit, the norm of the exact solutiondiverges at a similar rate — so the relative error in the expansion solution truncatedat order 2 k is O ( ε k +2 ), with ρ k r k serving as an effective radius of convergence.The framework we have chosen for this paper is intended to be general enoughto cover many interesting applications (such as a crawling gastropod [25] or an “un-wrapped” slipper bearing) but simple enough to obtain explicit detailed estimatesthat reveal the dependence of the error on the geometry h ( x ). We also wanted todetermine whether there might exist geometries for which the asymptotic expansionyields a convergent series. Although we do not have a rigorous proof, the answerappears to be negative even for the simplest case of a real analytic function such as h ( x ) = + sin 2 πx , for which the r k in (1.5) are bounded away from zero. It ishoped that this work will serve as a useful first step toward obtaining similar errorestimates for three-dimensional problems that include more general boundary condi-tions, incorporate end effects near the lateral edges of the domain (which we avoid bystudying the periodic case), and include the effect of inertia or viscoelasticity. In section 2, we derive Reynolds’ lubrication approximation in itsprimitive and stream function formulations. In section 3, we show how to computesuccessive terms in an asymptotic expansion of the stream function. In section 3.2,we prove a structure theorem describing the dependence of these terms on h ( x ) andits derivatives.In section 4, we formulate the problem weakly and analyze the truncation errorequation using weighted Sobolev spaces and a uniform Poincar´e–Friedrichs argument.The first challenge is to find the right weighted norms on the lower and upper bound-aries (equivalent to H / (Γ ) and H / (Γ ) for fixed ε ) to yield manageable errorestimates in terms of h when we change variables to straighten out the boundaries.In section 4.4, we reduce the problem of bounding the truncation errors to that ofbounding the second and fourth derivatives of the two highest order terms retainedin the asymptotic expansion, namely, k ψ (2 k ) xx k and k h ψ (2 k − xxxx k . We then use thestructure theorem of section 3.2 to compute these norms in order to obtain the con-stants ρ k and θ k in (1.4) for 0 ≤ k ≤
25. In section 4.5, we show how to computethe error in velocity, vorticity, and pressure from that of the stream function. Thisrequires that we determine how the Babu˘ska–Brezzi inf-sup constant β depends on h ( x ); see [24]. RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION r k ρ k is within a factor of 3 ofoptimal for k = 5, k = 10, and perhaps all k ≥
5. These calculations also suggest thateven when h ( x ) is real analytic, the expansion solution is an asymptotic series ratherthan a convergent series. This is because the constants ρ k converge to zero as k → ∞ .Fortunately, ρ k initially increases and does not become smaller than ρ = 0 .
197 until2 k = 26, which is already outside of the practical range of k . Finally, in Appendix A,we present our numerical algorithm for computing the expansion solutions, which canbe performed symbolically using a computer algebra system such as Mathematica orin floating point arithmetic, e.g., in C ++ .
2. Reynolds’ approximation.
Consider the Stokes equations on a periodicdomain of width ¯ W bounded below by a flat wall moving with constant speed ¯ V andabove by an inextensible sheet moving with constant speed ¯ V along a fixed curveΓ ,ε = { (¯ x, ¯ h (¯ x )) : 0 ≤ ¯ x ≤ ¯ W } ; see Figure 2.1. A bar is used to distinguish a physicalvariable from its dimensionless counterpart. We nondimensionalize the variables bychoosing a characteristic speed ¯ U and height ¯ H for the problem, and set ¯ x = ¯ W x ,¯ y = ¯ Hy , ¯ h (¯ x ) = ¯ Hh ( x ), ¯ V i = ¯ U V i , (¯ u, ¯ v ) = ¯ u = ( ¯ U u, ¯ U ¯ H ¯ W v ), and ¯ p = ¯ µ ¯ U ¯ W ¯ H p .The stream function ψ , flux Q , and vorticity ω introduced below satisfy ¯ ψ = ¯ U ¯ Hψ ,¯ Q = ¯ U ¯ HQ , and ¯ ω = ¯ U ¯ H ω .We have in mind a situation where the aspect ratio ε = ¯ H/ ¯ W of the physicaldomain is small. By scaling the x - and y -axes differently, we map the problem ontoa nicer geometry, which introduces terms in the equations that vanish in the singularlimit ε →
0. Specifically, we wish to find x -periodic functions u, v, p defined on therescaled domain(2.1) Ω = { ( x, y ) : 0 ≤ x ≤ , < y < h ( x ) } such that(2.2) p x = ε u xx + u yy , p y = ε v xx + ε v yy , v y = − u x (in Ω)subject to periodic boundary conditions on the left and right sides of Ω and(2.3) ( u, v ) | Γ = ( g , , ( u, v ) | Γ = ( g , h x g ) ¯ x = 0 ¯ x = ¯ W − e e ¯u = ¯ V e ¯u = ¯ V t Γ ,ε Γ ,ε tn ¯ y = y = ¯ H Ω ε ¯ h (¯ x ) − ¯ µ ¯∆ ¯u + ¯ ∇ ¯ p = 0¯ ∇ · ¯u = 0 − e e Γ Γ u = V e u = V q h x ε h x t x = 0 x = 1 y = y = Ω tn h ( x ) − ∆ ε u + p x = 0 − ε ∆ ε v + p y = 0 ∇ · u = 0 Fig. 2.1 . Geometry commonly encountered in lubrication-type problems. Left: Physical coor-dinate system. Right: Dimensionless coordinate system ( ∆ ε = ε ∂ x + ∂ y ). JON WILKENING on the bottom and top boundaries. Here(2.4) g ( x ) = V , g ( x ) = V (cid:2) ε h ′ ( x ) (cid:3) − / , i.e., g ( x ) = V cos θ ( x ), where θ = arctan( εh x ) is the angle of the curve ¯ h (¯ x ) relativeto the horizontal. Reynolds’ lubrication approximation is obtained by setting ε = 0in the equations and solving(2.5) p x = u yy , p y = 0 , v y = − u x , u | Γ = ( V ; 0) , u | Γ = (1; h x ) V . If we write (2.2) in the form L ( u ; p ) = (0; 0; 0), where L = L (0) + ε L (2) + ε L (4) isgiven by(2.6) L = − ∂ y ∂ x ∂ y ∂ x ∂ y + ε − ∂ x − ∂ y
00 0 0 + ε − ∂ x
00 0 0 , then (2.5) is just the zeroth order system L (0) ( u ; p ) = (0; 0; 0) with zeroth orderboundary conditions (expanding g and g in (2.3) in powers of ε ). The equation for v decouples from the others, and we find that p is independent of y and(2.7) u ( x, y ) = (cid:18) y − h ( x ) y (cid:19) p x ( x ) + (cid:18) − yh ( x ) (cid:19) V + yh ( x ) V . Integrating from 0 to h and solving for p x , we obtain(2.8) p x = 6 h ( V + V ) − h Q, where Q = R h u ( x, y ) dy is the volume flux through any cross section of the fluid.( Q is constant since ∇ · u = 0 and u is tangent to Γ and Γ ). Since p is periodic, R p x dx = 0, and we find that(2.9) Q = V + V I I (cid:18) I m = Z h ( x ) − m dx (cid:19) . Substituting (2.9) and (2.8) into (2.7) and using v y = − u x , v ( x,
0) = 0, we obtain thesolution p x = 6( V + V ) h (cid:18) − I I h (cid:19) ,u = ( V + V ) (cid:18) I I h − (cid:19) (cid:18) yh − y h (cid:19) + (cid:16) − yh (cid:17) V + yh V , (2.10) v = ( V + V ) (cid:18) I I h − (cid:19) (cid:18) y h − y h (cid:19) h x + V y h h x . The vertical component v of the velocity field is customarily omitted from zerothorder lubrication theory as ¯ v = ε ¯ U v is O ( ε ) on the thin geometry Ω ε of Figure 2.1.We may also derive (2.10) using a stream function formulation of the problem.Our procedure for computing higher order corrections to the lubrication approxima-tion and our method for estimating the error of these expansion solutions are bothdone in the stream function formulation. Let us define(2.11) ∆ ε = ε ∂ x + ∂ y , u ε = (cid:18) uε v (cid:19) . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION L ( u ; p ) = ( F ; F ; 0) with boundary conditions (2.3), i.e.,(2.12) − ∆ ε u ε + ∇ p = F , ∇ · u = 0 , u | Γ = ( g ; 0) , u | Γ = ( g ; h x g ) . Since u is incompressible, there is a stream function ψ such that(2.13) u = ∇ × ψ = ( ψ y , − ψ x ) , ∇ × u ε = ε v x − u y = − ∆ ε ψ. It follows from (2.12) that ψ satisfies the rescaled biharmonic equation(2.14) ∆ ε ψ = ψ yyyy + 2 ε ψ xxyy + ε ψ xxxx = ∇ × F , with periodic boundary conditions in the x -direction and(2.15) ( ψ = 0 ψ y = g ) on Γ , ( ψ = Qψ y ( x, h ( x )) = g ) on Γ , where Q = R h (0)0 u (0 , y ) dy . Since p is periodic, R p x ( x, dx = 0, i.e.,(2.16) Z ψ yyy ( x,
0) + F ( x, dx = 0 . Conversely, suppose we are able to find a flux Q and a classical solution ψ of (2.14)and (2.15) such that (2.16) holds. Then we define u = ∇ × ψ and note that (2.14)implies ∇ × (∆ ε u ε + F ) ≡
0, i.e., the integral(2.17) p ( x, y ) = Z γ (∆ ε u ε + F ) · t ds, (cid:18) t = unit tangent vector alongpath γ joining (0 ,
0) to ( x, y ) (cid:19) is independent of the path γ . A canonical choice for γ is(2.18) p ( x, y ) = Z x (cid:2) ε u xx + u yy + F (cid:3) ( ξ, dξ + Z y (cid:2) ε v xx + ε v yy + F (cid:3) ( x, η ) dη. Condition (2.16) is equivalent to requiring p (1 ,
0) = p (0 , p (1 , y ) = p (0 , y ) for 0 ≤ y ≤ h (0), since the integrand of the second integral in(2.18) is periodic in x . By construction, the variables u , p satisfy (2.12), where theboundary condition on Γ follows from the fact that ψ x + h x ψ y = 0 there; hence,classical solutions of the rescaled biharmonic equation yield classical solutions of therescaled Stokes equations and vice versa. Reynolds’ approximation (2.10) is recoveredif F and ε are set to zero in (2.14)–(2.16) when solving for ψ and Q ; see section 3.1.
3. Higher order corrections.
In this section we show how to compute succes-sive terms in the formal expansion of the solution of the rescaled biharmonic equation(2.14) in powers of ε = ¯ H/ ¯ W . For this purpose, it is convenient to manipulate theequations assuming they are satisfied classically. Once we obtain formulas for thehigher order approximations, we will show (in section 4) that they satisfy a weakformulation of the problem that makes it possible to obtain error estimates. See [15]for background on perturbation methods in partial differential equations. JON WILKENING
Matching like powers of ε in the expansion(3.1) (cid:2) ∂ y + 2 ε ∂ x ∂ y + ε ∂ x (cid:3) h ψ (0) + ε ψ (2) + ε ψ (4) + · · · i = 0 , we obtain the recursion ψ (0) yyyy = 0 ,ψ (2) yyyy = − ψ (0) xxyy ,ψ (2 k ) yyyy = − ψ (2 k − xxyy − ψ (2 k − xxxx , k = 2 , , , . . . . (3.2)The boundary conditions (2.15) become Bψ (2 k ) = (cid:16) , g (2 k )0 , Q (2 k ) , g (2 k )1 (cid:17) , k = 0 , , , , . . . , (3.3)where Bψ = ( ψ | Γ , ψ y | Γ , ψ | Γ , ψ y | Γ ) and g ( x ), g ( x ) were defined in (2.4):(3.4) g (2 k )0 ( x ) = ( V , k = 0 , , k > , g (2 k )1 ( x ) = V (cid:18) − / k (cid:19) h ′ ( x ) k . Condition (2.16) (with F = 0) becomes(3.5) Z ψ (2 k ) yyy ( x, dx = 0 , k = 0 , , , . . . . If F were nonzero in (2.14) and depended on ε in such a way that ∇ × F had anexpansion in even powers of ε , we could incorporate these terms into (3.2) and (3.5)as well; however, we will assume F = except in section 4, where we consider thegeneral case only to derive error estimates for the F = case. Let us denote theright-hand side of (3.2) by f (2 k ) ( x, y ) for k ≥
0. The terms ψ (2 k ) , Q (2 k ) in (3.2) and(3.3) may be computed via(3.6) (cid:16) ψ (2 k ) , Q (2 k ) (cid:17) = G (cid:16) f (2 k ) , g (2 k )0 , g (2 k )1 (cid:17) , k = 0 , , , . . . , where G is defined by Algorithm 3.1 in Figure 3.1. In this algorithm, we solve ψ yyyy = f by integrating four times in the y -direction and then correct the bound-ary conditions with a cubic polynomial. The formula for Q in the algorithm maybe derived from the one for ψ as follows. As ψ ,yyy ( x,
0) = 0, the requirement that R ψ yyy ( x, dx = 0 is equivalent to the condition(3.7) 0 = 6 Z − Q + 2 ψ − ψ ,y h + g h + g hh dx. Solving for Q and using R h − dx = I gives the result.The formulas ( u, v ) = ( ψ y , − ψ x ), ω = ε v x − u y , p x = u yy + ε u xx , and p y = ε v yy + ε v xx allow us to compute the expansions of u , ω , and p in terms of ψ :(3.8) u (2 k ) = ψ (2 k ) y , v (2 k ) = − ψ (2 k ) x , k ≥ ,ω (0) = − ψ (0) yy , ω (2 k ) = − ψ (2 k − xx − ψ (2 k ) yy , k ≥ ,p (0) x = ψ (0) yyy , p (2 k ) x = ψ (2 k − xxy + ψ (2 k ) yyy , k ≥ ,p (0) y = 0 , p (2) y = − ψ (0) xyy , p (2 k ) y = − ψ (2 k − xxx − ψ (2 k − xyy , k ≥ ,p (2 k ) ( x, y ) = Z x p (2 k ) x ( ξ, dξ + Z y p (2 k ) y ( x, η ) dη, k ≥ . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION Algorithm 3.1. ( ψ, Q ) = G ( f, g , g ): ψ = V f (cid:18) V = Volterra operator: V f ( x, y ) = Z y f ( x, η ) dη (cid:19) Q = 12 I Z ψ ( x, h ( x )) h ( x ) + − ψ ,y ( x, h ( x )) + g + g ( x ) h ( x ) dxψ ( x, y ) = ψ ( x, y ) + ( g h ( x )) yh ( x ) + (3 Q − ψ ( x, h ( x )) + ψ ,y ( x, h ( x )) h ( x ) − g h ( x ) − g ( x ) h ( x )) y h ( x ) + ( − Q + 2 ψ ( x, h ( x )) − ψ ,y ( x, h ( x )) h ( x ) + g h ( x ) + g ( x ) h ( x )) y h ( x ) return ( ψ, Q ) Fig. 3.1 . Algorithm to solve ψ yyyy = f , Bψ = (0 , g , Q, g ) , R ψ yyy ( x, dx = 0 . Equation (3.2) implies that ∂ x p (2 k ) y = ∂ y p (2 k ) x for k ≥
0; hence, differentiating underthe integral sign in (3.8), we see that p (2 k ) x and p (2 k ) y actually are the partial derivativesof p (2 k ) . Finally, our choice of Q (2 k ) ensures R p (2 k ) x ( ξ, dξ = R ψ (2 k ) yyy ( x, dx = 0so that p (2 k ) is periodic.Using Algorithm 3.1 to evaluate ( ψ (0) , Q (0) ) = G (0 , V , V ) yields Q (0) = V + V I I , (3.9) ψ (0) = ( V h ) yh + (cid:16) Q (0) − (2 V + V ) h (cid:17) y h + (cid:16) − Q (0) + ( V + V ) h (cid:17) y h , (3.10)which agrees with Reynolds’ approximation (2.10) when u (0) , p (0) are computed from ψ (0) . To compute higher order terms in the expansion, we need to study the recursion(3.6) more closely to determine how h will enter into the formulas for Q (2 k ) and ψ (2 k ) . In this sec-tion, we show how the terms ψ (2 k ) and Q (2 k ) in the stream function expansion de-pend on h . The key result of this section is that these higher order corrections havea structure similar to the zeroth order formulas (3.9) and (3.10), but the coefficienton each y n h n now belongs to a more complicated polynomial algebra in the symbols V , V , h , the derivatives of h , and certain weighted averages of the products of h and its derivatives. We also present a concise representation for the correction termsusing matrices of rational numbers that are independent of any particular choice ofshape function h . By splitting the analysis into one part that holds universally andanother that depends on h in a simple way, we are able to derive useful error estimatesgoverning the expansion solution truncated at any order in section 4.Let P = Q [ h, h x , h xx , . . . ] denote the algebra of polynomials in h and its deriva-tives over the rationals. A typical element of P might be 3 + h h xx h xxx . In P , thegenerators h , h x , etc., are treated as symbols rather than functions. Thus, if h ( x )happens to equal 1 identically, the polynomials 1 − h and h x are nonzero in P eventhough they are mapped to zero when P is (noninjectively) embedded in C ∞ ( T ),the space of C ∞ functions on the periodic interval T = [0 , p . If h is not smooth,its derivatives can still be manipulated symbolically and various subspaces (involving0 JON WILKENING
Algorithm 3.2. ( basis generation ) for k = 0 , . . . , k Φ k = { t k } , ( or d k = 1) for j = 2 , . . . , k for k = j, . . . , k Φ k = Φ k ∪ t j Φ k − j , ( or d k = d k + d k − j ) return { Φ , . . . , Φ k } ⋆ • •• ••• ⋆ ⋆ ⋆ k j ⋆ Fig. 3.2 . Algorithm to find a canonical basis Φ k for each space H k in the range ≤ k ≤ k .Here t ↔ h x , . . . , t k ↔ k ! h k − ∂ kx h . terms with few enough derivatives) can still be embedded in actual function spacessuch as C k ( T ).For any monomial α = Ch i h i x h i xx · · · ∈ P with C = 0, we define its superdegree to be the number of derivatives present:(3.11) sdeg( α ) = i + 2 i + 3 i + · · · . If α ∈ P , we define its superdegree to be the maximal superdegree of any of its terms,and set sdeg(0) = −∞ . Since Q is a field, sdeg( αβ ) = sdeg( α ) + sdeg( β ) for any α, β ∈ P . We say that α is homogeneous of superdegree k if each of its terms hassuperdegree k .Let H ⊂ P denote the subalgebra generated by the set { h k − ∂ kx h : k ≥ } , i.e.,(3.12) H = Q [ { h x , hh xx , h h xxx , . . . } ] , and for k ≥
0, let H k ⊂ H denote the subspace(3.13) H k = { } ∪ { α ∈ H : α is homogeneous of superdegree k } . Note that H k is finite-dimensional for all k , and H = Q is the set of constantpolynomials. We will denote the dimension of H k by(3.14) d k = dim( H k ) . Given an integer k ≥
0, we can use Algorithm 3.2 in Figure 3.2 to construct acanonical basis Φ k = { ϕ ( k )1 , . . . , ϕ ( k ) d k } for each H k with k in the range 0 ≤ k ≤ k .For notational convenience, let t j stand for j ! h j − ∂ jx h . As the outer loop (on j )progresses, Φ k contains a basis for the subspace of H k that involves only the symbols t , . . . , t j . Let us denote these auxiliary sets by(3.15) Φ kj = { t i . . . t i j j : i + 2 i + · · · + ji j = k } , ≤ j ≤ k. Then Φ k = { t k } , Φ kk = Φ k , and Φ kj = Φ k,j − ∪ t j Φ k − j,j for 2 ≤ j ≤ k . In otherwords, Φ kj consists of Φ k,j − together with all products of the variables t , . . . , t j ofsuperdegree k that contain at least one power of t j . The first several Φ k and d k are RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION = { } , Φ = { h x } , Φ = (cid:26) h x , hh xx (cid:27) , Φ = (cid:26) h x , hh x h xx , h h xxx (cid:27) , (3.16)Φ = (cid:26) h x , hh x h xx , h h xx , h h x h xxx , h h xxxx (cid:27) , ( d , . . . , d ) = { , , , , , , , , , , } , d = 627 , d = 204226 . We have found empirically that the first 75000 terms satisfy ( √ k k +1 ) < d k < √ k k +1 .In fact, we have recently learned of the Hardy–Ramanujan formula(3.17) d k ∼ exp (cid:16) π p k/ (cid:17) k √ k → ∞ for the number of partitions of the integer k . Thus, rather than 13, the base is in fact e π √ / = 13 . h . In the following theorem, V H k is the tensor product of V and H k , where(3.18) V = { } ∪ { α ∈ Q [ V , V ] : α is homogeneous of degree 1 } is the space of rational linear combinations of V and V . Recall from (2.9) above that I m = R h ( x ) − m dx . Theorem 3.3.
The terms Q (2 k ) , ψ (2 k ) in the stream function expansion definedby the recursion (3.6) and Algorithm have the form (3.19) Q (2 k ) = I I a (2 k ) + k − X ℓ =0 Q (2 ℓ ) b (2 k − ℓ ) , ψ (2 k ) = I I α (2 k ) + k X ℓ =0 Q (2 ℓ ) β (2 k − ℓ ) , where (3.20) α (2 k ) ( x, y ) = k +3 X n =1 α (2 k ) n ( x ) y n h ( x ) n , β (2 k ) ( x, y ) = k +3 X n =1 β (2 k ) n ( x ) y n h ( x ) n , and (3.21) α (2 k ) n ∈ I I h V H k , β (2 k ) n ∈ H k . Moreover, a (2 k ) = I R α (2 k )3 ( x ) h ( x ) dx and b (2 k ) = I R β (2 k )3 ( x ) h ( x ) dx .Remark In addition to pinning down the way in which h appears in theformulas for the stream function expansion, this theorem allows us to represent ψ (2 k ) and Q (2 k ) using matrices of rational numbers. Explicitly, (3.20) and (3.21) hold iffthere are matrices A (2 k ) , B (2 k ) with entries in V and Q , respectively, with rowsindexed from 0 to (2 k + 3) and columns indexed from 1 to d k , and containing2 JON WILKENING only zeros in row 0, such that(3.22) α (2 k ) ( x, y ) = ( Y k ( x, y )) T A (2 k ) (cid:18) I I h ( x )Φ k ( x ) (cid:19) ,β (2 k ) ( x, y ) = ( Y k ( x, y )) T B (2 k ) Φ k ( x ) , where Y k = (1 , yh , . . . , ( yh ) k +3 ) T and Φ k = ( ϕ (2 k )1 , . . . , ϕ (2 k ) d k ) T are treated as columnvectors. The purpose of the zeroth row is to make it easy to convert to orthogonalpolynomials in y/h if desired. The final statement of the theorem asserts that theformulas for a (2 k ) and b (2 k ) are also encoded in the matrices A (2 k ) and B (2 k ) . If weadopt Matlab notation and denote row i of A (2 k ) by A (2 k ) ( i, :), then(3.23) a (2 k ) = 12 A (2 k ) (3 , :) E (2 k )2 , b (2 k ) = 12 B (2 k ) (3 , :) E (2 k )3 , where(3.24) E (2 k ) m = (cid:16) E (2 k ) m, , . . . , E (2 k ) m,d k (cid:17) T = 1 I m Z Φ k ( x ) h ( x ) m dx. Note that E (2 k ) m,j is the weighted average of ϕ (2 k ) j with weight function I − m h − m . Forexample, E (0) m = (1), E (2) m = ( I m R h x h m dx, I m R hh xx h m dx ) T , etc.; see (3.16) above. Example
We can now represent Q (0) and ψ (0) in (3.9) and (3.10) by(3.25) a (0) = V + V , A (0) = V − + V − , B (0) = − . The second order terms Q (2) and ψ (2) involve these as well as a (2) = 12 (cid:2) V (cid:0) , (cid:1) + V (cid:0) , − (cid:1)(cid:3) E (2)2 , b (2) = 12 (cid:0) − , − (cid:1) E (2)3 , (3.26) A (2) = V − /
15 2 / /
15 2 / / − / − / / + V − /
30 7 / / − / / − / − / / , B (2) = / − / − / − / − / − / . For k ≥ A (2 k ) and B (2 k ) are both (2 k + 4) × d k matrices with rows 0 and 1containing only zeros. These matrices are universal: the shape function h enters intothe formulas only through Y k , Φ k , and E (2 k ) m in (3.22) and (3.23). In Appendix A,we show how to compute A (2 k ) and B (2 k ) directly from the lower order matrices A (2 ℓ ) and B (2 ℓ ) with 0 ≤ ℓ < k . Proof of Theorem Q (0) and ψ (0) have thedesired form. Suppose k ≥
1, and the theorem holds for 0 ≤ k < k . We must showthat it is also true for k = k . By (3.6),(3.27) (cid:16) ψ (2 k ) , Q (2 k ) (cid:17) = G (cid:16) − ψ (2 k − xxyy , , g (2 k )1 (cid:17) , k = 1 ,G (cid:16) − ψ (2 k − xxyy − ψ (2 k − xxxx , , g (2 k )1 (cid:17) , k ≥ . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION ψ ( − should be replaced by zero. The first step of Algorithm 3.1 is to compute ψ (2 k )0 .Using the induction hypothesis, we obtain ψ (2 k )0 = I I (cid:16) − V α (2 k − xxyy − V α (2 k − xxxx (cid:17) (3.28) + k − X ℓ =0 Q (2 ℓ ) (cid:16) − V β (2 k − ℓ − xxyy (cid:17) + k − X ℓ =0 Q (2 ℓ ) (cid:16) − V β (2 k − ℓ − xxxx (cid:17) . The upper limit of the last sum can be replaced by k −
1, since we interpret β ( − as zero. We would like to rewrite this in the form(3.29) ψ (2 k )0 = I I k +3 X n =4 α (2 k ) n ( x ) y n h n ! + k − X ℓ =0 Q (2 ℓ ) k − ℓ +3 X n =4 β (2 k − ℓ ) n ( x ) y n h n ! . If we use the induction hypothesis and substitute (3.20) into (3.28), the operator V ∂ y annihilates a single power of y and antidifferentiates higher powers of y twice.Similarly, V antidifferentiates all powers of y four times. Thus, for k = k and4 ≤ n ≤ k + 3, we should define(3.30) α (2 k ) n ( x ) = − h n ∂ x (cid:16) α (2 k − n − h − n +2 (cid:17) n ( n −
1) + − h n ∂ x (cid:16) α (2 k − n − h − n +4 (cid:17) n ( n − n − n − β (2 k ) n in terms of β (2 k − n − and β (2 k − n − . The second termshould be omitted when k = 1 or n = 4, and is zero when n = 5. As part of theinduction hypothesis, we may assume that (3.30) and its β version hold for 1 ≤ k < k as well, so that each term in the sum over ℓ in (3.28) also has the form described in(3.29). Note that for n ≥ ϕ ( x ),(3.31) ∂ x ( h − n ϕ ) = h − ( n +1) ( h∂ x − nh x ) ϕ. By Lemmas 3.6 and 3.7 below, h∂ x and multiplication by h x both map H k to H k +1 for all k ≥
0. Thus h n ∂ x (cid:16) α (2 k − n − h − n +2 (cid:17) (3.32) = h [ h∂ x − ( n − h x ][ h∂ x − ( n − h x ] (cid:16) h − α (2 k − n − (cid:17) ∈ I I V h H k ,h n ∂ x (cid:16) β (2 k − n − h − n +2 (cid:17) = [ h∂ x − ( n − h x ][ h∂ x − ( n − h x ] (cid:16) β (2 k − n − (cid:17) ∈ H k , with similar formulas for h n ∂ x ( α (2 k − n − h − n +4 ) and h n ∂ x ( β (2 k − n − h − n +4 ). We con-clude that α (2 k ) n and β (2 k ) n have the form claimed in (3.21) when k = k and 4 ≤ n ≤ k +3. Finally, we obtain Q (2 k ) and ψ (2 k ) from ψ (2 k )0 in (3.29) using Algorithm 3.1.4 JON WILKENING
They satisfy (3.19) and (3.20) if we set k = k and define α (2 k )1 = 0, β (2 k )1 = 0,(3.33) α (2 k )2 ( x ) = k +3 X n =4 ( n − α (2 k ) n ( x ) − V (cid:18) − / k (cid:19) h kx ,β (2 k )2 ( x ) = k +3 X n =4 ( n − β (2 k ) n ( x ) ,α (2 k )3 ( x ) = k +3 X n =4 (2 − n ) α (2 k ) n ( x ) + V (cid:18) − / k (cid:19) h kx ,β (2 k )3 ( x ) = k +3 X n =4 (2 − n ) β (2 k ) n ( x ) ,a (2 k ) = I R α (2 k )3 ( x ) h ( x ) dx , and b (2 k ) = I R β (2 k )3 ( x ) h ( x ) dx . As part of the inductionhypothesis, we may assume (3.33) also holds for 1 ≤ k < k . The factors of n in(3.33) are due to the terms ± ψ ,y ( x, h ) h in the formula for ψ (2 k ) in Algorithm 3.1.The terms 3 Q (2 k ) y h and − Q (2 k ) y h in the formula for ψ (2 k ) are accounted for in(3.19) by extending the upper limit of the sum over ℓ from k − k and noting that β (0) ( x, y ) = 3 y h − y h . Thus, ψ (2 k ) and Q (2 k ) have the desired form, and α (2 k ) n , β (2 k ) n belong to the appropriate spaces, as claimed.To complete this proof, we need two simple lemmas about the spaces H k (whichalso serve as the foundation for our numerical algorithm described in Appendix A). Lemma 3.6. If k ≥ and ϕ ∈ H k , then h x ϕ ∈ H k +1 .Proof . This follows easily from the definition of H k in (3.13). Lemma 3.7. If k ≥ and ϕ ∈ H k , then h∂ x ϕ ∈ H k +1 .Proof . If k = 0, then h∂ x ϕ = 0 ∈ H k +1 . Suppose k ≥
1, and the result holds for k < k . Let ϕ ∈ H k be a monomial. Then there is a k ∈ { , . . . , k } and a monomial β ∈ H k − k such that ϕ = ( h k − ∂ kx h ) β . But then(3.34) h∂ x ϕ = ( k − h x ϕ + (cid:0) h k ∂ k +1 x h (cid:1) β + (cid:0) h k − ∂ kx h (cid:1) ( h∂ x β ) . Evidently, all three terms belong to H k +1 , the third due to the induction hypothesis.This result can now be applied term by term for any polynomial ϕ ∈ H k .
4. Error analysis.
To estimate the error of the expansion of ψ and Q throughorder 2 k , we show that the truncation error satisfies a weak form of the rescaled bihar-monic equation (2.14) with data ( F , g , g ) of order ε k +2 . We also prove a uniformcoercivity result for the family of bilinear forms involved in the weak formulation,which allows us to bound the truncation error in terms of the data.Throughout this section, we will treat Ω and T = [0 , p as C ∞ manifolds byidentifying the points(4.1) Ω : (0 , y ) ∼ (1 , y ) , < y < h (0) ,T : 0 ∼ C k (Ω) or C k ( T ) is understood to have k continuous periodic derivatives; ∂ Ω =Γ ∪ Γ ; ∂T = ∅ ; the support of a function φ ∈ C kc (Ω) vanishes near Γ and Γ butnot necessarily at x = 0 and x = 1; and the Sobolev spaces H k (Ω) and H k (Ω) are thecompletions of C k (Ω) and C kc (Ω) in the k · k k norm and thus contain only x -periodicfunctions with appropriate smoothness at x = 0 , RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION An interest-ing difference between the biharmonic equation and the Poisson equation is that theboundary conditions in the latter are completely specified in the problem statement,whereas one of them (the flux Q ) in the former problem must be determined as partof the solution. The integral condition (2.16), which uniquely determines Q , mustalso be reformulated weakly, since it involves more than two derivatives of ψ . Thiscan be done [14] by slightly enlarging the space of test functions to include functionsthat are constant along Γ (rather than equal to 0 there). To this end, we define(4.2) Ψ = n φ ∈ H (Ω) : ( φ, ∂ y φ ) | Γ = (0 , , ( φ, ∂ y φ ) | Γ = (const , o . For φ , ψ in H (Ω), we define the bilinear form(4.3) a ε ( ψ, φ ) = Z Ω ψ yy φ yy + 2 ε ψ xy φ xy + ε ψ xx φ xx dA = a (0) ( ψ, φ ) + ε a (2) ( ψ, φ ) + ε a (4) ( ψ, φ ) . To obtain estimates that hold uniformly in ε , it will be useful to work with theweighted norms and seminorms(4.4) k ψ k = Z Ω ψ dA, | ψ | ,ε = Z Ω ψ y + ( εψ x ) dA, | ψ | ,ε = a ε ( ψ, ψ ) , k ψ k ,ε = q k ψ k + | ψ | ,ε , k ψ k ,ε = q k ψ k + | ψ | ,ε + | ψ | ,ε . For fixed ε , these norms are equivalent to the usual Sobolev norms in which ε is setto 1 in these expressions. We use x to parametrize functions defined on Γ or Γ anddefine the weighted boundary norm(4.5) k g k / ,ε = ∞ X k = −∞ (cid:2) πkε ) (cid:3) / | c k | , c k = Z g ( x ) e − πikx dx. We equip the dual spaces Ψ ′ and H − (Ω) = [ H (Ω) ] ′ with the weighted norms(4.6) k l k − ,ε = sup k ψ k ,ε =1 |h l, ψ i| , k F k − ,ε = sup k u k ,ε + k εv k ,ε =1 |h F , ( u, v ) i| . Since k ψ k ,ε ≥ k ψ y k ,ε + k εψ x k ,ε , the linear functional h l, ψ i = h F , ∇ × ψ i on Ψsatisfies k l k − ,ε ≤ k F k − ,ε . Definition 4.1 ( weak solutions ). Suppose (4.7) h ∈ C , ( T ) , F ∈ H − (Ω) , g ∈ H / (Γ ) , g ∈ H / (Γ ) . We say that ( ψ, Q ) ∈ H (Ω) × R is a weak solution of (2.14)–(2.16) if (4.8) a ε ( ψ, φ ) = h F , ∇ × φ i for all φ ∈ Ψ and the boundary conditions (4.9) Bψ = (0 , g , Q, g ) hold in the trace sense, where Bψ := ( ψ | Γ , ψ y | Γ , ψ | Γ , ψ y | Γ ) . Proposition 4.2.
Every classical solution is a weak solution. JON WILKENING
Proof . We assume ψ ∈ C (Ω) and (2.14)–(2.16) hold classically; this requiresadditional regularity for F , g , g , of course. If we multiply (2.14) by a test function φ ∈ C (Ω) ∩ Ψ and use the identity χ ( ∇ × v ) = ∇ × ( χ v ) + ( ∇ × χ ) · v , we obtain0 = Z Ω φ ( − ∆ ε ψ + ∇ × F ) dA (4.10)= Z Ω ( ∇ × [ φ (∆ ε u ε + F )] + ( ∇ × φ ) · [∆ ε u ε + F ]) dA = Z Γ − Γ φ (∆ ε u ε + F ) · t ds + Z Ω [( ∇ × φ ) ε · ( ∇ × ∆ ε ψ ) + ( ∇ × φ ) · F ] dA = Z Γ − Γ [ φ ( · · · ) − (∆ ε ψ )( ∇ × φ ) ε ] · t ds + Z Ω [ − (∆ ε φ )(∆ ε ψ ) + ( ∇ × φ ) · F ] dA, where ( ∇ × φ ) ε = ( φ y , − ε φ x ) and the curves Γ and Γ are both oriented from leftto right as in Figure 2.1. The conditions(4.11) φ | Γ = 0 , ∂ y φ | Γ = 0 ,φ | Γ = const , ∂ y φ | Γ = 0 , ensure that the boundary terms are zero: the first boundary term is equal to(4.12) ( φ (cid:12)(cid:12) Γ )[ p (1 , h (1)) − p (0 , h (0))] = 0(with p as in (2.17), where it was shown to be periodic), and the second is zero since ∇ × φ = 0 on Γ and Γ . One more integration by parts gives R Ω (∆ ε φ )(∆ ε ψ ) dA = a ε ( ψ, φ ), so (4.8) holds. Since C (Ω) ∩ Ψ is dense in Ψ and both sides of (4.8) arebounded linear functionals of φ ∈ Ψ, this formula holds for all φ ∈ Ψ. The following two theorems are the key to obtainingerror estimates for the expansion solutions of section 3.
Theorem 4.3.
The bilinear form a ε ( · , · ) is coercive on Ψ (uniformly in ε ) withrespect to the weighted norm k·k ,ε , i.e., there is a constant α > such that α k ψ k ,ε ≤ a ε ( ψ, ψ ) for all ε > , ψ ∈ Ψ .Proof . Without loss of generality, we may assume the characteristic height ¯ H ofthe domain was chosen so that 0 < h ( x ) ≤ ≤ x ≤
1. We now use a standardPoincar´e–Friedrichs argument [5]. Suppose ψ ∈ C (Ω) ∩ Ψ. Then | ψ ( x, y ) | = (cid:12)(cid:12)(cid:12)(cid:12)Z y ψ y ( x, η ) dη (cid:12)(cid:12)(cid:12)(cid:12) ≤ y Z h ( x )0 | ψ y ( x, η ) | dη, (4.13) | ψ ( x, y ) | = (cid:12)(cid:12)(cid:12)(cid:12)Z y ( y − η ) ψ yy ( x, η ) dη (cid:12)(cid:12)(cid:12)(cid:12) ≤ y Z h ( x )0 | ψ yy ( x, η ) | dη. (4.14)Integrating over Ω,(4.15) k ψ k ≤ h Z Ω | ψ y | dA ≤ | ψ | ,ε , k ψ k ≤ h Z Ω | ψ yy | dA ≤ | ψ | ,ε . Repeating this argument on the derivatives of ψ yields(4.16) | ψ | ,ε = k ψ y k + k εψ x k ≤ (cid:0) | ψ y | ,ε + | εψ x | ,ε (cid:1) = 12 | ψ | ,ε = 12 a ε ( ψ, ψ ) RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION k ψ k ,ε ≤ a ε ( ψ, ψ ). Since C (Ω) ∩ Ψ is dense in Ψ, we conclude that(12 / k ψ k ,ε ≤ a ε ( ψ, ψ ) for all ψ ∈ Ψ as claimed.
Theorem 4.4.
A weak solution ψ of the boundary value problem (2.14)–(2.16) exists and is unique. Moreover, the following estimate holds: k ψ k ,ε ≤ k F k − ,ε (4.17) +
72 + 860 ε k h x k ∞ + ε k h x k ∞ + 4 ε (cid:13)(cid:13)(cid:13)(cid:13) hh xx (cid:13)(cid:13)(cid:13)(cid:13) ∞ / (cid:18)(cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε + (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε (cid:19) . In particular, if ε ≤ r with r − = max( k h x k ∞ , k hh xx k / ∞ ) , then (4.18) k ψ k ,ε ≤ k F k − ,ε + 15 (cid:18)(cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε + (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε (cid:19) . Proof . We begin by constructing a function ψ ∈ H (Ω) that satisfies the bound-ary conditions (2.15) with Q = 0. First, we map the domain Ω to the x -periodic unitsquare R = T × (0 ,
1) via the transformation(4.19) e ψ ( x, y ) = h ( x ) − / ψ ( x, h ( x ) y ) , ≤ x ≤ , < y < . We include h − / here to avoid powers of h − in (4.27), where h = min ≤ x ≤ h ( x ).We require e ψ ( x,
0) = 0, e ψ ,y ( x,
0) = h ( x ) − / g ( x ), e ψ ( x,
1) = 0, and e ψ ,y ( x,
1) = h ( x ) − / g ( x ). To construct such a function, we define ζ ∈ C ( R ) via ζ ( y ) = , y ≤ − y + 2 y + y , − ≤ y ≤ y − y + y , ≤ y ≤ , ≤ y −1 −0.5 0 0.5 1−0.200.2 ζ (y) slope = 1 and set(4.20) e ψ ( x, y ) = ∞ X k = −∞ (cid:18) c k ζ ( h k i y ) h k i + d k ζ ( h k i ( y − h k i (cid:19) e πikx , where(4.21) h k i = (cid:2) πkε ) (cid:3) / , [ c k , d k ] = Z [ g , g ]( x ) h ( x ) − / e − πikx dx. The value and slope of ζ at y = 0 and | y | ≥ e ψ satisfies the desiredboundary conditions. Assume for the moment that each d k is zero (i.e., g ≡ S be the strip T × R . We may use (4.20) to define e ψ on all of S and take its Fouriertransform(4.22) (cid:16) e ψ (cid:17) ∧ ( k, η ) = Z Z ∞−∞ e ψ ( x, y ) e − πi ( kx + ηy ) dy dx = c k ˆ ζ ( η/ h k i ) h k i . JON WILKENING
Since ζ is antisymmetric and supported on [ − , (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) ,ε,R = (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) ,ε,S ≤ (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) ,S + 2 (cid:12)(cid:12)(cid:12) e ψ (cid:12)(cid:12)(cid:12) ,ε,S + (cid:12)(cid:12)(cid:12) e ψ (cid:12)(cid:12)(cid:12) ,ε,S = X k Z ∞−∞ (cid:2) πkε ) + (2 πη ) (cid:3) (cid:12)(cid:12)(cid:12)(cid:12)(cid:16) e ψ (cid:17) ∧ ( k, η ) (cid:12)(cid:12)(cid:12)(cid:12) dη = X k h k i| c k | Z ∞−∞ (cid:2) πt ) (cid:3) (cid:12)(cid:12)(cid:12) ˆ ζ ( t ) (cid:12)(cid:12)(cid:12) dt = 898105 (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε . A similar argument works if we assume g ≡
0, but g
0. Thus, on R , we have(4.24) (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) ,ε ≤ s (cid:18)(cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε + (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε (cid:19) . Next we use the formula ψ ( x, y ) = h ( x ) / e ψ ( x, yh ( x ) ) to obtain(4.25) ψ ,y = h / e ψ ,y , ψ ,yy = h − / e ψ ,yy ,ψ ,x = h / e ψ ,x − yh − / h x e ψ ,y + 32 h / h x e ψ ,ψ ,xy = h / e ψ ,xy − yh − / h x e ψ ,yy + 12 h − / h x e ψ ,y ,ψ ,xx = h / e ψ ,xx + 3 h / h x e ψ ,x − yh − / h x e ψ ,xy + 32 h / h xx e ψ + y h − / h x e ψ ,yy − yh − / h x e ψ ,y − yh − / h xx e ψ ,y + 34 h − / h x e ψ . Using Lemma 4.5 below and 0 < h ( x ) ≤
1, we find that Z Ω ψ dA = Z Z ψ ( x, h ( x ) y ) h ( x ) dx dy = Z R h e ψ dA ≤ Z R e ψ dA, (4.26) Z Ω ψ ,y dA ≤ Z R e ψ ,y dA, Z Ω ψ ,yy dA = Z R e ψ ,yy dA, Z Ω ε ψ ,x dA ≤ Z R ε e ψ ,x dA + 4 ε k h x k ∞ (cid:20)Z R e ψ ,y dA + 94 Z R e ψ dA (cid:21) , Z Ω ε ψ ,xy dA ≤ Z R ε e ψ ,xy dA + 4 ε k h x k ∞ (cid:20) Z R e ψ ,yy dA + 12 Z R e ψ ,y dA (cid:21) , Z Ω ε ψ ,xx dA ≤ Z R ε e ψ ,xx dA + 30 ε k h x k ∞ Z R ε e ψ ,x dA + 30 ε k h x k ∞ Z R ε e ψ ,xy dA + 30 ε k h h xx k ∞ Z R e ψ dA + 30 ε k h x k ∞ Z R e ψ ,yy dA + 30 ε k h x k ∞ Z R e ψ ,y dA + 30 ε k h h xx k ∞ Z R e ψ ,y dA + 30 ε k h x k ∞ Z R e ψ dA. Note that the inverse powers of h in (4.25) are canceled when we change variablesfrom Ω to R by the factors of h that arise from the substitution y → hy and from the RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION k ψ k ,ε ≤ (cid:0) / ε k h x k ∞ + 30 ε k h x k ∞ + 30 ε k hh xx k ∞ (cid:1) / (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) ,ε . Finally, we correct ψ by a function in Ψ to obtain the weak solution ψ , which mustsatisfy(4.28) ψ − ψ ∈ Ψ , a ε ( ψ − ψ , φ ) = h l, φ i := h F , ∇× φ i− a ε ( ψ , φ ) for all φ ∈ Ψ . Since l is a bounded linear functional on Ψ, the Lax–Milgram theorem implies exis-tence and uniqueness of the solution ψ of (4.28) and gives the error bound(4.29) k ψ − ψ k ,ε ≤ α − k l k ,ε ≤ k F k − ,ε + k ψ k ,ε ) . Combining this with (4.24) and (4.27) and using the triangle inequality gives (4.17),where we note that ( ) ( ) ≤
72 and 30( ) ( ) ≤ Lemma 4.5.
For any a , . . . , a ∈ R , ( a + a + a ) ≤ a + 4 a + 4 a , (4.30)( a + · · · + a ) ≤ a + 103 a + 15 a + 403 a + 30 (cid:0) a + a + a (cid:1) + 1603 a . Proof . In general, given positive real numbers γ , . . . , γ n such that P n γ − j ≤ a ∈ R n , we have ( P n a j ) ≤ P n γ j a j . This is a consequence of theCauchy–Schwarz inequality:(4.31) X j a j = X j (cid:16) γ − / j (cid:17) (cid:16) γ / j a j (cid:17) ≤ X j γ − j X j γ j a j . One readily checks that + + = 1 and ( + · · · + ) = ≤ In section 3, we showed how to constructsuccessive terms in the stream function expansion by solving the recursion (3.2)–(3.5).Theorem 3.3 guarantees that derivatives of h higher than 2 k do not appear in theformulas for ψ (0) , . . . , ψ (2 k ) ; hence, if h ∈ C k ( T ), these functions satisfy (3.2)–(3.5)in the classical sense (with k replaced by ℓ and running from 0 to k instead of 0 to ∞ ).Thus, if h ∈ C k +4 , ψ (2 k ) approx = ψ (0) + ε ψ (2) + · · · + ε k ψ (2 k ) satisfies(4.32) ∆ ε ψ (2 k ) approx = ε k +2 (cid:16) ψ (2 k ) xxyy + ψ (2 k − xxxx (cid:17) + ε k +4 ψ (2 k ) xxxx . The truncation error ψ (2 k ) err = ψ exact − ψ (2 k ) approx then satisfies ∆ ε ψ err = − ∆ ε ψ approx with O ( ε k +2 ) boundary data. Since the right-hand side of (4.32) and the boundary dataare known in terms of h , we are able to estimate the size of ψ (2 k ) err using Theorem 4.4above. However, to use this theorem, we need to formulate (4.32) weakly.0 JON WILKENING
We begin by showing that the ψ (2 ℓ ) satisfy a weak version of the recursion (3.2).Suppose k ≥ h ∈ C k ( T ). Let φ ∈ Ψ and denote the constant value of φ on Γ by q . We multiply both sides of (3.2) by φ and integrate by parts using Z Ω φ ψ (2 ℓ ) yyyy dA = Z Ω φ yy ψ (2 ℓ ) yy dA + q Z ψ (2 ℓ ) yyy ( x, h ( x )) dx, (4.33)2 Z Ω φ ψ (2 ℓ ) xxyy dA = 2 Z Ω φ xy ψ (2 ℓ ) xy dA + q Z h ψ (2 ℓ ) xxy ( x, h ( x )) − ψ (2 ℓ ) xyy ( x, h ( x )) h x i dx, Z Ω φ ψ (2 ℓ ) xxxx dA = Z Ω φ xx ψ (2 ℓ ) xx dA − q Z ψ (2 ℓ ) xxx ( x, h ( x )) h x dx to obtain the recursion a (0) (cid:16) ψ (0) , φ (cid:17) = 0 ,a (0) (cid:16) ψ (2) , φ (cid:17) = − a (2) (cid:16) ψ (0) , φ (cid:17) , (4.34) a (0) (cid:16) ψ (2 ℓ ) , φ (cid:17) = − a (2) (cid:16) ψ (2 ℓ − , φ (cid:17) − a (4) (cid:16) ψ (2 ℓ − , φ (cid:17) , ℓ = 2 , , . . . , k. By (3.8), the boundary terms in (4.33) combine to form(4.35) q Z h p (2 ℓ ) x ( x, h ( x )) + p (2 ℓ ) y ( x, h ( x )) h x i dx = 0 , ≤ ℓ ≤ k, when substituted into (3.2). Other boundary terms do not arise in (4.33), since φ = 0on Γ and φ x = φ y = 0 on Γ and Γ .Now suppose h ∈ C k +1 , ( T ), i.e., h ( x ) has 2 k + 1 continuous periodic derivativesand ∂ k +1 x h is Lipschitz continuous so that ∂ k +2 x h ∈ L ∞ ( T ). Let ( ψ exact , Q exact ) bethe weak solution of(4.36) ∆ ε ψ = 0 , Bψ = (0 , g , Q, g ) , with g , g given in (2.4), and define the truncation errors and approximate solutions(4.37) ψ (2 k ) err = ψ exact − ψ (2 k ) approx , ψ (2 k ) approx = ψ (0) + ε ψ (2) + · · · + ε k ψ (2 k ) ,Q (2 k ) err = Q exact − Q (2 k ) approx , Q (2 k ) approx = Q (0) + ε Q (2) + · · · + ε k Q (2 k ) . Since ψ (0) , . . . , ψ (2 k ) satisfy (4.34) and (3.3) while a ε ( ψ exact , φ ) = 0 for every φ ∈ Ψ,we may expand a ε ( ψ (2 k ) err , φ ) in powers of ε to obtain the truncation error equation(4.38) a ε (cid:16) ψ (2 k ) err , φ (cid:17) = − ε k +2 h F k , ∇ × φ i , φ ∈ Ψ ,Bψ (2 k ) err = (cid:16) , , Q (2 k ) err , ε k +2 γ k (cid:17) , where γ k = ε − k − (cid:16) g − h g (0)1 + ε g (2)1 + · · · + ε k g (2 k )1 i(cid:17) (4.39)and (cid:28) F k , (cid:18) φ y − φ x (cid:19)(cid:29) = ( a (2) (cid:0) ψ (0) , φ (cid:1) + ε a (4) (cid:0) ψ (0) , φ (cid:1) , k = 0 ,a (2) (cid:0) ψ (2 k ) , φ (cid:1) + ε a (4) (cid:0) ψ (2 k ) , φ (cid:1) + a (4) (cid:0) ψ (2 k − , φ (cid:1) , k ≥ . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION F k ∈ H − (Ω) that have this action on the subspace(4.40) V = {∇ × φ : φ ∈ Ψ } = (cid:8) ( u, v ) ∈ H (Ω) : u x + v y = 0 (cid:9) . For example, F k = ε − k − [ ∇ p (2 k ) approx − (∆ ε u (2 k ) approx ) ε ] satisfies ∆ ε ψ (2 k ) approx = ε k +2 ∇× F k classically and, using (3.8), may be shown to have following action on H (Ω) :(4.41) h F k , ( u ; v ) i = Z Ω (cid:16) ψ (0) xx (cid:17) (cid:0) u y − ε v x (cid:1) − (cid:16) ε − ψ (0) xy (cid:17) ( εv y ) dA, k = 0 , h F k , ( u ; v ) i = Z Ω (cid:16) ψ (2 k ) xx (cid:17) (cid:0) u y − ε v x (cid:1) − (cid:16) ε − ψ (2 k ) xy (cid:17) ( εv y ) − (cid:16) ε − ψ (2 k − xx (cid:17) (cid:0) ε v x (cid:1) dA, k ≥ . This choice is suboptimal because the terms ε − ψ (2 k ) xy and ε − ψ (2 k − xx diverge as ε → ε with v y and ε with v x due to the definition (4.6) of k F k k − ,ε . Instead,we will use the following functional, which agrees with (4.41) on V : h F k , ( u ; v ) i = Z Ω (cid:16) ψ (0) xx (cid:17) (cid:0) u y − ε v x (cid:1) dA, k = 0 , h F k , ( u ; v ) i = Z Ω (cid:16) ψ (2 k ) xx (cid:17) (cid:0) u y − ε v x (cid:1) + h ψ (2 k − xxxx e φ [ u ] dA, k ≥ , (4.42) e φ [ u ]( x, y ) := − Z h ( x ) y u ( x, η ) h ( x ) dη = Z h ( x ) y ( η − y ) u y ( x, η ) h ( x ) dη. Note that if φ ∈ Ψ and q = φ | Γ , then e φ [ φ y ] = ( φ − q ) h − . The purpose of the h − here is to be able to include an h with ψ (2 k − xxxx in the error estimates below (toproperly consolidate terms). Another alternative to (4.42) that would lead to similarestimates below is h F k , ( u ; v ) i = R Ω [ − ε ψ (2 k ) xx v x − ψ (2 k +2) yy u y ] dA . Let us assume from now on that ε ≤ r with r =max( k h x k ∞ , k hh xx k / ∞ ) − . Then by (4.38) and Theorem 4.4, we have(4.43) (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε ≤ ε k +2 (cid:18) k F k k − ,ε + 15 (cid:13)(cid:13)(cid:13) h − / γ k (cid:13)(cid:13)(cid:13) / ,ε (cid:19) . It remains to bound the norms of F k and γ k in terms of h . From (4.42), we have(4.44) (cid:12)(cid:12)(cid:12) e φ [ u ]( x, y ) (cid:12)(cid:12)(cid:12) ≤ Z h ( x ) y ( η − y ) h ( x ) dη ! Z h ( x )0 | u y ( x, η ) | dη ! , where the first integral is ( h − y ) h and the second is independent of y . Hence(4.45) (cid:13)(cid:13)(cid:13) e φ [ u ] (cid:13)(cid:13)(cid:13) = Z Z h ( x )0 (cid:12)(cid:12)(cid:12) e φ [ u ]( x, y ) (cid:12)(cid:12)(cid:12) dy dx ≤ k u y k . From (4.42), we then have |h F k , ( u ; v ) i| ≤ (cid:18) a + b √ (cid:19) k u y k + a (cid:13)(cid:13) ε v x (cid:13)(cid:13) ≤ (cid:18) a + 4 ab √
12 + b (cid:19) / (cid:0) k u k ,ε + k εv k ,ε (cid:1) / , (4.46)2 JON WILKENING where a = k ψ (2 k ) xx k and b = k h ψ (2 k − xxxx k . Using ab √ ≤ a + b , we find that(4.47) k F k k − ,ε ≤ r a + 720 b ≤ a + 35 b, k F k k − ,ε ≤ a + b. Finally, by (4.43), we obtain(4.48) (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε ≤ ε k +2 (cid:18) (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13) h ψ (2 k − xxxx (cid:13)(cid:13)(cid:13) + 15 (cid:13)(cid:13)(cid:13) h − / γ k (cid:13)(cid:13)(cid:13) / ,ε (cid:19) , where the fourth derivative term should be omitted when k = 0. In section 4.5, thefollowing bound will also prove useful:(4.49) (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + ε k +2 (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) ≤ ε k +2 (cid:18) (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13) h ψ (2 k − xxxx (cid:13)(cid:13)(cid:13) + 15 (cid:13)(cid:13)(cid:13) γ k h / (cid:13)(cid:13)(cid:13) / ,ε (cid:19) . The truncation error in the flux expansion satisfies(4.50) Q (2 k ) err = ψ (2 k ) err ( x, h ( x )) = Z h ( x )0 ( h ( x ) − η ) ∂ ψ (2 k ) err ∂y ( x, η ) dη for any x ∈ T . Using estimates similar to (4.44) and (4.45), we find that(4.51) (cid:12)(cid:12)(cid:12) Q (2 k ) err (cid:12)(cid:12)(cid:12) ≤ √ (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε . Thus, we have reduced the problem of bounding the truncation errors to that ofchecking the norms of three quantities that can be computed explicitly in closedform.We begin by attacking the boundary term k h − / γ k k / ,ε . This norm was definedin (4.5) above. Recall that(4.52) γ k = V (cid:2) ε h x (cid:3) − / − k X ℓ =0 (cid:18) − / ℓ (cid:19) (cid:0) ε h x (cid:1) ℓ ! ε − k − , where (cid:0) − / (cid:1) = 1 and for ℓ ≥ (cid:18) − / ℓ (cid:19) = ( − / − / · · · ( − [ ℓ − / · · · ( ℓ ) = ( − ℓ · · · · · ℓ − ℓ . Taking the logarithm of this product and its inverse, one may show that(4.54) 1 √ ℓ + 1 ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:18) − / ℓ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ≤ √ ℓ + 1 , ℓ = 0 , , , . . . . Since h ∈ C k +1 , ( T ) ⊂ C , ( T ), we know γ k is at least Lipschitz continuous on T and so belongs to H ( T ). As a result, (cid:13)(cid:13)(cid:13) γ k h / (cid:13)(cid:13)(cid:13) / ,ε ≤ (cid:13)(cid:13)(cid:13) γ k h / (cid:13)(cid:13)(cid:13) ,ε = Z h − γ k + ε (cid:18) − h − / h x γ k + h − / γ k,x (cid:19) dx ≤ Z h − γ k + 54 ε h − h x γ k + 54 ε h − γ k,x dx. (4.55) RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION ε k h x k ∞ ≤ / <
1, the binomial expansion of [1+ ε h x ] − / converges uniformlyon T = [0 , p . As the terms in this expansion alternate in sign, the error in truncatingthe series is smaller (pointwise) than the first omitted term. Therefore, for each x ∈ T ,(4.56) | γ k ( x ) | ≤ V (cid:12)(cid:12)(cid:12)(cid:12)(cid:18) − / k + 1 (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) h x ( x ) k +2 ≤ V √ k + 4 h x ( x ) k +2 . Since ( − ℓ ) (cid:0) − / ℓ (cid:1) = (cid:0) − / ℓ − (cid:1) , by differentiating (4.52) we obtain(4.57) γ k,x = V (cid:2) ε h x (cid:3) − / − k − X ℓ =0 (cid:18) − / ℓ (cid:19) (cid:0) ε h x (cid:1) ℓ ! (cid:0) − ε h x h xx (cid:1) ε − k − . The terms in the expansion of [1 + ε h x ] − / also alternate in sign, so(4.58) | γ k,x ( x ) | ≤ V (cid:12)(cid:12)(cid:12)(cid:12)(cid:18) − / k (cid:19) h x ( x ) k +1 h xx ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ V k + 1) √ k + 4 (cid:12)(cid:12) h x ( x ) k +1 h xx ( x ) (cid:12)(cid:12) . Combining (4.55), (4.56), and (4.58) and using ( k +1) k +4 ≤ + k , we find that (cid:13)(cid:13)(cid:13) γ k h / (cid:13)(cid:13)(cid:13) / ,ε ≤ V k + 4 Z h k +4 x h + 54 ε h x h " h k +4 x + 16( k + 1) h kx (cid:18) hh xx (cid:19) dx (4.59) ≤ V I e E (2 k +2)1 , + V ε k h x k ∞ (cid:20) I e E (2 k +2)3 , + (cid:18) k (cid:19) I e E (2 k +2)3 , (cid:21) ≤ V I (cid:20)
14 + ε k h x k ∞ I /I (cid:18) k (cid:19)(cid:21) max ( m,j ) ∈{ (1 , , (3 , , (3 , } e E (2 k +2) m,j , where(4.60) I m = Z h ( x ) m dx, e E (2 k ) m,j = 1 I m Z (cid:16) ϕ (2 k ) j ( x ) (cid:17) h ( x ) m dx. Recall that Φ k = { ϕ (2 k )1 , . . . , ϕ (2 k ) d k } = { h kx , hh k − x h xx , . . . , k )! h k − ∂ kx h } is abasis for the space H k defined in (3.13). e E (2 k ) m,j is the square of the 2-norm (or secondmoment) of ϕ (2 k ) j ( x ) with respect to the probability measure I − m h ( x ) − m dx , whereas E (2 k ) m,j in (3.24) is the expected value of ϕ (2 k ) j ( x ) with respect to this measure. Using(4.59) in (4.48) gives a bound on the error caused by failing to satisfy the boundaryconditions in the stream function expansion. It is perhaps surprising that this boundcan be expressed in terms of three simple integrals involving h and its derivatives. Remark I and I are dimensionless quantities in (4.59) — if h and x stillcarried dimensions of length, an extra length scale (e.g., h max , which is currently set to1) would need to be included in the Sobolev norms to allow k ψ k , | ψ | ,ε , and | ψ | ,ε tobe added together; this length scale would also appear in (4.59) to nondimensionalize I /I in the denominator.Finally, we estimate k ψ (2 k ) xx k and k h ψ (2 k − xxxx k in (4.48) in terms of h . It will beuseful in our analysis to define(4.61) r k = max ≤ ℓ ≤ k +2 ((cid:13)(cid:13)(cid:13)(cid:13) ℓ ! h ℓ − ∂ ℓx h (cid:13)(cid:13)(cid:13)(cid:13) /ℓ ∞ )! − JON WILKENING so that for 0 ≤ ℓ ≤ k + 1, 1 ≤ j ≤ d ℓ , and m = 1 , ,
3, we have(4.62) (cid:12)(cid:12)(cid:12) ϕ (2 ℓ ) j ( x ) (cid:12)(cid:12)(cid:12) ≤ r − ℓk , (cid:12)(cid:12)(cid:12) E (2 ℓ ) m,j (cid:12)(cid:12)(cid:12) ≤ r − ℓk , ≤ e E (2 ℓ ) m,j ≤ r − ℓk . If h is real analytic as well as periodic, a standard contour integral argument showsthat there is an r > k ∂ kx h k ∞ ≤ k ! r − k for all k ≥
0. Such an r serves asa common lower bound for r k in (4.61) for all k ≥
0. If h is a constant function, theresults below hold with r k = ∞ and r − k = 0 (i.e., the lubrication approximation isexact).Our first task will be to bound the growth of the terms Q (2 k ) in the expansion ofthe flux. By Theorem 3.3 and Remark 3.4, there are rational matrices A (2 k )0 , A (2 k )1 , B (2 k ) with rows indexed from 0 to 2 k + 3 and columns indexed from 1 to d k suchthat(4.63) Q (2 k ) = I I a (2 k ) + k − X ℓ =0 Q (2 ℓ ) b (2 k − ℓ ) , k ≥ , where a (2 k ) = [ V A (2 k )0 (3 , :) + V A (2 k )1 (3 , :)] E (2 k )2 and b (2 k ) = B (2 k ) (3 , :) E (2 k )3 . SeeExample 3.5 above for a reminder of how this works. We now use (4.62) togetherwith the fact that | v · w | ≤ k v k k w k ∞ for v, w ∈ R d to conclude that(4.64) (cid:12)(cid:12)(cid:12) a (2 ℓ ) (cid:12)(cid:12)(cid:12) ≤ (cid:16) | V | κ (2 ℓ )0 + | V | κ (2 ℓ )1 (cid:17) r − ℓk , (cid:12)(cid:12)(cid:12) b (2 ℓ ) (cid:12)(cid:12)(cid:12) ≤ κ (2 ℓ )2 r − ℓk , ≤ ℓ ≤ k, where(4.65) κ (2 k ) i = 12 d k X j =1 (cid:12)(cid:12)(cid:12) A (2 k ) i (3 , j ) (cid:12)(cid:12)(cid:12) , i = 0 , , κ (2 k )2 = 12 d k X j =1 (cid:12)(cid:12)(cid:12) B (2 k ) (3 , j ) (cid:12)(cid:12)(cid:12) , k ≥ . It follows from (4.63) that if we increase κ (2 k )0 , κ (2 k )1 via the loop for k = 1 , , , . . .κ (2 k ) i = κ (2 k ) i + k − X ℓ =0 κ (2 ℓ ) i κ (2 k − ℓ )2 , ( i = 0 , , (4.66)then(4.67) (cid:12)(cid:12)(cid:12) Q (2 ℓ ) (cid:12)(cid:12)(cid:12) ≤ I I (cid:16) | V | κ (2 ℓ )0 + | V | κ (2 ℓ )1 (cid:17) r − ℓk , ( k ≥ , ≤ ℓ ≤ k ) . The constants κ (2 k )0 , κ (2 k )1 , κ (2 k )2 do not depend on h and may be computed once andfor all along with the matrices A (2 k ) i and B (2 k ) ; see Table 4.1.Now that the terms Q (2 ℓ ) have been bounded, we are ready to estimate k ψ (2 k ) xx k and k h ψ (2 k − xxxx k . Recall that(4.68) ψ (2 k ) ( x, y ) = I I α (2 k ) ( x, y ) + k X ℓ =0 Q (2 ℓ ) β (2 k − ℓ ) ( x, y ) , where α (2 k ) ( x, y ) and β (2 k ) ( x, y ) have the matrix representations (3.22). A slightmodification of the proof of Theorem 3.3 shows that there are also matrices ˙ A (2 k )0 , RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION Table 4.1 κ (2 k ) i before and after loop (4.66) . k κ (2 k )0 before κ (2 k )1 before κ (2 k )2 κ (2 k )0 after κ (2 k )1 after0 5 . × − . × − . × +00 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × +00 . × +00 . × +00 . × +00 . × +00 . × +00 . × +00 . × +01 . × +01 . × +01 . × +01 . × +01 . × +01 . × +02 . × +02 . × +02 . × +02 . × +02 . × +03 . × +03 . × +03 . × +03 . × +04 . × +04 . × +05 . × +05 . × +05 . × +05 . × +06 . × +06 . × +06 . × +06 . × +06 . × +08 . × +08 . × +08 . × +08 . × +08
10 6 . × +09 . × +10 . × +09 . × +10 . × +10 ˙ A (2 k )1 , ˙ B (2 k ) and ¨ A (2 k )0 , ¨ A (2 k )1 , ¨ B (2 k ) of dimension (2 k + 4) × d k +2 and (2 k + 4) × d k +4 ,respectively, such that I I α (2 k ) xx ( x, y ) = h ( x ) − ( Y k ( x, y )) T h V ˙ A (2 k )0 + V ˙ A (2 k )1 i Φ k +2 ( x ) ,β (2 k ) xx ( x, y ) = h ( x ) − ( Y k ( x, y )) T ˙ B (2 k ) Φ k +2 ( x ) ,I I α (2 k ) xxxx ( x, y ) = h ( x ) − ( Y k ( x, y )) T h V ¨ A (2 k )0 + V ¨ A (2 k )1 i Φ k +4 ( x ) ,β (2 k ) xxxx ( x, y ) = h ( x ) − ( Y k ( x, y )) T ¨ B (2 k ) Φ k +4 ( x ) , (4.69)where Y k = (1 , yh , . . . , ( yh ) k +3 ) T . We can achieve significantly sharper estimates of k ψ (2 k ) xx k and k h ψ (2 k − xxxx k by expressing the dependence of ψ on y using orthogonalpolynomials. Let(4.70) e Y k = (cid:16) e P (cid:16) yh (cid:17) , e P (cid:16) yh (cid:17) , . . . , e P k +3 (cid:16) yh (cid:17)(cid:17) T be the vector of shifted Legendre polynomials [1], which satisfy(4.71) Z e P m ( x ) e P n ( x ) dx = δ mn n + 1 , m, n ≥ . The first several are(4.72) e P = 1 , e P = 2 x − , e P = 6 x − x + 1 , e P = 20 x − x + 12 x − . The well-known recurrence [1](4.73) e P n ( x ) = 2 n − n (2 x − e P n − ( x ) − n − n e P n − ( x ) , n ≥ , can be used to construct a nested family of lower triangular matrices R k of dimension(2 k + 4) × (2 k + 4) with indices starting at zero such that(4.74) e Y k = R k Y k , Y T k = e Y T k R − T k . For example,(4.75) R = − − − −
30 20 , R − T = / / /
40 1 / / /
200 0 1 / /
40 0 0 1 / . JON WILKENING
The entries of R − T k are nonnegative and have unit column sums for all k ≥
0. Nextwe renormalize the shifted Legendre polynomials and define˙ P n ( x, y ) = h ( x ) − / √ n + 1 e P n ( y/h ( x )) , (4.76)˙ Y k = (cid:16) ˙ P , . . . , ˙ P k +3 (cid:17) T = h − / D k e Y k , D k = diag (cid:16) √ , √ , . . . , √ k + 7 (cid:17) so that (4.69) becomes I I α (2 k ) xx = (cid:16) ˙ Y k ( x, y ) (cid:17) T D − k R − T k h V ˙ A (2 k )0 + V ˙ A (2 k )1 i (cid:16) h ( x ) − / Φ k +2 ( x ) (cid:17) , (4.77) β (2 k ) xx = (cid:16) ˙ Y k ( x, y ) (cid:17) T D − k R − T k ˙ B (2 k ) (cid:16) h ( x ) − / Φ k +2 ( x ) (cid:17) ,I I h α (2 k ) xxxx = (cid:16) ˙ Y k ( x, y ) (cid:17) T D − k R − T k h V ¨ A (2 k )0 + V ¨ A (2 k )1 i (cid:16) h ( x ) − / Φ k +4 ( x ) (cid:17) ,h β (2 k ) xxxx = (cid:16) ˙ Y k ( x, y ) (cid:17) T D − k R − T k ¨ B (2 k ) (cid:16) h ( x ) − / Φ k +4 ( x ) (cid:17) . Example
From (3.9) and (3.10), we have Q (0) = V + V I I ,ψ (0) xx = (cid:20)(cid:18) V + V (cid:19) (cid:18) − h x + 2 hh xx h (cid:19) + Q (0) (cid:18) h x − hh xx h (cid:19)(cid:21) y h (4.78) + (cid:20) ( V + V ) (cid:18) h x − hh xx h (cid:19) + Q (0) (cid:18) − h x + 6 hh xx h (cid:19)(cid:21) y h , which has the form described in (4.68) and (4.69) with(4.79) h ˙ A (0)0 , ˙ A (0)1 , ˙ B (0)0 i = − − − − − −
24 12 and the form described in (4.68) and (4.77) with D − R − T h ˙ A (0)0 , ˙ A (0)1 , ˙ B (0)0 i =
16 26 56 − − √ √ √ − √ − √ − √ √ − √ √ − √ − √ √ √ − √ √ − √ − √ √ , where we recall that Φ ( x ) = ( h x , hh xx ) T . The matrices ¨ A (0)0 , ¨ A (0)1 , ¨ B (0) representing ψ (0) xxxx are each 4 × ( x ) was given in (3.16).To compute k ψ (2 k ) xx k and k h ψ (2 k − xxxx k , we note that each of the expressions in(4.77) is of the form P k +3 n =0 ˙ P n ( x, y ) w n ( x ), where w ( x ) = Sz ( x ), S is a constantmatrix, z ( x ) = h − m Φ k +2 j ( x ), j = 1 or 2, and m = 1 or 3. For fixed x , we have R h ( x )0 ˙ P m ( x, y ) ˙ P n ( x, y ) dy = δ mn . It follows that(4.80) Z Z h ( x )0 X n ˙ P n ( x, y ) w n ( x ) ! dy dx = Z X n w n ( x ) dx. RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION P n w n ( x ) = k w ( x ) k ≤ k S k F k z ( x ) k , where k · k F and k · k are theFrobenius and 2-norms of a matrix and vector, respectively. Integrating, we have R k z ( x ) k dx = I m k e E (2 k +2 j ) m k . Since k e E (2 ℓ ) m k ≤ d ℓ r − ℓk for 0 ≤ ℓ ≤ k + 1, we define(4.81) K (2 k ) i = p d k +2 (cid:13)(cid:13)(cid:13) D − k R − T k ˙ A (2 k ) i (cid:13)(cid:13)(cid:13) F , i = 0 , , k ≥ ,K (2 k )2 = p d k +2 (cid:13)(cid:13)(cid:13) D − k R − T k ˙ B (2 k ) (cid:13)(cid:13)(cid:13) F , k ≥ , e K (2 k ) i = p d k +4 (cid:13)(cid:13)(cid:13) D − k R − T k ¨ A (2 k ) i (cid:13)(cid:13)(cid:13) F , i = 0 , , k ≥ , e K (2 k )2 = p d k +4 (cid:13)(cid:13)(cid:13) D − k R − T k ¨ B (2 k ) (cid:13)(cid:13)(cid:13) F , k ≥ , so that(4.82) (cid:13)(cid:13)(cid:13)(cid:13) I I α (2 ℓ ) xx (cid:13)(cid:13)(cid:13)(cid:13) ≤ p I (cid:16) | V | K (2 ℓ )0 + | V | K (2 ℓ )1 (cid:17) r − ℓ − k , ≤ ℓ ≤ k, (cid:13)(cid:13)(cid:13) β (2 ℓ ) xx (cid:13)(cid:13)(cid:13) ≤ p I K (2 ℓ )2 r − ℓ − k , ≤ ℓ ≤ k, (cid:13)(cid:13)(cid:13)(cid:13) I I h α (2 ℓ ) xxxx (cid:13)(cid:13)(cid:13)(cid:13) ≤ p I (cid:16) | V | e K (2 ℓ )0 + | V | e K (2 ℓ )1 (cid:17) r − ℓ − k , ≤ ℓ ≤ k − , (cid:13)(cid:13)(cid:13) h β (2 ℓ ) xxxx (cid:13)(cid:13)(cid:13) ≤ p I e K (2 ℓ )2 r − ℓ − k , ≤ ℓ ≤ k − . From the bound (4.67) on | Q (2 ℓ ) | and the formula (4.68) for ψ (2 k ) in terms of α (2 k ) and β (2 k − ℓ ) , we see that after increasing K (2 k )0 , K (2 k )1 , e K (2 k )0 , e K (2 k )1 via(4.83) K (2 k ) i = K (2 k ) i + k X ℓ =0 κ (2 ℓ ) i K (2 k − ℓ )2 , i = 0 , , k ≥ , e K (2 k ) i = e K (2 k ) i + k X ℓ =0 κ (2 ℓ ) i e K (2 k − ℓ )2 , i = 0 , , k ≥ , we have(4.84) (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) ≤ p I (cid:16) | V | K (2 k )0 + | V | K (2 k )1 (cid:17) r − k − k , k ≥ , (cid:13)(cid:13)(cid:13) h ψ (2 k − xxxx (cid:13)(cid:13)(cid:13) ≤ p I (cid:16) | V | e K (2 k − + | V | e K (2 k − (cid:17) r − k − k , k ≥ . In each term Q (2 ℓ ) β (2 k − ℓ ) , we have used I /I I ≤ I / √ I by √ I . The constants K (2 k ) i , e K (2 k ) i do notdepend on h and may be computed once and for all along with the constants κ (2 k ) i and the matrices A (2 k ) i and B (2 k ) ; see Tables 4.2 and 4.3.Finally, we combine the boundary estimate (4.59) with the interior estimate (4.84)to bound the truncation error via (4.48). In terms of r k , (4.59) gives(4.85) (cid:13)(cid:13)(cid:13) γ k h / (cid:13)(cid:13)(cid:13) / ,ε ≤ | V | p I "
12 + εr k r I I r k r − k − k , k ≥ . JON WILKENING
Table 4.2 K (2 k ) i before and after loop (4.83) . k K (2 k )0 before K (2 k )1 before K (2 k )2 K (2 k )0 after K (2 k )1 after0 9 . × − . × +00 . × +00 . × +00 . × +00 . × +00 . × +00 . × +00 . × +00 . × +01 . × +00 . × +01 . × +01 . × +01 . × +01 . × +01 . × +01 . × +02 . × +02 . × +02 . × +02 . × +02 . × +03 . × +03 . × +03 . × +03 . × +04 . × +04 . × +04 . × +04 . × +05 . × +05 . × +05 . × +05 . × +05 . × +06 . × +07 . × +07 . × +07 . × +07 . × +08 . × +08 . × +08 . × +08 . × +08 . × +10 . × +10 . × +10 . × +10 . × +10
10 1 . × +12 . × +12 . × +12 . × +12 . × +12 Table 4.3 e K (2 k ) i before and after loop (4.83) . k e K (2 k )0 before e K (2 k )1 before e K (2 k )2 e K (2 k )0 after e K (2 k )1 after0 1 . × +02 . × +02 . × +02 . × +02 . × +02 . × +02 . × +03 . × +03 . × +03 . × +03 . × +03 . × +03 . × +04 . × +04 . × +04 . × +04 . × +04 . × +04 . × +04 . × +05 . × +05 . × +05 . × +05 . × +05 . × +05 . × +06 . × +06 . × +06 . × +06 . × +07 . × +08 . × +08 . × +08 . × +08 . × +08 . × +09 . × +10 . × +10 . × +10 . × +10 . × +11 . × +11 . × +11 . × +11 . × +12 . × +13 . × +13 . × +13 . × +13 . × +13
10 2 . × +15 . × +15 . × +15 . × +15 . × +15 We now define(4.86) ρ k = h max (cid:16) K (2 k )0 , K (2 k )1 + (cid:17)i − k +2 , k = 0 , h max (cid:16) K (2 k )0 + e K (2 k − , K (2 k )1 + e K (2 k − + , ρ − kk − (cid:17)i − k +2 , k ≥ , so that (4.49), (4.84), and (4.85) imply(4.87) (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + ε k +2 (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) ≤ p I ( | V | + | V | ) " ρ − k − k + 15 εr k r I I r k εr k (cid:19) k +2 . To simplify this expression, we define(4.88) θ k = 15 ρ k +2 k r k and summarize the main result of this section as a theorem. RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION Table 4.4 ρ k and θ k . k ρ k θ k . × +00 . × − . × − . × − . × − . × − . × − . × − . × − . × −
10 0.232 1 . × −
11 0.218 1 . × −
12 0.204 1 . × −
13 0.193 1 . × −
14 0.183 1 . × −
15 0.173 5 . × −
16 0.164 3 . × −
17 0.157 1 . × −
18 0.149 6 . × −
19 0.143 2 . × −
20 0.137 1 . × −
21 0.131 2 . × −
22 0.126 8 . × −
23 0.122 2 . × −
24 0.117 5 . × −
25 0.113 1 . × − Theorem 4.8.
Suppose k ≥ , h ∈ C k +1 , ( T ) , < h ( x ) ≤ for ≤ x ≤ , and ε ≤ r / . Then the truncation errors ψ (2 k ) err and Q (2 k ) err in (4.37) satisfy the bound (4.89) √ (cid:12)(cid:12)(cid:12) Q (2 k ) err (cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε ≤ (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + ε k +2 (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) ≤ p I ( | V | + | V | ) " θ k εr k r I I ερ k r k (cid:19) k +2 , where I m , r k , ρ k , and θ k were defined in (2.9), (4.61), (4.86), and (4.88). Remark
The constants in this estimate have been organized to be either(1) given in the problem statement or easily computable from h ; or (2) difficult tocompute but universal (independent of h ). The first 26 constants in the latter category( ρ k and θ k ) are given in Table 4.4. We have therefore identified the features of h thatare most likely to affect the validity of the lubrication approximation. In particular,higher derivatives are allowed to be large in regions where h is small (since r k dependson the uniform norms of the products ℓ ! h ℓ − ∂ ℓx h rather than on ∂ ℓx h alone). Remark
The term ρ − kk − in the definition of ρ k ensures that ρ k +2 k is anonincreasing function of k . This assumption is useful in the following section whenderiving a bound on the truncation error of the pressure. Note that ρ k itself is allowedto increase as long as ρ k +2 k does not. It is probably not necessary to include this termin the definition of ρ k since it is not the argmax for 1 ≤ k ≤
25, and by that point θ k (and hence ρ k +2 k ) appears to be decreasing rapidly without it; see Figure 4.1. We now show how to use the errorbound we have obtained for the stream function to bound the error in the velocity,0
JON WILKENING ρ k−1 and log θ k versus k ρ k − l og θ k k ρ k−1 log θ k Fig. 4.1 . Plot of ρ − k and log θ k versus k . Note that ρ − k initially decreases but eventuallygrows almost linearly, indicating that the lubrication expansion is probably an asymptotic seriesrather than a convergent series. The term involving θ k in (4.89) is only important when k is small,since θ k converges rapidly to zero as k → ∞ . vorticity, and pressure. We define u (2 k ) , v (2 k ) , ω (2 k ) , and p (2 k ) in terms of ψ (2 k ) as in(3.8) and define, e.g., ω (2 k ) err = ω exact − ω (2 k ) approx , ω (2 k ) approx = ω (0) + ε ω (2) + · · · + ε k ω (2 k ) . (4.90)From (3.8), we then have (cid:16) u (2 k ) err , v (2 k ) err (cid:17) = ∇ × ψ (2 k ) err , ω (2 k ) err = − ∆ ε ψ (2 k ) err − ε k +2 ψ (2 k ) xx , (4.91) ∂ x p (2 k ) err = − ∂ y ω (2 k ) err , ∂ y p (2 k ) err = ∂ x (cid:0) ε ω exact (cid:1) , k = 0 ,∂ x (cid:16) ε ω (2 k − err (cid:17) , k ≥ , (4.92)which immediately gives bounds on the error in velocity and vorticity:(4.93) (cid:18)(cid:13)(cid:13)(cid:13) u (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + (cid:13)(cid:13)(cid:13) εv (2 k ) err (cid:13)(cid:13)(cid:13) ,ε (cid:19) / ≤ (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε ≤ ( ∗ ) , k ≥ , (cid:13)(cid:13)(cid:13) ω (2 k ) err (cid:13)(cid:13)(cid:13) ≤ (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + ε k +2 (cid:13)(cid:13)(cid:13) ψ (2 k ) xx (cid:13)(cid:13)(cid:13) ≤ ( ∗ ) , k ≥ , where ( ∗ ) is the right-hand side of (4.89). Obtaining a bound on the error in thepressure is somewhat more difficult, as it relies on the fact that the gradient is anisomorphism from L (Ω) (the space of square integrable functions with zero mean)onto the polar set(4.94) V = (cid:8) f ∈ H − (Ω) : h f , u i = 0 whenever u ∈ V (cid:9) , where V = { u ∈ H (Ω) : ∇ · u = 0 } . Given f ∈ V , there is a unique p ∈ L (Ω)such that ∇ p = f ; moreover, p satisfies(4.95) k p k ≤ β − | f | − , β = inf p sup u | ( p, ∇ · u ) |k p k | u | . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION f . More precisely, as k · k and | · | are equivalent on H (Ω), the negative norms(4.96) | f | − = sup | u | + | v | =1 |h f , ( u ; v ) i| , k f k − = sup k u k + k v k =1 |h f , ( u ; v ) i| are equivalent on H − (Ω) . Since | u | ≤ k u k for all u ∈ H (Ω), we have k f k − ≤| f | − for all f ∈ H − (Ω) .Explicit estimates [7, 23, 10] for the LBB constant β in (4.95) have been obtainedfor rectangular domains (with no periodicity), e.g.,(4.97) 1 ℓ sin π ≤ β ≤ π √ ℓ , ℓ = max (cid:18) L L , L L (cid:19) , Ω = (0 , L ) × (0 , L ) . The lower bound here also works for an x -periodic rectangle as the condition that u | x =0 = 0 and u | x = L = 0 is more restrictive than u | x =0 = u | x = L . Explicit estimatesare also known for domains that are star-shaped with respect to each point in a ballof radius R contained inside Ω; see [13, 23]. Our interest in the present work is in x -periodic domains with the upper boundary given by a function h ( x ). Such domainsare not in general star-shaped, so the previously known results do not apply. In[24], we improve the estimate for the lower bound on β in (4.97) for an x -periodicrectangle by a factor of about 3.5 and show how to avoid invoking Rellich’s theoremin the change of variables to the case that Ω is x -periodic with the upper boundarygiven by h ( x ). It is shown that β − ≤ (cid:0) r − (cid:1) (cid:18) h h (cid:19) / max (cid:18) , Lh , h h (cid:19) , h = min ≤ x ≤ L h ( x ) ,h = max ≤ x ≤ L h ( x ) , where r = max( k h x k ∞ , k hh xx k / ∞ ) − and L is the period of h ( x ). In the currentcase, the length scales ¯ H and ¯ W were chosen so that L = 1 and h ≤ ∇ p = f yields(4.98) k p k ≤ β − | f | − , β − ≤ max (cid:18) h − / , h − / (cid:19) (cid:0) r − (cid:1) . The dependence on gap thickness h occurs because p can change rapidly in the gapwithout a large penalty from u . In [24], an example is given to show that the factorof h − / in the formula (4.98) for β − cannot be improved.We have reduced the problem of bounding p (2 k ) err to that of bounding the functionalon the right-hand side of ∇ p (2 k ) err = f k in (4.92), namely,(4.99) h f k , ( u ; v ) i = (R Ω ω (0) err u y − ε ω exact v x dA, k = 0 , R Ω ω (2 k ) err u y − ε ω (2 k − err v x dA, k ≥ . First, we check that f k belongs to V . If u, v ∈ H (Ω) and u x + v y = 0, the function φ ( x, y ) = R y u ( x, η ) dη satisfies φ y = u , φ x = − v and so belongs to Ψ. As a result,(4.100) h f , ( u ; v ) i = Z Ω − ω (0) φ yy + ω exact (cid:0) φ yy + ε φ xx (cid:1) dA = Z Ω ψ (0) yy φ yy + (∆ ε ψ exact )(∆ ε φ ) dA = 0 , JON WILKENING and for k ≥ h f k , ( u ; v ) i = z }| { h f k − , ( u ; v ) i + Z Ω − ε k ω (2 k ) φ yy − ε k ω (2 k − φ xx dA = ε k Z Ω ψ (2 k ) yy φ yy + ψ (2 k − xx φ yy + ψ (2 k − yy φ xx + ψ (2 k − xx φ xx | {z } omit if k = 1 dA = 0 . Next, we bound the norm of f k . From (4.99), we see that | f | − ≤ (cid:13)(cid:13)(cid:13) ω (0) err (cid:13)(cid:13)(cid:13) + (cid:13)(cid:13) ε ω exact (cid:13)(cid:13) , | f k | − ≤ (cid:13)(cid:13)(cid:13) ω (2 k ) err (cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13) ε ω (2 k − err (cid:13)(cid:13)(cid:13) , k ≥ . Denoting the right-hand side of (4.89) by ( ∗ ), we claim that(4.101) ε r k ω exact k ≤ ( ∗ ) , ( k = 0) , ε r k (cid:13)(cid:13)(cid:13) ω (2 k − err (cid:13)(cid:13)(cid:13) ≤ ( ∗ ) , ( k ≥ . Once this is shown to be true, we will have the following bound on | f k | − :(4.102) | f k | − ≤ (cid:0) r k (cid:1) ( ∗ ) , ( k ≥ . Together with (4.98) and the fact that r k ≤ r for k ≥
0, this will give(4.103) (cid:13)(cid:13)(cid:13) p (2 k ) err (cid:13)(cid:13)(cid:13) ≤ max (cid:18) h − / , h − / (cid:19) (cid:0) r k + r − k (cid:1) ( ∗ ) . Note that there is at least a power of r − k in ( ∗ ) to prevent this bound from divergingas r k → ∞ . It may be possible to improve the bound in this regime by replacing k ε ω exact k in (4.101) by k ε ∂ x ω exact k − , but this seems very difficult. At any rate,if r k = ∞ , then ( ∗ ) = 0, h ( x ) is a constant function, the exact and approximatevorticity are constants, f k is the zero functional, and p exact = p (2 k ) err = 0. Let us nowprove (4.101). For k ≥
1, this follows from(4.104) " ρ − kk − + 15 εr k − r I I r k − εr k − (cid:19) k (cid:18) εr k (cid:19) ≤ " ρ − k − k + 15 εr k r I I r k εr k (cid:19) k +2 , which holds because ρ k +2 k and r k are nonincreasing functions of k ; see Remark 4.10and the definition of r k in (4.61). For k = 0, we use Theorem 4.4 to conclude that(4.105) k ω exact k ≤ k ψ k ,ε ≤ (cid:18)(cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε + (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε (cid:19) , where g ( x ) = V and g ( x ) = (1 + ε h x ( x ) ) − / . Now(4.106) (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε ≤ V (cid:13)(cid:13)(cid:13) h − / (cid:13)(cid:13)(cid:13) ,ε = V Z h − + 14 ( εh x ) h − dx ≤ V " I + 14 (cid:18) εr (cid:19) I RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION (cid:13)(cid:13)(cid:13) h − / g (cid:13)(cid:13)(cid:13) / ,ε ≤ V (cid:13)(cid:13)(cid:13) h − / (cid:0) ε h x (cid:1) − / (cid:13)(cid:13)(cid:13) ,ε (4.107)= V Z (cid:0) h − (cid:1) ( · ) − + (cid:20) − h − / ( εh x )( · ) − / − (cid:16) h − / (cid:17) ( · ) − / ε h x h xx (cid:21) dx ≤ V Z h − + 12 h − ( εh x ) + 8 h − ( εh x ) (cid:18) ε hh xx (cid:19) dx ≤ V " I + (cid:18)
12 + 8 ε r (cid:19) (cid:18) εr (cid:19) I . (4.108)Since we have assumed that ε ≤ r /
3, we conclude that(4.109) ε r k ω exact k ≤ p I ( | V | + | V | ) "
15 + r
12 + 881 ! εr r I I ε r . Comparing this to (4.87) with k = 0 and noting from Table 4.4 that ρ − ≥
15, weobtain (4.101) as claimed. Thus, we have proved the following theorem.
Theorem 4.11.
Suppose k ≥ , h ∈ C k +1 , ( T ) , < h ≤ h ( x ) ≤ for x ∈ T ,and ε ≤ r / . Then the truncation errors of the stream function, flux, velocity,vorticity, and pressure satisfy the bounds (cid:13)(cid:13)(cid:13) ψ (2 k ) err (cid:13)(cid:13)(cid:13) ,ε ≤ ( ∗ ) , | Q (2 k ) err | ≤ ( ∗ ) √ , (cid:18)(cid:13)(cid:13)(cid:13) u (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + (cid:13)(cid:13)(cid:13) εv (2 k ) err (cid:13)(cid:13)(cid:13) ,ε (cid:19) / ≤ ( ∗ ) , (cid:13)(cid:13)(cid:13) ω (2 k ) err (cid:13)(cid:13)(cid:13) ≤ ( ∗ ) , (cid:13)(cid:13)(cid:13) p (2 k ) err (cid:13)(cid:13)(cid:13) ≤ max (cid:18) h − / , h − / (cid:19) (cid:0) r k + r − k (cid:1) ( ∗ ) , (4.110) where T = [0 , p is the periodic unit interval ( ∗ ) = p I ( | V | + | V | ) " θ k εr k r I I ερ k r k (cid:19) k +2 , (4.111) r k = max ≤ ℓ ≤ k +2 ((cid:13)(cid:13)(cid:13)(cid:13) ℓ ! h ℓ − ∂ ℓx h (cid:13)(cid:13)(cid:13)(cid:13) /ℓ ∞ )! − , I m = Z h ( x ) − m dx, (4.112) and ρ k , θ k are constants independent of h that can be computed once and for all asdescribed in section and listed in Table
5. Finite element validation.
In this section, to test the error bounds of The-orem 4.11, we compute u (2 k ) err , v (2 k ) err , p (2 k ) err , ω (2 k ) err numerically for the simple geometrydescribed by(5.1) h ( x ) = 1 + a − a πx ) , Case 1: a = 1 / , Case 2: a = 1 / , with boundary conditions V = − . V = 1. We do this by comparing u (2 k ) approx , v (2 k ) approx , p (2 k ) approx , ω (2 k ) approx in (4.90) to finite element solutions of the Stokes equationson appropriately rescaled geometries.4 JON WILKENING −9 −7 −5 −3 −1 ε no r m o f e rr o r Actual errors k = , s l ope s = . , . , . k = , s l o p e s = . , . , . k = , s l o p e s = . , . , . k = , s l o p e s = . , . , . u p ω FE error up ω −7 −5 −3 −1 ε e rr o r Error bounds from Theorem 4.11 k = k = k = k = k = Fig. 5.1 . Top: plot of k u (2 k ) err k ,ε = ( k u (2 k ) err k ,ε + k εv (2 k ) err k ,ε ) / , k p (2 k ) err k , and k ω (2 k ) err k for a = 1 / , . ≤ ε ≤ . , and k ∈ { , , , , } . The slopes of the lines were computed via linearregression using the smallest values of ε for which the finite element solution is trusted ( ε ≥ . for k = 10 and ε ≥ . for k = 20 ). As expected, for fixed k , the error is O ( ε k +2 ) . Bottom: plotof the error bound ( ∗ ) in (4.111) , using V = − . , V = 1 , I = 2 . , I = 24 . , and r k = 0 . for k ≥ , as appropriate for h ( x ) in (5.1) with a = 1 / . The results are summarized in Figures 5.1–5.5. For 21 values of ε spaced expo-nentially between ε = 0 .
04 and ε = 0 .
3, we set up a logically rectangular, M × N finite element mesh on the domain(5.2) Ω ε = { ( x, y ) : 0 ≤ x ≤ , < y < εh ( x ) } . RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION x y ω -1.2-2-2.8-3.6-4.4-5.2-6-6.8-7.6-8.4-9.2-10 x y ω0 x y ω4 x y ω10 x y ω20 Fig. 5.2 . Contour plots of ω exact , ω (0) err , ω (4) err , ω (10) err , and ω (20) err for h ( x ) in (5.1) with a = 1 / , V = − . , V = 1 . , and ε = 0 . . Each of these plots corresponds to one of the markers inFigure . JON WILKENING x y p x y p0 x y p4 x y p10 x y p20 Fig. 5.3 . Contour plots of p exact , p (0) err , p (4) err , p (10) err , and p (20) err for h ( x ) in (5.1) with a = 1 / and ε = 0 . . The “exact” solution was computed using a least squares finite element method with node quartic triangular elements on a × grid. RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION −6 −4 −2 ε no r m o f e rr o r Error bounds and actual errors u p ω k = , s l ope s = . , . , . k = , s l o p e s = . , . , . k = , s l o p e s = . , . , . k = , s l o p e s = . , . , . a−priori estimatesactual errors k = k = k = k = up ω Fig. 5.4 . Error estimates and actual errors with a = 1 / . x y ω x y p x y ω10 x y p10 Fig. 5.5 . Plots of ω exact , p exact , ω (10) err , and p (10) err with a = 1 / , N = 64 , M = 4800 . JON WILKENING
The mesh points are aligned vertically with equal spacing ∆ y = h ( x ) /N , while thegrid spacing in the x -direction is chosen to keep the aspect ratios of the grid cells asclose to 1 as possible; we do this by solving an ODE to enforce ∆ x ≈ h ( x ) /N , whichalso determines M . For a = 1 /
5, we use N = 96 with M ranging from 768 to 5376as ε ranges from 0 . .
04; for a = 1 / N = 64 with M ranging from1600 to 10368. Four-by-four blocks of neighboring grid cells are merged and cut intotwo 15 node triangles. Interior nodes of the triangles are adjusted to keep the edgesstraight except on the top boundary, where we use quartic isoparametric elements.We solve the Stokes equations on this mesh using a least squares finite element methodsimilar to [6] but using quartic elements to model the velocity components u and v ,the pressure p , the vorticity ω = v x − u y , and two strain rates τ = u y + v x and γ = v y − u x . We use multigrid to solve the resulting system of equations, which takesfrom 3 to 15 minutes on a 2.4 GHz desktop machine with 16 GB RAM.Once the finite element solution is known at the grid points, we normalize thevelocity, pressure, and vorticity as described in section 2 and rescale the domainfrom Ω ε to Ω. We then use the method described in Appendix A to compute ψ (0) , ψ (2) , . . . , ψ (20) and their derivatives through order 3 at the grid points. Next, we usethe formulas in (3.8) to obtain u (2 k ) , v (2 k ) , ω (2 k ) , and p (2 k ) for k = 0 , . . . ,
10. Forpressure, we use 20 point Gaussian quadrature to integrate p (2 k ) x along the x -axis todetermine p (2 k ) ( x,
0) at the mesh nodes. The integration of p (2 k ) y in the y -direction isdone analytically. With the expansion coefficients in hand, we evaluate u (2 k ) err = u exact − u (2 k ) approx , u (2 k ) approx = u (0) + ε u (2) + · · · + ε k u (2 k ) , (5.3)etc., at the grid nodes, where we use the finite element solution for u exact . We thenrun through the triangles and sum up the local contributions to the errors(5.4) (cid:13)(cid:13)(cid:13) u (2 k ) err (cid:13)(cid:13)(cid:13) ,ε + (cid:13)(cid:13)(cid:13) εv (2 k ) err (cid:13)(cid:13)(cid:13) ,ε , (cid:13)(cid:13)(cid:13) p (2 k ) err (cid:13)(cid:13)(cid:13) , (cid:13)(cid:13)(cid:13) ω (2 k ) err (cid:13)(cid:13)(cid:13) by interpolating the values at the grid nodes and integrating the resulting polynomialson the triangle; this step is very similar to the assembly of the stiffness matrix. Finally,we store the results in a file for visualization (see Figures 5.2, 5.3, and 5.5) andrecord the norms of the truncation errors for comparison with the error bounds ofTheorem 4.11.The results of this comparison are shown in Figures 5.1 and 5.4. As expected,for fixed k , the actual errors decay as O ( ε k +2 ). The a priori error bounds eventuallydecrease like O ( ε k +2 ) as well, but the term involving θ k in (4.111) is significant overthis range of ε in some of the cases, causing the slopes to be larger: θ k r k r I I = k = 0 k = 1 k = 2 k = 5 k = 10 a = 1 / . .
94 0 . . . × − a = 1 /
100 257 19 . . .
020 2 . × − . This effect is much more pronounced when a = 1 /
100 in (5.1) due to(5.5) r I I = 12 r
32 + 1 a + 32 a = ( . , a = 1 / , . , a = 1 / . The deviation from linearity in the plots of “actual error” for small ε and large k isdue to error in the finite element solutions, which are accurate to about 9 digits. This RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION a = 1 /
100 since the pressure and vorticity of the exact solution inthe vicinity of the narrow gap increases as a decreases, and also because we were forcedto use a coarser mesh with a = 1 /
100 to avoid running out of computer memory in thefinite element simulations. The data points with ε = 0 .
099 in Figure 5.1 correspondto the contour plots in Figures 5.2 and 5.3, where we plot ω exact , ω (2 k ) err , p exact , and p (2 k ) err for 2 k = 0 , , ,
20. The data points with ε = 0 .
099 in Figure 5.4 correspond tothe contour plots in Figure 5.5. We remark that the apparently large value of p (10) err inthe narrow gap in Figure 5.5 is due to smoothing in the least squares finite elementsolver; the expansion solution is more accurate than the finite element solution in thisregion of the domain. The error patterns that emerge in all these cases are ratherinteresting, indicating that the spaces H k in Theorem 3.3 (the structure theorem)can be quite complicated even for simple curves h ( x ).Although our estimates for the error in pressure include an additional factor of h − / ( r k + r − k ) , all our numerical experiments (including complicated geometries inwhich the inf-sup constant β − does exhibit h − / behavior) indicate that k p (2 k ) err k iscomparable to k ω (2 k ) err k . In fact, for large k , pressure seems to be the most accuratelycomputed variable; see Figures 5.1 and 5.4. We do not know how to explain this asthe pressure is determined by solving (4.92), which involves inverting the operator ∇ : L (Ω) → H − (Ω) . For some reason, in lubrication-type problems, the right-hand side f k belongs to a subspace of H − (Ω) that is not amplified by ( ∇ ) − whensolving ∇ p (2 k ) err = f k .The following table shows the minimum ratio of the a priori error estimate to theactual error k u (2 k ) err k ,ε for the data points in Figures 5.1 and 5.4 that were used tocompute the slopes of the best-fit lines: k = 0 1 2 5 10(min ratio, a = 1 / / (2 k +2) . . . . . a = 1 / / (2 k +2) . . . . − . For example, in the 10 calculations (with ε ranging from 0 . ≤ ε ≤ . k = 4 line in Figure 5.1, the ratios of the a priorierrors to the exact errors ranged between 1 . × and 1 . × , so we recorded √ . × ≈ .
0. This table gives information on how far the values ρ k in Table 4.4are from their optimal values. For example, if we increased ρ by more than a factorof 2 . θ fixed, the estimate (4.110) would fail to hold for this geometry.Since r − k in (4.61) is used as a convenient upper bound on all the integrals | E (2 ℓ ) m,j | / ℓ and | e E (2 ℓ ) m,j | / ℓ that arise in the definition of Q (2 k ) and also in the bounds for k ψ (2 k ) xx k and k h ψ (2 k − xxxx k , it is remarkable that the values of ρ k we computed are within afactor of 3 of optimal for k = 5, k = 10, and perhaps all k ≥
6. Discussion.
Although we are able to estimate the effective radius of conver-gence ρ k r k quite closely, our estimates of k ψ (2 k ) err k ,ε , k ω (2 k ) err k , etc., are likely to beseveral orders of magnitude too large. One shouldn’t expect an a priori bound thatholds for all geometries alike to provide an exceptionally sharp bound for any spe-cific geometry. Instead, our analysis provides a clear picture of the features of h ( x )that cause the effective radii of curvature r k ρ k to become small, namely, large valuesof h k − ∂ kx h . No previous study has ever described how the constant hidden in the O ( ε k +2 ) depends on h ; instead, h has always been fixed at the outset and only thelimit as ε → JON WILKENING
Approximation of velocity at a point ε u ( / , / ) FE solution
Fig. 6.1 . Comparison of u (2 k ) approx (solid lines) to u exact (dots) at the point ( x, y ) = ( , h ) for k = 0 , , , , , , , , . Here h ( x ) = + sin(2 πx ) , V = − . , and V = 1 . This function h ( x ) is real analytic and periodic, yet the expansion solution appears to be an asymptotic seriesrather than a convergent series. Another feature of this analysis is that it separates the constants into two types:those that are (1) given in the problem statement or easily computable from h ; or(2) difficult to compute but universal (independent of h ). We listed the first severalconstants in the latter category ( ρ k and θ k ) in Table 4.4. It is interesting that ρ k actually increases until 2 k = 10 and doesn’t get as bad as ρ again until 2 k = 26.However, at that point it seems to be decreasing steadily like 1 /k , indicating thatthe effective radius of curvature in our a priori error bound will shrink to zero as k → ∞ . The reason for this is that the recurrences (A.3) and (A.4) relating thematrices A (2 k ) i and B (2 k ) to their lower order counterparts cause the norms of thesematrices to grow like k !. Thus, although ρ k involves k th roots of these constants,these k th roots still grow linearly in k . On the other hand, if h is real analytic as wellas periodic, a standard contour integral argument shows that there is an r > k ∂ kx h k ∞ ≤ k ! r − k for all k ≥
0; thus, the constants r k will remain bounded awayfrom zero. For example, if h ( x ) is of the form (5.1), one may show that if a ∈ (0 , / k ℓ ! h ℓ − ∂ ℓx h k /ℓ ∞ occurs when ℓ = 2, so all the r k are equalto r = ( π √ − a ) − . It is conceivable that when h is real analytic, the norms of thefunctions ψ (2 k ) grow slowly enough that the stream function expansion converges inspite of the fact that the matrices A (2 k ) i and B (2 k ) in their representation (3.22) blowup like k !. This would simply mean that we chose a bad basis in terms of which torepresent ψ . We used orthogonal polynomials in (4.77) to improve this basis, but theremay be other improvements. Figure 6.1 shows that this is not the case. Even when h varies sinusoidally, the expansion solution appears to be an asymptotic series ratherthan a convergent series: all the variables, including the flux terms Q (2 k ) , appear togrow like k ! as k becomes large.Nevertheless, the expansion solutions can be extremely accurate (almost exact) aslong as they are used for a geometry that falls within the effective radius of convergence RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION
Appendix A. Implementation.
We have developed two methods for comput-ing the higher order corrections described in section 3.2 using a computer. In the first,we use Mathematica to evaluate the derivatives and antiderivatives in recursion (3.6)and Algorithm 3.1 symbolically. With this approach, the main challenge occurs atthe step where Q (2 k ) is defined as a definite integral. We do this through patternmatching and symbol replacement. At the stage where the definite integral is to beevaluated, we replace all instances of ∂ jx h in the integrand by j ! t j /h j − . Each term inthe result (call it R ) will contain a factor of h − or h − , with no other dependence on h . For each k = k , . . . , j = 1 , . . . , d k , we find the terms in R that contain ϕ (2 k ) j (left in the form t i · · · t i k k described in Algorithm 3.2) as a factor. These terms areremoved from R while their symbolic integrals (with ϕ (2 k ) j /h m replaced by I m E (2 k ) m,j )are divided by 2 I and added to the desired flux Q (2 k ) . By running through the ϕ (2 k ) j in decreasing k order, we convert higher order products (e.g., t t /h ) into symbols(e.g., I E (4)3 , ) before one of their lower order factors can be converted incorrectly (e.g.,into I E (2)3 , t ). This approach is effective through 6th or 8th order but becomes ratherslow as the complexity of the expansion increases.The second approach is much faster and can be implemented in any modern pro-gramming language. We have written a version in C ++ and a version in Mathematica.Instead of representing the basis functions ϕ ( k ) j for H k using a computer algebra sys-tem, we represent them as ( k + 1)-tuples of integers. For example, the functions 1, h x , h h x h xxx , and h h x h xx h xxxx in H , H , H , and H are represented by (0),(0 , , , , , , , , , , , , , , i , . . . , i k ) represents a basisfunction for H k iff(A.1) i + 2 i + · · · + ki k = k, i = i + 2 i + · · · + ( k − i k . We begin by constructing the basis sets Φ k for 0 ≤ k ≤ k and storing them as( k + 1) × d k integer matrices with columns corresponding to the ϕ ( k ) j . This is doneusing Algorithm 3.2, which returns the columns sorted lexicographically from the lastslot to the first slot (e.g., (3 , , , T < (2 , , , T < (3 , , , T ). Sorted columnsallow us to find the column index corresponding to a given tuple in log d k time.Next, for 0 ≤ k ≤ k −
1, we compute the operators h∂ x and h x · from H k to H k +1 and store them as sparse integer matrices of dimension d k +1 × d k . If column J of Φ k contains the tuple ( i , . . . , i k ), we define i k +1 = 0 and compute h x · : ( i , . . . , i k ) ( i , i + 1 , i , . . . , i k +1 ) , (A.2) h∂ x : ( i , . . . , i k ) X { r : i r =0 } i r ( r + 1)( i + 1 , . . . , i r − , i r +1 + 1 , . . . , i k +1 ) , where the omitted indices are unmodified and the +1 and − r = 0 in the sum. The factor of ( r + 1) is due to the factorials in the definition ofthe ϕ ( k ) j . The column index l of each ( k + 2)-tuple in the result is found in Φ k +1 , andthe corresponding coefficient (1 or i r ( r + 1)) is added to the l th row and J th columnof the sparse matrix representing h∂ x or h x · . The entries of these sparse matricesare positive, and the column sums (i.e., 1-norms) are all equal to 1 for h x · and to i + 2 i + · · · + ( k + 1) i k = 2 k for h∂ x (by (A.1)).2 JON WILKENING
Once the operators h∂ x and h x · are known, we use them to recursively computethe matrices A (2 k ) = V A (2 k )0 + V A (2 k )1 and B (2 k ) in (3.22). We start by setting A (0)0 = (0 , , − , T , A (0)1 = (0 , , − , T , and B (0) = (0 , , , − T as in Example 3.5.For 1 ≤ k ≤ k , we mimic the proof of Theorem 3.3 to build up A (2 k ) and B (2 k ) rowby row. For 4 ≤ n ≤ k + 3 and i = 0 ,
1, we use sparse matrix–vector multiplicationto define the rows(A.3) A (2 k ) i ( n, :) = − h∂ x − ( n − h x ][ h∂ x − ( n − h x ] h A (2 k − i ( n, :) T i n ( n − T ,B (2 k ) ( n, :) = − h∂ x − ( n − h x ][ h∂ x − ( n − h x ] (cid:2) B (2 k − ( n, :) T (cid:3) n ( n − ! T . If k ≥
2, then for 6 ≤ n ≤ k + 3, we add the following vectors to A (2 k ) i ( n, :) and B (2 k ) ( n, :), respectively:(A.4) (cid:18) − [ h∂ x − ( n − h x ][ h∂ x − ( n − h x ][ h∂ x − ( n − h x ][ h∂ x − ( n − h x ] h A (2 k − i ( n, :) T i n ( n − n − n − (cid:19) T , (cid:18) − [ h∂ x − ( n − h x ][ h∂ x − ( n − h x ][ h∂ x − ( n − h x ][ h∂ x − ( n − h x ] [ B (2 k − ( n, :) T ] n ( n − n − n − (cid:19) T . Next we zero out rows 0 and 1 of A (2 k )0 , A (2 k )1 , B (2 k ) , and set(A.5) A (2 k ) i (2 , :) = k +3 X n =4 ( n − A (2 k ) i ( n, :) , A (2 k ) i (3 , :) = k +3 X n =4 (2 − n ) A (2 k ) i ( n, :) ,B (2 k ) (2 , :) = k +3 X n =4 ( n − B (2 k ) ( n, :) , B (2 k ) (3 , :) = k +3 X n =4 (2 − n ) B (2 k ) ( n, :) . Finally, we subtract (cid:0) − / k (cid:1) from A (2 k )1 (2 ,
1) and add it to A (2 k )1 (3 ,
1) to account forthe boundary data, where we recall that the rows and columns are indexed starting at0 and 1, respectively. Using this approach, our C ++ code can compute these matricesthrough order 2 k = 50 using floating point arithmetic in a few seconds, while ourMathematica code can compute through order 2 k = 30 in exact rational arithmeticin about an hour. This allows us to explore the properties of the stream functionexpansion and test our error estimates to quite a high order. REFERENCES[1]
M. Abramowitz and I. A. Stegun , Handbook of Mathematical Functions with Formulas,Graphs, and Mathematical Tables , Dover, New York, 1964.[2]
A. Assemien, G. Bayada, and M. Chambat , Inertial effects in the asymptotic behavior of athin film flow , Asympt. Anal., 9 (1994), pp. 177–208.[3]
G. Bayada and M. Chambat , The transition between the Stokes equations and the Reynoldsequation: A mathematical proof , Appl. Math. Optim., 14 (1986), pp. 73–93.[4]
G. Bayada and M. Chambat , Mod´elisation de la jonction d’un ´ecoulement tridimensionnelet d’un film mince bidimensionnel , C. R. Acad. Sci. Paris, Ser. I, 309 (1989), pp. 81–84.[5]
D. Braess , Finite Elements—Theory, Fast Solvers, and Applications in Solid Mechanics ,Cambridge University Press, Cambridge, UK, 1997.RROR ESTIMATES FOR REYNOLDS’ APPROXIMATION [6] Z. Cai, T. A. Manteuffel, and S. F. McCormick , First-order system least squares for theStokes equations, with application to linear elasticity , SIAM J. Numer. Anal., 34 (1997),pp. 1727–1741.[7]
E. V. Chizhonkov and M. A. Olshanskii , On the domain geometry dependence of the LBBcondition , M2AN Math. Model. Numer. Anal., 34 (2000), pp. 935–951.[8]
G. Cimatti , How the Reynolds equation is related to the Stokes equations , Appl. Math. Optim.,10 (1983), pp. 267–274.[9]
I. Ciuperca, I. Hafidi, and M. Jai , Singular perturbation problem for the incompressibleReynolds equation , Electron. J. Differential Equations, 2006 (2006), pp. 1–19.[10]
M. Dobrowolski , On the LBB constant on stretched domains , Math. Nachr., 254–255 (2003),pp. 64–67.[11]
A. Duvnjak and E. Maru˘si´c-Paloka , Derivation of the Reynolds equation for lubrication ofa rotating shaft , Arch. Math., 36 (2000), pp. 239–253.[12]
H. G. Elrod , A derivation of the basic equations for hydrodynamic lubrication with a fluidhaving constant properties , Quart. Appl. Math., XVII (1960), pp. 349–359.[13]
G. P. Galdi , An Introduction to the Mathematical Theory of the Navier-Stokes Equations,Vol. Linearized Steady Problems , Springer–Verlag, New York, 1994.[14]
V. Girault and P.-A. Raviart , Finite Element Methods for Navier–Stokes Equations ,Springer–Verlag, Berlin, 1986.[15]
J. Kevorkian and J. D. Cole , Multiple Scale and Singular Perturbation Methods , Springer-Verlag, New York, 1996.[16]
W. E. Langlois , Slow Viscous Flow , Macmillan, New York, 1964.[17]
I. Moise, R. Temam, and M. Ziane , Asymptotic analysis of the Navier-Stokes equations inthin domains , Topol. Methods Nonlinear Anal., 10 (1997), pp. 249–282.[18]
S. A. Nazarov , Asymptotic solution of the Navier-Stokes problem on the flow of a thin layerof fluid , Siberian Math. J., 31 (1990), pp. 296–307.[19]
S. A. Nazarov , Asymptotics of the Stokes system solutions at a surfaces contact point , C. R.Acad. Sci. Paris, S´er. I, 312 (1991), pp. 207–211.[20]
C. Pozrikidis , Introduction to Theoretical and Computational Fluid Dynamics , Oxford Uni-versity Press, New York, 1997.[21]
G. Raugel and G. R. Sell , Navier-Stokes equations on thin D domains. I: Global attractorsand global regularity of solutions , J. Amer. Math. Soc., 6 (1993), pp. 503–568.[22]
O. Reynolds , On the theory of lubrication and its applications to Mr. Beauchamp Tower’sexperiments, including an experimental determination of the viscosity of olive oil , Philos.Trans. R. Soc. Lond., 177 (1886), pp. 157–234.[23]
G. Stoyan , Iterative Stokes solvers in the harmonic Velte subspace , Computing, 67 (2001),pp. 13–33.[24]
J. Wilkening , Inf-sup estimates for the Stokes problem in a periodic channel , arXiv:0706.4082.[25]