Qualitatively accurate spectral schemes for advection and transport
QQUALITATIVELY ACCURATE SPECTRAL SCHEMES FORADVECTION AND TRANSPORT
HENRY O. JACOBS & RAM VASUDEVAN
Abstract.
The transport and continuum equations exhibit a number of con-servation laws. For example, scalar multiplication is conserved by the trans-port equation, while positivity of probabilities is conserved by the contin-uum equation. Certain discretization techniques, such as particle based meth-ods, conserve these properties, but converge slower than spectral discretizationmethods on smooth data. Standard spectral discretization methods, on theother hand, do not conserve the invariants of the transport equation and thecontinuum equation. This article constructs a novel spectral discretizationtechnique that conserves these important invariants while simultaneously pre-serving spectral convergence rates. The performance of this proposed methodis illustrated on several numerical experiments. Introduction
Let M be a compact n -manifold with local coordinates x , . . . , x n and let X be asmooth vector-field on M , whose local components are given by ( X ( x ) , . . . , X n ( x )).This paper is concerned with the following pair of partial differential equations(PDEs): ∂ t f + X i ∂ i f = 0(1) ∂ t ρ + ∂ i ( ρX i ) = 0(2)for a time dependent function f and a time dependent density ρ on M . In theabove PDEs we are following the Einstein summation convention, and summing overthe index “ i .” Equation (1), which is sometimes called the “transport equation,”describes how a scalar quantity is transported by the flow of X [45]. Equation(2), which is sometimes called the “continuum equation” or “Liouville’s equation”describes how a density (e.g. a probability distribution) is transported by the flowof X . Such PDEs arises in a variety of contexts, ranging from mechanics [4, 45] tocontrol theory [24], and can be seen as zero-noise limits of the forward and backwardKolmogorov equations [35].The solution to (1) takes the form f ( x ; t ) = f ((Φ tX ) − ( x ); 0) ≡ (Φ tX ) ∗ f ( · ; 0)where Φ tX : M → M is the flow map of X at time t [30, Chapter 18]. From thisobserve that (1) exhibits a variety of conservation laws. For example, if f and g are solutions to (1), then so is their product, f · g , and their sum, f + g . Similarly,the solution to (2) takes the form ρ ( x ; t ) = | det( D (Φ tX ) − ( x )) | ρ ((Φ tX ) − ( x ); 0) :=(Φ tX ) ∗ ρ ( · ; 0). One can deduce that the L -norm of ρ ( x ; t ) is conserved in time [30,Theorem 16.42]. Finally, (2) is the adjoint evolution equation to (1) in the sense Date : January 27, 2016. a r X i v : . [ c s . S Y ] J a n HENRY O. JACOBS & RAM VASUDEVAN that the integral (cid:104) f, ρ (cid:105) := (cid:82) f ρ is constant in time . This motivates the followingdefinition of qualitative accuracy: Definition 1.1.
A numerical method for (1) and (2) is qualitatively accurate if itconserves discrete analogs of scalar multiplication/addition, the L -norm and thetotal mass for densities and the sup-norm for functions. Both (1) and (2) can be numerically solved by a variety of schemes. For acontinuous initial condition, f , for example, the method of characteristics [15, 1]describes a solution to (1) as a time-dependent function f ( x t ; t ) = f ( x ) where x t is the solution to ˙ x = X ( x ). This suggests using a particle method to solve for f at a discrete set of points [31]. In fact, a particle method would inherit manydiscrete analogs of the conservation laws of (1), and would as a result be qualitativelyaccurate . For example, given the input h = f · g , the output of a particle methodis identical to the (componentwise) product of the outputs obtained from inputing f and g separately. However, particle methods converge much slower than theirspectral counterparts when the function f is highly differentiable [7, 19].In the case where M is the unit circle, S , a spectral method can be obtained byconverting (2) to the Fourier domain where it takes the form of an Ordinary Dif-ferential Equation (ODE): d ˆ ρ k dt + 2 πik (cid:98) X k − (cid:96) ˆ f (cid:96) where ˆ ρ k and (cid:98) X k denote the Fouriertransforms of ρ and X [44]. In particular, this transformation converts (2) into anODE on the space of Fourier coefficients. A standard spectral Galerkin discretiza-tion is obtained by series truncation.Such a numerical method is good for C k -data, in the sense that the convergencerate, over a fixed finite time T >
0, is faster than O ( N − k ) where N is the order oftruncation [7, 19, 20]. In particular, spectral schemes converge faster than particlemethods when the initial conditions have some degree of regularity. Unfortunatelythe spectral algorithm given above is not qualitatively accurate , as is demonstratedby several examples in Section 8.The goal of this paper is to find a numerical algorithm for (1) and (2) which issimultaneously stable, spectrally convergent, and qualitatively accurate. Previous work.
Within mechanics, spectral methods for the continuum andtransport equation are a common-place where they are viewed as special cases offirst order hyperbolic PDEs [7, 19]. Various Galerkin discretizations of the Koop-man operator have been successfully used for generic dynamic systems [9, 34],most notably fluid systems [39] where such discretizations serve as a generaliza-tion of dynamic mode decomposition [42]. Dually, Ulam-type discretizations of theFrobenius-Perron operator [29, 46] have been used to find invariant manifolds ofsystems with uniform Gaussian noise [16, 17]. In continuous time, Petrov-Galerkindiscretization of the infinitesimal generator of the Frobenius Perron operator con-verge in the presence of noise [6] and preserve positivity in a Haar basis [28].In this article, we consider a unitary representation of the diffeomorphismsof M known to representation theorists [27, 47] and quantum probability theo-rists [33, 36]. To be more specific, we consider the action of diffeomorphisms on To see this compute ddt (cid:104) f, ρ (cid:105) = (cid:82) ( ∂ t f ) ρ + f ( ∂ t ρ ). One finds that the final integral vanishesupon substitution of (1) and (2) and applying integration by parts. The Koopman operator is a linear operator “ K ” which yields the solution f = K · f to (1).We refer the reader to [9] for a survey of recent applications. UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT3 the Hilbert space of half-densities [5, 22]. Half densities can be abstractly summa-rized as an object whose “square” is a density or, alternatively, can be understoodas a mathematician’s nomenclature for a physicist’s “wave functions.” One of thebenefits of working with half-densities, over probability densities, is that the spaceof half-densities is a Hilbert space, while the space of probability densities is aconvex cone [22]. This tactic of inventing the square-root of an abstract objectin order to simplify a problem has been used throughout mathematics. The mostfamilar example would be the invention of the complex numbers to find the rootsof polynomials [43]. A more modern example within applied mathematics can befound in [3] where the (conic) space of positive semi-definite tensor fields whichoccur in non-Newtonian fluids is transformed into the (vector) space of symmetrictensors [3]. Similarly, an alternative notion of half-densities is invoked in [13] totransform the mean-curvature flow PDE into a better behaved one.1.2.
Main contributions.
In this paper we develop numerical schemes for (1) and(2). First, we derive an auxiliary PDE, (8), on the space of half-densities in Section3. We relate solutions of (8) to solutions of (1) and (2) in Theorem 4.1. Second,we pose an auxiliary spectral scheme for (8) in Section 5. Our auxiliary schemeinduces numerical schemes for (1) and (2) via Theorem 4.1. Third, we derivea spectral convergence rate for our auxiliary scheme in Section 6. The spectralconvergence rate for our auxiliary scheme induce spectral convergence rates fornumerical schemes for (1) and (2). Finally, we prove our schemes are qualitativelyaccurate, as in Definition 1.1, in Section 7. We end the paper by demonstratingthese findings in numerical experiments in Section 8. We observe our algorithm for(2) to be superior to a standard spectral Galerkin discretization, both in terms ofnumerical accuracy and qualitative accuracy.1.3.
Notation.
Throughout the paper M denotes a smooth compact n -manifoldwithout boundary. The space of continuous complex valued functions is denoted C ( M ) and has a topology induced by the sup-norm, (cid:107) · (cid:107) ∞ (see [44, 40, 1, 12]).Given a Riemannian metric, g , the resulting Sobolev spaces on M are denoted W k,p ( M ; g ) (see [23]). The tangent bundle to M is denoted by T M , and the n thiterated Whitney sum is denoted by (cid:76) n T M (see [30, 1]). A (complex) densityis viewed as a continuous map ρ : (cid:76) n T M → C which satisfies certain geometricproperties which permit a notion of integration. We denote the space of densitiesby Dens( M ) and the integral of ρ ∈ Dens( M ) is denoted by (cid:82) ρ [30, Chapter16]. By completion of Dens( M ) with respect to the norm (cid:107) ρ (cid:107) := (cid:82) | ρ | we obtain aBanach space, L ( M ). We should note that L ( M ) is homeomorphic to the space ofdistributions up to choosing a partition of unity of M . Given a function f ∈ C ( M ),we denote the multiplication of ρ ∈ L ( M ) by f ρ , and we denote the dual-pairingby (cid:104) f, ρ (cid:105) := (cid:82) f ρ . We let W s, denote the closed subspace of L ( M ) whose elementsexhibit s > H we denote the Banach-algebra of boundedoperators by B ( H ) and topological group of unitary operators by U ( H ). The adjointof an operator L : H → H is denoted by L † . The trace of a trace class operator, L ,is denoted by Tr( L ). The commutator bracket for operators A, B on H is denotedby [ A, B ] := A · B − B · A (see [12]). HENRY O. JACOBS & RAM VASUDEVAN Insights from Operator Theory
Before we describe our algorithms, we take a moment to reflect on the virtueof pursuing qualitative accuracy. If one knows that some entity is conserved un-der evolution, then one can reduce the search for solution by scanning a smallerspace of possibilities. For some, this might be justification enough to proceed aswe are. However, more can be said in this case. The properties that are preservedhave a special relationship to (1) and (2). It is a result known to algebraic ge-ometers, at least implicitly, that algebra-preservation characterizes (1) completelywithout any extra “baggage.” In other words, the only linear evolution PDE thatpreserves the algebra of C ( M ) is (1) . As a corollary, the only linear evolution PDEon densities which preserves duality with functions is of the form (2). Because ofthis fact, many nice properties held by (1) and (2) (such as bounds) are also beheld by a qualitatively accurate integration scheme. Therefore, qualitatively accu-rate schemes leverage the defining aspects of the (1) and (2) to produce numericalapproximations with the same qualitative characteristics.In the remainder of this section, we illustrate how (1) is the unique PDE whichpreserves C ( M ) as an algebra. In the interest of space, we provide references inplace of proofs. To begin, recall the following definitions: Definition 2.1 ([12]) . A Banach-algebra is a Banach space, A , which is equippedwith a multiplication-like operation, “ ( a, b ) ∈ A × A (cid:55)→ ab ∈ A ,” that is bounded,“ (cid:107) ab (cid:107) ≤ (cid:107) a (cid:107)(cid:107) b (cid:107) ,” and associative “ ( ab ) c = a ( bc ) for any a, b, c ∈ A .”A C ∗ -algebra is a Banach algebra, A , over the field C with an involution, “ a ∈ A (cid:55)→ a ∗ ∈ A ,” that satisfies: ( ab ) ∗ = b ∗ a ∗ , ( λa ) ∗ = ¯ λa ∗ , (cid:107) a (cid:107) = (cid:107) a ∗ (cid:107) , (3) for all a, b ∈ A and λ ∈ C . A is unital if A has a multiplicative identity. A is commutative if ab = ba for all a, b ∈ A .Finally, a map T : A → A is called a ∗ -automorphism if T is a bounded linearautomorphism that preserves products, i.e. T ( ab ) = T ( a ) T ( b ) . We denote the spaceof ∗ -automorphisms of A by Aut( A ) . The notion of a C ∗ -algebra may appear abstract so we provide two importantexamples: Example 2.2.
Let X be a topological space. The space of complex valued contin-uous functions with compact support, C ( X ) , is a commutative C ∗ -algebra underthe sup-norm and the standard addition/multiplication/conjugation operations ofcomplex valued functions. If X is compact, then C ( X ) ≡ C ( X ) is unital becausethe constant function, “ f ( x ) = 1 ” is a multiplicative identity. Example 2.3.
For a Hilbert space, H , the space of bounded operators, B ( H ) , isa (non-commutative) C ∗ -algebra under operator multiplication and addition, withthe involution given by the adjoint mapping “ L (cid:55)→ L † ”, and the norm given by theoperator-norm. While the notion of a general C ∗ -algebra is, a priori , more abstract than theexamples above, this feeling of abstraction is an illusion. One of the cornerstonesof operator theory is that all C ∗ -algebras are contained within these examples: Theorem 2.4 (Theorem 1 [18]) . Any C ∗ -algebra is isomorphic to a sub-algebra of B ( H ) for some Hilbert space, H . UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT5
For commutative C ∗ -algebras a stronger result holds if one considers the spaceof character: Definition 2.5 (Space of Characters [21]) . A character of a C ∗ -algebra, A , is anelement of the dual space, x ∈ A ∗ , such that x ( ab ) = x ( a ) x ( b ) for all a, b ∈ A .We denote the space of characters by X A . For each a ∈ A there is a function ˆ a : X A → C given by ˆ a ( x ) = x ( a ) . The map a ∈ A (cid:55)→ ˆ a ∈ C ( X A ) is called the Gelfand Transform . Proposition 2.6 (Lemma 1.1 [21]) . X A ⊂ A ∗ is a compact Hausdorff space withrespect to the relative topology if we impost the weak topology on A ∗ . For example, if A is a space of continuous complex functions on a manifold M ,then X A is the space of Dirac-delta distributions, which is homeomorphic to M itself. The following is a Corollary to 2.4: Theorem 2.7 (Lemma 1 [18]) . Any commutative C ∗ -algebra, A , is canonicallyhomeomorphic to C ( X A ) . The result of Theorem 2.7 is that all commutative C ∗ -algebras are effectivelyrepresented by Example 2.2. Historically, Theorems 2.4 and 2.7 have been valuedbecause they turn abstract C ∗ -algebras (as described by the Definition 2.1) intoExamples 2.2 and 2.3. In this paper, we go the opposite direction. We start withan evolution equation, (1), on the space C ( M ) for a compact manifold M and wetransform it into an equation on a commutative sub-algebra of B ( H ) for a suitablychosen Hilbert space, H by finding an embedding from C ( M ) into B ( H ). Thatwe seemingly transform “concrete” objects into “abstract” objects is one possiblereason that the algorithms in this paper were not constructed earlier. However,“abstract” does not necessarily imply difficult, with respect to numerics. In fact,as this paper shows, it is easier to represent advection equations in this operator-theoretic form. In essence, this is related to a corollary to Theorem 2.7: Corollary 2.8 (follows from Corollary 1.7 of [21]) . Let M be a manifold. If T : C ( M ) → C ( M ) is ∗ -automorphism then there is a unique homeomorphism Φ T ∈ Diff ( M ) such that T [ a ]( x ) = a (Φ T ( x )) . That is, Diff ( X A ) ≡ Aut( A ) as a topological group. Moreover, a linear evolution equation on C ( M ) given by ∂ t f + D [ f ] = 0 for some differential operator, D , preserves the algebra of C ( M ) if and only if D [ f ] = X i ∂f∂x i for some vector-field X . The dual operator is thennecessarily of the form D ∗ [ ρ ] = ∂∂x i ( ρX i ) . Said more plainly, (1) and (2) are the generators of all algebra-preserving au-tomorphisms of C ( M ). Thus, conservation of sums of products is more than justa fundamental property of (1) and (2). Conservation of sums and products is thedefining property of (1) and (2) . As a result, it is natural for this to be reflectedin a discretization . 3. Half densities and other spaces
At the core of any Galerkin scheme, including spectral Galerkin, is the use of aHilbert space upon which everything can be approximated via least squares pro-jections. The methods we present are no exception. In this section, we define a Note: By aiming for qualitative accuracy without sacrificing spectral convergence, we reducethe coefficient of convergence. Therefore this pursuit makes sense from the standpoint of numericalaccuracy as well.
HENRY O. JACOBS & RAM VASUDEVAN canonical L -space associated to a compact manifold M , denoted by L ( M ) forlater use in a spectral discretization . We also define the Sobolev spaces H s ( M ; g )which arise from equipping M with a Riemannian metric, g .For a smooth compact n -manifold, M , let Dens( M ) denote the space of smoothdensities, which we view as anti-symmetric multilinear functions on (cid:76) n T M . Definition 3.1.
A half-density is a smooth complex-valued function ψ : (cid:76) n T M → C such that | ψ | ∈ Dens( M ) . The space of half densities is denoted by (cid:112) Dens( M ) . The following proposition immediately follows from this definition.
Proposition 3.2 (see Appendix A [5]) . If ψ , ψ ∈ (cid:112) Dens( M ) then the scalarproduct ψ · ψ : (cid:76) n T M → C is a complex valued density. This definition is an equivalent reformulation of the half densities defined inthe context of geometric quantization (see [22, Chapter 4] or [5, Appendix A]). Inphysical terms, half densities are a geometric manifestation of the wave functionsused in quantum mechanics. It is unfortunate that physicists call these “wave-functions” given that they are not functions. To test this assertion, observe howelements of (cid:112)
Dens( M ) transform. Under a C -automorphism, Φ : M → M , ahalf density ψ ∈ (cid:112) Dens( M ) transform to a new half density Φ ∗ ψ according to theformula (Φ ∗ ψ ) x ( v , . . . , v n ) := ψ Φ − ( x ) ( D Φ − ( x ) · v , . . . , D Φ − ( x ) · v n )(4)for any x ∈ M and any v , . . . , v n ∈ T x M . This transformation law is inferredby substituting the transformation law for a density into the definition of a half-density. In other words, this is the unique transformation law such that squaringboth sides yields the transformation law for a density. Notably, this is in contrastto the transformation law for functions, which sends f ∈ C ( M ) to the function f ◦ Φ − .In local coordinates, x , . . . , x n , on an open set U ⊂ M , it is common to writea (complex) density ρ : (cid:76) n T M → C as function “ ρ ( x )” for x ∈ U ⊂ M . Thisconvention is permissible as long as one realizes that what is really meant is that ρ ( x ) = f ρ ( x ) | dx ∧ · · · ∧ dx n | for some complex valued function f ρ : U → C . There-fore, when one writes “ ρ ( x ) transforms like ρ (Φ − ( x )) (cid:12)(cid:12) det( D Φ − ( x )) (cid:12)(cid:12) ”, what theyare really describing is how f ρ is transformed. The same notational convention canbe used to represent half-densities locally as “functions” with a different transfor-mation law. In this case the transformation law for half-densities is locally givenby: ˜ ψ ( x ) = (cid:12)(cid:12) det (cid:0) D Φ − ( x ) (cid:1)(cid:12)(cid:12) / ψ (cid:0) Φ − ( x ) (cid:1) . (5)As | ψ | ∈ Dens( M ) for any ψ ∈ (cid:112) Dens( M ), | ψ | can be integrated and weobserve that half densities are naturally equipped with the norm: (cid:107) ψ (cid:107) := (cid:18)(cid:90) M | ψ | (cid:19) / which we call the 2 -norm . We urge the reader familiar with the space L µ ( M ), with respect to some measure µ , notread this section nonetheless. The L -space we use is slightly different, and this fact permeatesthe entire article. UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT7
Definition 3.3. L ( M ) is defined as the completion of (cid:112) Dens( M ) with respect tothe -norm. The space L ( M ) is equipped with a complex inner-product given by (cid:104) ψ | φ (cid:105) = (cid:90) M ¯ ψφ (6) through polar decomposition, and so L ( M ) is a Hilbert space. Lastly, given the transformation law for half-densities, (4) and (5), one candescribe how half-densities are transported by the flow, Φ tX , of the vector field, X . The Lie derivative of a half-density with respect to X is defined as £ X [ ψ ] = − ddt (cid:12)(cid:12) t =0 (Φ tX ) ∗ ψ and is given in local coordinates by: £ X [ ψ ] = 12 X i ∂ψ∂x i + 12 ∂∂x i (cid:0) ψX i (cid:1) . (7)The advection equation can then be written as: ∂ t ψ + £ X [ ψ ] = 0 . (8)Despite the Lie derivative being unbounded, a unique solution is defined for alltime: Proposition 3.4 (Stone’s Theorem [12, 40]) . The unique solution to (8) is of theform ψ ( t ) = U ( t ) · ψ (0) where U ( t ) is the one-parameter semigroup generated by theoperator £ X . Explicitly, U ( t ) is the operator “ (Φ tX ) ∗ ” in the sense that the solutionto (8) is ψ ( t ) = (Φ tX ) ∗ ψ (0) where Φ tX is the time flow map of X at time t .Proof. By inspection we can observe that(Φ tX ) ∗ ( ¯ ψ ψ )) = (Φ tX ) ∗ ψ (Φ tX ) ∗ ψ . By proposition 3.2, ¯ ψ ψ is a density, and so we can integrate it. The integral of adensity is invariant under C transformations [30, Proposition 16.42] and we find0 = ddt (cid:12)(cid:12)(cid:12)(cid:12) t =0 (cid:90) M (Φ tX ) ∗ ( ¯ ψ ψ ) = ddt (cid:12)(cid:12)(cid:12)(cid:12) t =0 (cid:90) M (cid:0) £ X [ ¯ ψ ] ψ ) + ¯ ψ £ X [ ψ ] (cid:1) = (cid:104) £ X [ ψ ] | ψ (cid:105) + (cid:104) ψ | £ X [ ψ ] (cid:105) . Therefore, the operator, £ X is anti-Hermitian. We can see that £ X is denselydefined, as it is well defined on (cid:112) Dens( M ), which is dense in L ( M ) by construc-tion. Stone’s theorem implies that there is a one-to-one correspondence betweendensely defined anti-Hermitian operators on L ( M ) and one-parameter groups U ( t )consisting of unitary operators on L ( M ). Observe that ψ ( t ) = (Φ tX ) ∗ ψ (0) solve(8) directly, by taking its time-derivative. Thus U ( t ) = (Φ tX ) ∗ is the unique one-parameter subgroup we are looking for. (cid:3) The relationship with classical L spaces. To understand the relation-ship to classical Lebesgue spaces, recall that for any manifold M (possibly non-orientable) one can assert the existence of a smooth non-negative reference density µ [30, Chapter 16]. Upon choosing such a µ ∈ Dens( M ), the 2-norm of a continuouscomplex function f : M → C with respect to µ is (cid:107) f (cid:107) µ, = (cid:18)(cid:90) M | f | µ (cid:19) / . (9) HENRY O. JACOBS & RAM VASUDEVAN and L ( M ; µ ) is the completion of the space of continuous functions with respectto this norm. The relationship between L ( M ) and L ( M ; µ ) is that they areequivalent as topological vector-spaces: Proposition 3.5.
Choose a non-vanishing positive density µ : (cid:76) n T M → R + .Let √ µ denote the square root of µ . For any ψ ∈ L ( M ) there exists a unique f ∈ L ( M ; µ ) such that ψ = f · √ µ . This yields an isometry between L ( M ) and L ( M ; µ ) .Proof. It suffices to prove that (cid:112)
Dens( M ) is isomorphic to the space of square inte-grable (w.r.t. µ ) continuous functions, because the later space is dense in L ( M ; µ ).Let ψ ∈ (cid:112) Dens( M ). Then ψ is a continuous density and there exists a uniquefunction g ∈ C ( M ) such that ψ = g · µ . By taking the square root of both sideswe can obtain a unique function f ∈ C such that ψ = f √ µ . The function f isunique with respect to ψ and the map ψ ∈ (cid:112) Dens( M ) (cid:55)→ f ∈ C ( M ; C ) sends (cid:107) · (cid:107) to (cid:107) · (cid:107) µ, by construction. Thus the map is continuous. The inverse of the map isgiven by f ∈ C ( M ; C ) (cid:55)→ f √ µ ∈ (cid:112) Dens( M ). (cid:3) If the spaces are nearly identical the reader may wonder why L ( M ) matters.In fact, the pair are not identical in all aspects. As described earlier, under changeof coordinates or advection, the elements of each space transform differently. Moreimportantly, L ( M ) is not canonically contained within the space of square inte-grable functions, and functions and densities are not contained in L ( M ). Such anembedding may only be obtained by choosing a non-canonical “reference density”,as in Proposition 3.5. This has numerous consequences in terms of what we canand can not do. For example, an operator with domain on L ( M ) can not gen-erally be applied to objects in L ( M ) in the same way. These limitations can behelpful, since they permit vector fields to act differently on objects in L ( M ) thanon objects in L ( M ). These prohibitions serve as safety mechanisms, analogousto the use of overloaded functions in object oriented programs, which due to theirargument type distinctions, effectively banish certain bugs from arising.3.2. Sobolev spaces.
While the “canonicalism” of L ( M ) is useful for this discus-sion, the canonical Sobolev spaces are not. Since the algorithms proposed in thispaper are proven to converge in a Sobolev space, we must still choose a norm andwe rely upon traditional metric dependent definitions. To begin, equip M with aRiemannian metric g : T M ⊕ T M → R . The metric, g , induces a positive density µ g , known as the metric density and an inner-product on C ∞ ( M ) given by: (cid:104) f , f (cid:105) g = (cid:90) f · f µ g . (10)The metric also induces an elliptic operator, known as the Laplace-Beltrami oper-ator ∆ : C ∞ ( M ) → C ∞ ( M ), which is negative-semidefinite (i.e. (cid:82) f ∆ f µ g ≤ f ∈ C ( M )). If M is compact, then L ( M ; µ g ) ∼ = L ( M ) is a separable Hilbertspace and the Helmholtz operator, 1 − ∆, is a positive definite operator with adiscrete spectrum [44]. For any s ≥ Sobolev norm : (cid:107) ψ (cid:107) s, = ( (cid:104) f, (1 − ∆) s · f (cid:105) g ) / (11) Explicitly, if √· : R + → R + is the standard square-root function. Then √ µ := √· ◦ µ : (cid:76) n T M → R + ⊂ C is the half-density which we are considering. UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT9 where f and ψ are related by ψ = f √ µ g . Then we define H s ( M ; g ) as the comple-tion of (cid:112) Dens( M ) with respect to the (cid:107)·(cid:107) s, norm. Such a definition is isomorphic,in the category of topological vector-spaces, to the one provided in [23]. In orderto prove this claim, observe that it holds for bounded sets in R n , and then apply apartition of unity argument to obtain the desired equivalence on manifolds. In par-ticular, note that H ( M ; g ) = L ( M ). It is notable that H s ( M ; g ) as a topologicalvector-space is actually not metric dependent [23, Proposition 2.2]. However, thenorm (cid:107) · (cid:107) s, is metric dependent. Proposition 3.6 (Sobelev Embedding Theorem [44]) . Let ( M, g ) be a compactRiemmanian manifold. If s > t ≥ then H s ( M ; g ) is compactly embedded within H t ( M, g ) .Proof. Let e , e , . . . be the Hilbert basis for L ( M ; µ g ) which diagonalizes ∆ inthe sense that ∆ e i = λ i e i for a sequence 0 = λ ≤ λ ≤ λ ≤ · · · . The operator(1 + ∆) s is given by (1 + ∆) s · f = (cid:88) i e i (1 + λ i ) s (cid:104) e i , f (cid:105) g , (12)and so { (1 + λ i ) − s e i } ∞ i =1 is a Hilbert basis for H s ( M ; g ).Let us call e ( s ) i = (1 + λ i ) − s e i . The embedding of H s ( M ; g ) into H t ( M, g ) isthen given in terms of the respective basis elements by e ( s ) i (cid:55)→ (1 + λ i ) − ( s − t ) e ( t ) i .As s > t and λ i → + ∞ , we see that this embedding is a compact operator [12,Proposition 4.6]. (cid:3) Quantization
In physics, “quantization” refers to the process of substituting certain physicallyrelevant functions with operators on a Hilbert space, while attempting to preservethe symmetries and conservation laws of the classical theory [5, 14, 22]. In thissection, we quantize (1) and (2) by replacing functions and densities with boundedand trace-class operators on L ( M ). This is useful in Section 5 when we discretize.To begin, let us quantize the space of continuous real-valued functions C ( M ).For each f ∈ C ( M ), there is a unique bounded Hermitian operator, H f : L ( M ) → L ( M ) given by scalar multiplication. That is to say ( H f · ψ )( x ) = f ( x ) ψ ( x ) for any ψ ∈ L ( M ). By inspection one can observe that the map “ f (cid:55)→ H f ” is injective andpreserves the algebra of C ( M ) because H f · g + h = H f · H g + H h and (cid:107) H f (cid:107) op = (cid:107) f (cid:107) ∞ .Similarly, (and in the opposite direction) for any trace class operator A there isa unique distribution ρ A ∈ C ( M ) (cid:48) such that: (cid:90) f ρ A = Tr( H † f · A )(13)for any f ∈ C ( M ). More generally, for any A in the dual-space B ( L ( M )) ∗ , thereis a ρ A ∈ C ( M ) (cid:48) such that (cid:104) f, ρ A (cid:105) = (cid:104) H f , A (cid:105) . The map “ A (cid:55)→ ρ A ” is merely theadjoint of the injection “ f (cid:55)→ H f ”. Therefore “ A (cid:55)→ ρ A ” is surjective.We can now convert the evolution PDEs (1) and (2) into ODEs of operators on L ( M ). Theorem 4.1.
Let X ( t ) ∈ X ( M ) be a time-dependent vector-field. Then f satisfies (1) if and only if H f satisfies dH f dt + [ H f , £ X ] = 0 . (14) If A is trace-class and satisfies dAdt + [ A, £ X ] = 0 , (15) then ρ A satisfies (2) . Finally, if ψ satisfies (8) , then ρ A = ψ satisfies (2) and ψ ⊗ ψ † satisfies (15) .Proof. Let f satisfy (1). For an arbitrary ψ ∈ L ( M ) we observe that [ H f , £ X ] · ψ is given in coordinates by:([ H f , £ X ] · ψ ) ( x ) = f ( x ) (cid:18) X j ∂ψ∂x j + 12 ∂∂x j ( ψ X j ) (cid:19) ( x )(16) − X j ∂∂x j (cid:12)(cid:12)(cid:12)(cid:12) x ( f ψ ) + 12 ∂∂x j (cid:12)(cid:12)(cid:12)(cid:12) x ( f ψ X j )(17)where we have used (7). Application of the product rule to each of these termsyields a number of cancellations and we find:[ H f , £ X ] · ψ = − X j ∂f∂x j ψ = ( ∂ t f ) ψ = dH f dt · ψ. (18)As ψ is arbitrary, we have shown that H f satisfies (14). Each line of reasoning isreversible, and so we have proven the converse as well.In order to handle densities note that (cid:104) f, ρ (cid:105) is constant in time when f and ρ satisfy (1) and (2), respectively. By the definition of ρ A , Tr( H † f · A ) = (cid:104) f, ρ A (cid:105) .Therefore: 0 = ddt (cid:16) Tr( H † f · A ) (cid:17) = Tr (cid:32) dH † f dt · A + H † f · dAdt (cid:33) . (19)As was just shown, dH f dt = − [ H f , £ X ] so:0 = Tr (cid:18) − [ H f , £ X ] † · A + H † f · dAdt (cid:19) = Tr( − £ † X H † f · A + H † f £ † X · A + H † f · dAdt ) . (20)Upon noting that £ † X = − £ X and that Tr( abc ) = Tr( bca ):0 = Tr (cid:18) H † f ([ˆ ρ, £ X ] + d ˆ ρdt ) (cid:19) . (21)As H f was chosen arbitrarily, the desired result follows. Again, this line of reasoningis reversible.Lastly, if ψ satisfies (8) and ρ = ψ then we see ∂ t ρ = ∂ t ( ψ ) = 2( ∂ t ψ ) ψ (22) = 2 (cid:18) − ∂ i ( ψX i ) − X i ∂ i ψ (cid:19) ψ = (cid:0) − X i ∂ i ψ − ψ∂ i X i (cid:1) ψ (23) = − ψ ( ∂ i ψ ) X i − ( ∂ i X i ) ψ = − ∂ i ( ρ ) X i − ( ∂ i X i ) ρ = − ∂ i ( ρX i ) . (24) (cid:3) UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT11
The benefit of using (14) and (15) to represent the PDEs of concern is that (14)and (15) may be discretized using a standard least squares projections on L ( M )without sacrificing qualitative accuracy.5. Discretization
This section presents the numerical algorithms for solving (1) and (2). Thebasic ingredient for all the algorithms in this section are a Hilbert basis and anODE solver. Denote a Hilbert basis by { e , e , . . . } for L ( M ). For example,for a Riemannian metric, g , if { f , f , . . . } denote eigen-functions of the Laplaceoperator, then { E k = f k √ µ g | k ∈ N } forms a smooth Hilbert basis for L ( M )where µ g denotes the Riemannian density. We call { E k } the Fourier basis. Toensure convergence, we assume: Assumption 5.1.
Our basis { e k } is such that there exists a metric g for which theunitary transformation which sends the basis { e k } to the Fourier basis is boundedwith respect to the (cid:107) · (cid:107) s, -norm for some s > . In this section we provide a semi-discretization of (1) and (2). Just as a noteto the reader, a “semi-discretization” of the PDE ∂ t φ + F ( φ ) = 0 for some par-tial differential operator, F , is just a discretization of F which converts the PDEinto an ODE [20]. In particular, we assume access to solvers of finite dimensionalODEs, denoted “OdeSolve.” In practice any ODE solver such as Euler’s method,Runge-Kutta, or even well tested software such as [8] could be used to computesuch solutions. Most notably, the method of [10] is specialized to isospectral flowssuch as (14) and (15) by using discrete-time isospectral flows. More explicitly,let OdeSolve( F, x , t ) denote the numerically computed solution x ( t ) to the ODE“ ˙ x = F ( x )” at time t ∈ R , with initial condition x ∈ M . Before constructing analgorithm to spectrally discretize (1) and (2) in a qualitatively accurate manner,we first solve (8) using a standard spectral discretization in Algorithm 1 [7, 38]. Algorithm 1:
A spectral discretization to solve (8) for half densities.
Input : ψ (0) ∈ L ( M ), t ∈ R , N ∈ N .initialize z (0) ∈ C N .;initialize X N ∈ C N × N ; for i = 1 , . . . , N do [ z (0)] i = (cid:82) M ¯ e i ( x ) ψ (0 , x ); for j = 1 , . . . , N do [ X N ] ij = (cid:82) M ¯ e i ( x )( X α ∂ α e j + ∂ α ( X α e j ))( x ) endend initialize the function F : C N → C N given by F ( z ) = X N · z .; z ( t ) = OdeSolve( F, z (0) , t ); Output : ψ N ( t ) = (cid:80) ni =1 [ z ( t )] i e i .To summarize, Algorithm 1 produces a half-density ψ N ( t k ) ∈ V N by projecting(8) to V N . This projection is done by constructing the operator X N = π N ◦ £ X | V N : V N → V N . In Section 6 we prove that ψ N ( t k ) converges to the solution of (8)as N → ∞ . We see that ψ N ( t ) evolves by unitary transformations, just as the exact solution to (8) does. This correspondence is key in providing the qualitativeaccuracy of algorithms that follow, so we formally state it here. Proposition 5.2.
The output of Algorithm 1 is given by U N ( t k ) · ψ N (0) when ψ (0) ∈ L ( M ) is the input to Algorithm 1 where ψ N (0) = π N ( ψ (0)) and U N ( t ) isthe unitary operator as in Proposition 3.4 generated by X N .Proof. The operator X N in Algorithm 1 is anti-Hermitian on V N . It thereforegenerates a unitary action on V N ⊂ L ( M ) when inserted into OdeSolve. (cid:3) Before continuing, we briefly state a sparsity result that aides in selecting a basis.We say an operator A : L ( M ) → L ( M ) is sparse banded diagonal with respectto a Hilbert basis { e , e , . . . } if there exists an integer W ∈ N such that A ( e i ) is afinite sum elements of the form e i + δ j for fewer than W offsets δ j for i = 0 , , , . . . . Theorem 5.3.
Let x , . . . , x n be a dense coordinate chart for M on some denseopen set , then e k = f k √ µ for functions f k ∈ L ( M ; µ ) where µ = | dx ∧ · · · ∧ dx n | (see Proposition 3.5). If ρ ( ∂∂x j ) and H f k are sparse banded diagonal, and if thevector-field X is given in local coordinates by X i = (cid:80) k c ik f k with fewer than W > of c ik ’s being non-zero for each i = 1 , . . . , n , then the matrix X N in Algorithm 1 issparse banded diagonal and the sparsity of X N is O ( W/N ) .Proof. The result follows directly from counting. (cid:3)
Theorem 5.3 suggests selecting a basis where W is small, or at least finite. Forexample, if M were a torus, and the vector-field was made up of a finite number ofsinusoids, then a Fourier basis would yield a W equal to the maximum number ofterms along all dimensions.By Theorem 4.1, the square of the result of Algorithm 1 is a numerical solutionto (2). We can use this to produce a numerical scheme to (2) by finding the squareroot of a density. Given a ρ ∈ Dens( M ), let ρ + denote the positive part and ρ − denote the negative part so that ρ = ρ + − ρ − , then ψ = (cid:112) ρ + − i (cid:112) ρ − is a squareroot of ρ since ρ = ψ . This yields Algorithm 2 to spectrally discretize (2) in aqualitatively accurate manner for densities which admit a square root. Algorithm 2:
A spectral discretization to solve (2) for densities
Data : ρ (0) ∈ L ( M ) , t ∈ R , N ∈ N .Initialize ψ (0) = (cid:112) ρ + (0) − i (cid:112) ρ − (0);Set ψ N ( t ) = Algorithm 1( ψ (0) , t, N ); Output : ρ N ( t, x ) = ψ N ( t, x ) .Alternatively, we could have considered the trace-class operator A N ( t k ) = ψ N ( t k ) ⊗ ψ N ( t k ) † as an output. This would be an numerical solution to (15), and would berelated to our original output in that ρ N ( t k ) = ρ A N ( t k ) . Finally, we present an Such a chart always exists on a compact manifold by choosing a Riemannian metric andextending a Riemannian exponential chart to the cut-locus [41, 26].
UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT13 algorithm to solve (14) (in lieu of solving (1)). This algorithm is presented fortheoretical interest at the moment.
Algorithm 3:
A spectral discretization to solve (14) for functions
Data : f (0) ∈ C ( M ) , t ∈ R , N ∈ N .initialize F N (0) , X n ∈ C N × N .;initialize the linear map B : C N × N → C N × N given by B ( H ) = − [ H, X N ].; for i, j = 1 , . . . , N do [ F N (0)] ij = (cid:82) M ¯ e i ( x ) f ( x ) e j ( x ) ;[ X N ] ij = (cid:82) M ¯ e i ( x )( X α ∂ α e j + ∂ α ( X α e j ))( x ) ; end F N ( t ) = OdeSolve( B, F N (0) , t ); Output : The (compact) operator H f,N ( t k ) = (cid:80) Ni,j =1 [ F N ( t k )] ij e i ⊗ e † j .We find that the ouput of Algorithm 3 bears algebraic similarities similarities tothe exact solution to the infinite dimensional ODE, (14) (which is isomorphic to(1) by Theorem 4.1). This is stated in a proposition analogous to Proposition 5.2. Proposition 5.4. H f,N ( t ) = U N ( t ) · H f,N (0) · U N ( t ) † for any t ∈ R . Moreover, U N ( t ) is identical to the unitary transformation of Proposition 5.2. Lastly, theexact solution of (14) is of the form H f ( t ) = U ( t ) · H f (0) · U ( t ) † as well.Proof. This follows from the fact that algorithm outputs the solution to an isospec-tral flow “ ˙ F N + [ F N , X N ]” where X N is anti-Hermitian and that the H f satisfiesthe isospectral flow (14). (cid:3) Error analysis
In this sections we derive convergence rates. We find that the error boundfor Algorithm 1 induces error bounds for the Algorithms 2 and 3. Therefore, wefirst derive a useful error bound for Algorithm 1. Our proof is a generalization ofthe convergence proof in [37], where (8) is studied (modulo a factor of two timerescaling) on the torus. We begin by proving an approximation bound. In all thatfollows, let π N : L ( M ) → V N denote the orthogonal projection. Proposition 6.1. If ψ ∈ H ¯ s ( M ) and ¯ s > s ≥ , then (cid:107) ψ − π N ( ψ ) (cid:107) s, < d C ¯ s,s ¯ s − s (cid:107) ψ (cid:107) ¯ s, N − s − s ) /n (25) for some constant C ¯ s,s and d = dim( M ) .Proof. We can assume that e , e , . . . is a Fourier basis. The results are unchangedupon applying Assumption 5.1 and converting to the Fourier basis. Any ψ ∈ H s ( M ; g ) can expanded as ψ = ˆ ψ k e k where ˆ ψ k = (cid:104) e k | ψ (cid:105) . As ψ ∈ H s ( M ; g ) itfollows that (cid:107) ψ (cid:107) s, = ∞ (cid:88) k =0 (cid:12)(cid:12)(cid:12) ˆ ψ k (cid:12)(cid:12)(cid:12) (1 + λ k ) ¯ s < ∞ . (26)A corollary of Weyl’s asymptotic formula is that λ k is O ( k /n ) for large k [11,page 155]. After substitution of this asymptotic result into (26) for large k , we see that | ˆ ψ k | is asymptotically dominated by Ck − − s/n for some constant C . Forsufficiently large N we find (cid:107) ψ − π N ( ψ ) (cid:107) s, = (cid:88) k>N (1 + λ k ) s | ˆ ψ k | ≤ C (cid:88) k>N (1 + λ k ) s k − − s/n (27)and by another application of the Weyl formula (cid:107) ψ − π N ( ψ ) (cid:107) s, ≤ ˜ C (cid:88) k>N k s − s ) /n ≤ C s, ¯ s d ¯ s − s N − s − s ) /n . (28)Where the last inequality is derived by bounding the infinite sum with an integral. (cid:3) With this error bound for the approximation error we can derive an error boundfor Algorithm 1:
Theorem 6.2.
Let ψ (0) ∈ H ¯ s ( M ) for ¯ s > s > . Let T > and t ∈ [0 , T ] . Let ψ ( t ) be denote the solution to (8) with initial condition ψ (0) . Finally, let ψ N ( t ) be the output of Algorithm 1 with respect to the inputs ψ (0) , t, N for some N ∈ N .Then the error ε N ( t ) := (cid:107) ψ ( t ) − ψ N ( t ) (cid:107) s, satisfies: ε N ( t ) ≤ (cid:107) ψ (0) (cid:107) ¯ s, K T (cid:18) N − s − t + n ¯ s − s N − s − s ) /n (cid:19) e C T t (29) where K T and C T are positive and constant with respect to N , s , and ¯ s . In particularfor s = (¯ s + 1) / : ε N ( t ) ≤ (cid:107) ψ (0) (cid:107) ¯ s, K T (cid:18) N − ¯ s t + n ¯ s − N (1 − ¯ s ) /n (cid:19) e C T t . (30)To prove Theorem 6.2, we need a perturbed version of Gronwall’s inequality: Lemma 6.3. If dudt ≤ Ku + (cid:15) for some K > then u ( t ) ≤ ( (cid:15)t + u (0)) e Kt .Proof. Let w ( t ) = u ( t ) e − Kt . Then for t ≥ dwdt = dudt e − Kt − Kw ≤ ( Ku + (cid:15) ) e − Kt − Kw = (cid:15)e − Kt ≤ (cid:15) (31)Thus w ( t ) ≤ (cid:15)t + w (0) = (cid:15)t + u (0). (cid:3) Now we can prove Theorem 6.2:
Proof of Theorem 6.2.
Note that dε N dt = ε N (cid:104) ψ − ˜ ψ | (1 + ∆) s ddt ( ψ − ˜ ψ ) (cid:105) By theCauchy-Schwarz inequality dε N dt ≤ (cid:107) £ X [ ψ ] − π N ( £ X [ ψ N ]) (cid:107) s, (32) = (cid:107) £ X [ ψ ] − π N ( £ X [ ψ − ( ψ − ψ N )]) (cid:107) s, . (33)By the triangle inequality and the definition of the operator norm: dε N dt ≤ (cid:107) (1 − π N ) (cid:107) H s − ,op (cid:107) X (cid:107) H s ,op (cid:107) ψ (cid:107) s, + (cid:107) π N (cid:107) op (cid:107) X (cid:107) H s ,op ε N (34)By Proposition 3.4 we observe that ψ ( t ) is related to ψ through the flow of X which is a C k -diffeomorphism if X is C k . From the local expression Proposition UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT15 (cid:107) ψ ( t ) (cid:107) s, is bounded by a scalar multiple of (cid:107) ψ (cid:107) s, .Thus we may write the above bound in the form dε N dt ≤ K (cid:48) (cid:107) − π N (cid:107) H s − ,op (cid:107) ψ (cid:107) s, + C T ε N (35)for constants C T and K (cid:48) . As s >
1, for sufficiently large N we can compute that (cid:107) − π N (cid:107) H s − ,op ≤ (1 + λ N +1 ) − ( s − where λ N denotes the N th eigenvalue ofthe Laplace operator. This is accomplished by observing the operator 1 − π N in aFourier basis and applying to appropriate norms. By Weyl’s asymptotic formula [11,Theorem B.2], λ N asymptotically behaves like N /n . Therefore by Lemma 6.3 with (cid:15) = C T n − s − /d (cid:107) ψ (cid:107) s, : ε N ( t ) ≤ ( K (cid:48) N − s − /n (cid:107) ψ (cid:107) s, t + ε N (0)) e C T t . (36)That ε N (0) behaves as K (cid:48)(cid:48) (cid:107) ψ (cid:107) ¯ s, n − s − s ) /d is a re-statement of Proposition 6.1.We then set K T = max( K (cid:48) , K (cid:48)(cid:48) ). (cid:3) Having derived an error bound for Algorithm 1, we can derive an error boundfor Algorithm 2.
Theorem 6.4.
Let ρ (0) be a distribution in W ¯ s, ( M ) for ¯ s > s > . Let T > and t ∈ [0 , T ] be fixed. Let ρ ( t ) be the solution of (2) at time t . Finally, let ρ N ( t ) be the output of Algorithm 2 with respect to the input ( ρ (0) , t, N ) for some N ∈ N .Then: (cid:107) ρ ( t ) − ρ N ( t ) (cid:107) ≤ (cid:107) ρ (0) (cid:107) ¯ s, K (cid:18) N − s − t + d ¯ s − s N − s − s ) /n (cid:19) e C T t (37) where K is constant with respect to N , and C T is the same constant as in Theorem6.2.Proof. Without loss of generality, assume that ρ is non-negative (otherwise split itinto its non-negative and non-positive components). Let ψ ∈ L ( M ) be such that ρ = ψ , as described in Algorithm 2. It follows that ψ ∈ H s ( M ) and we compute (cid:107) ρ ( t ) − ρ N ( t ) (cid:107) = (cid:90) M | ρ ( t ) − ρ N ( t ) | = (cid:90) M | ψ − ψ N | (38)If we let φ N = ψ − ψ N then we can re-write the above as (cid:107) ρ ( t ) − ρ n ( t ) (cid:107) = (cid:90) M | ψ − ( ψ − φ N ) | = (cid:90) M | ψφ N − φ N | (39) ≤ (cid:107) ψ (cid:107) (cid:107) φ N (cid:107) + (cid:107) φ N (cid:107) = 2 (cid:107) ρ (cid:107) / · (cid:107) φ N (cid:107) + (cid:107) φ N (cid:107) (40)Above we have applied Holder’s inequality to L ( M ), which still holds upon usingthe isometry in Proposition 3.5. Theorem 6.2 provides a bound for (cid:107) φ N (cid:107) . Substi-tution of this bound into the above inequality yields the theorem. (cid:3) Finally, we prove that Algorithm 3 converges to a solution of (14), which isequivalent to a solution of (1) courtesy of Theorem 4.1:
Proposition 6.5.
Let f ∈ C k ( M ) and let H f,N = π N ◦ H f ◦ π N . Then (cid:107) H f − H f,N (cid:107) H s ,op ≤ D ns N − s/n (cid:107) ˆ f (cid:107) op (41) where s > k ≥ , and D is constant. Proof.
Let π ⊥ N = 1 − π N . By Proposition 6.1 we know that (cid:107) π ⊥ N ( ψ ) (cid:107) ≤ ns N − s/n (cid:107) ψ (cid:107) s, (42)for s >
0, then: (cid:107) H f − H f,N (cid:107) H s ,op = sup (cid:107) ψ (cid:107) s, =1 (cid:104) ψ | H f − H f,N | ψ (cid:105) = sup (cid:107) ψ (cid:107) s, =1 (cid:0) (cid:104) ψ | H f | ψ (cid:105) − (cid:104) ψ − π ⊥ N ( ψ ) | H f | ψ − π ⊥ N ( ψ ) (cid:105) (cid:1) = sup (cid:107) ψ (cid:107) s, =1 (cid:0) (cid:60)(cid:104) π ⊥ N ( ψ ) | H f | ψ (cid:105) − (cid:104) π ⊥ N ( ψ ) | H f | π ⊥ N ( ψ ) (cid:105) (cid:1) ≤ sup (cid:107) ψ (cid:107) s, =1 ( (cid:107) π ⊥ N ( ψ ) (cid:107) − (cid:107) π ⊥ N ( ψ ) (cid:107) ) (cid:107) H f (cid:107) op By (42) the result follows. (cid:3)
Theorem 6.6.
Let
T > and t ∈ [0 , T ] be fixed. Let f ( t ) denote the solution to (1) at time t with initial condition f (0) ∈ C k ( M ) . Let H f,N ( t ) denote the outputof Algorithm 3 with respect to the inputs ( f (0) , t, N ) for some N ∈ N . Then: (cid:107) H f ( t ) − H f,N ( t ) (cid:107) H s ,op ≤ D ns N − s/n (cid:107) H f ( t ) (cid:107) op (43) + K T (cid:107) H f,N ( t ) (cid:107) op (cid:18) N − s + 2 ns − N (1 − s ) /n (cid:19) e C T t for the same constant D as in Proposition 6.5 and the same constants C T , K T asin Theorem 6.2.Proof. We find (cid:107) H f ( t ) − H f,N ( t ) (cid:107) H s ,op = sup (cid:107) ψ (cid:107) s, =1 (cid:104) ψ | H f ( t ) − H f,N ( t ) | ψ (cid:105) In light of Proposition 5.4 we find= sup (cid:107) ψ (cid:107) s, =1 (cid:104) ψ | U ( t ) · H f (0) · U ( t ) † − U N ( t ) · H f,N (0) · U N ( t ) † | ψ (cid:105) . The output of Algorithm 3 indicates that H f,N (0) = π ⊥ N ◦ H f (0) ◦ π ⊥ N . Therefore,the above inline equation becomes= sup (cid:107) ψ (cid:107) s, =1 (cid:104) U ( t ) † ψ | H f (0) | U ( t ) † ψ (cid:105) − (cid:104) U N ( t ) † ψ | π ⊥ N ◦ H f (0) ◦ π ⊥ N | U N ( t ) † ψ (cid:105) . and finally= sup (cid:107) ψ (cid:107) s, =1 (cid:104) U ( t ) † ψ | H f (0) − H f,N (0) | U ( t ) † ψ (cid:105) − (cid:104) φ ( t ) | H f,N (0) | φ ( t ) (cid:105) (44)where φ ( t ) = U N ( t ) † ψ − U ( t ) † ψ .The first term is bounded by Proposition 6.5. To bound the second term we mustbound φ . As U N ( t ) † ψ is the backwards time numerical solution to (8) and U ( t ) † ψ is the exact backward time solution to (8), Theorem 6.2 prescribes the existence ofconstants K and C such that: (cid:107) φ (cid:107) s, = (cid:107) U N ( t ) † ψ − U ( t ) † ψ (cid:107) s, ≤ K (cid:107) ψ (cid:107) s, (cid:18) N − s − + ns − s N − s − s ) /n (cid:19) e Ct UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT17 for any s < s . This expression can be simplified by noting that (cid:107) ψ (cid:107) s, = 1, setting s = (1 + s ) /
2, and noting that the H s norm is stronger than the L -norm to get: (cid:107) φ (cid:107) ≤ K (cid:18) N − s + 2 ns − N (1 − s ) /n (cid:19) e Ct . By applying the Cauchy-Schwarz inequality to (44) and our derived bound on φ : (cid:107) H f ( t ) − H f,N ( t ) (cid:107) H s ,op ≤ (cid:107) H f (0) − H f,N (0) (cid:107) H s ,op + K (cid:107) H f,N (cid:107) op (cid:18) N − s + 2 ns − N (1 − s ) /n (cid:19) e Ct Upon invoking Proposition 6.5 we get the desired result. (cid:3) Qualitative Accuracy
In this section, we prove that our numerical schemes are qualitatively accurate.We begin by illustrating the preservation of appropriate norms. Throughout thissection let ψ N ( t ), ρ N ( t ), and H f,N ( t ) denote the sequence of outputs of Algorithms1, 2, and 3 with respect to initial conditions ψ (0) ∈ H s ( M ; g ) , ρ (0) ∈ W s, and f (0) ∈ C s ( M ) for N = 1 , , . . . . Theorem 7.1.
Let ψ, ρ, f denote solutions to (8) , (2) , and (1) respectively. Let ψ N ( t ) , ρ N ( t ) , and H f,N ( t ) , denote outputs from algorithms 1, 2, and 3 respectivelyfor a time t < ∞ . Then (cid:107) ψ N ( t ) (cid:107) , (cid:107) ρ N ( t ) (cid:107) , and (cid:107) H f,N ( t ) (cid:107) op are constant withrespect to t for arbitrary N ∈ N . Moreover, lim N →∞ (cid:107) ψ N ( t ) (cid:107) = (cid:107) ψ ( · ; t ) (cid:107) , lim N →∞ (cid:107) ρ N ( t ) (cid:107) nuc = (cid:107) ρ ( · ; t ) (cid:107) , lim N →∞ (cid:107) H f,N ( t ) (cid:107) op = (cid:107) f ( · ; t ) (cid:107) sup . Proof.
To prove (cid:107) H f,N (cid:107) op is conserved note that the evolution is isospectral [10].We have already shown that H f,N ( t ) converges to H f ( t ) in the operator norm.Convergence of the norms follows from the fact that (cid:107) H f ( t ) (cid:107) op = (cid:107) f (cid:107) sup . Anidentical approach is able to prove the desired properties for ρ N ( t ) and ψ N ( t ) aswell. (cid:3) Theorem 7.1 is valuable because each of the norms is naturally associated to theentity which it bounds, and these quantities are conserved for the PDEs that thispaper approximates. For example, (cid:107) H f (cid:107) op = (cid:107) f (cid:107) sup for a function f , and this isconstant in time when f is a solution to (1). A discretization constructed accordingto Algorithm 3 according to Theorem 7.1 is constant for any N , no matter howsmall.The full Banach algebra C ( M ) is conserved by advection too. This property isencoded in our discretization as well. Theorem 7.2.
Let f ( x ; t ) , g ( x ; t ) , and h ( x ; t ) be solutions of (1) and let k = f · g + h .Let H f,N , H g,N and H h,N be numerical solutions constructed by Algorithm 3, then H k,N ( t ) = H f,N · H g,N + H h,N satisfies ddt H k,N = [ X N , H k,N ] . (45) Moreover, H k,N ( t ) strongly converges to H k as N → ∞ in the operator norm on H s ( M ) when f, g, h ∈ C s ( M ) for s > . Proof.
By construction, the output of Algorithm 3 is the result of an isospectralflow, and is therefore of the form H f,N ( t ) = U N ( t ) H f,N (0) U N ( t ) † (46) H g,N ( t ) = U N ( t ) H g,N (0) U N ( t ) † (47) H h,N ( t ) = U N ( t ) H h,N (0) U N ( t ) † . (48)We then observe H k,N ( t ) = U N ( t ) H k,N (0) U N ( t ) † = U ( t ) ( H f,N (0) H g,N (0) + H h,N (0)) U ( t ) † (49) = U ( t ) H f,N (0) U ( t ) † U ( t ) H g,N (0) U ( t ) † + U ( t ) H h,N (0) U ( t ) † (50) = H f,N ( t ) H g,N ( t ) + H h,N ( t ) . (51)Differentiation in time implies the desired result. Convergence follows from Theo-rem 6.6. (cid:3) Finally, the duality between functions and densities is preserved by advection.If f satisfies (1) and ρ satisfies (2) then (cid:82) f ρ is conserved in time. Algorithms 2and 3 satisfy this same equality: Theorem 7.3.
For each N ∈ N , Tr( H f,N A ρ,N ( t )) is constant in time where A ρ,N ( t ) = ψ N ( t ) ⊗ ψ N ( t ) † . Moreover, Tr( H f,N A ρ,N ) converges to the constant (cid:82) f ρ as N → ∞ .Proof. As H f,N ( t ) = U N ( t ) H f,N (0) U N ( t ) † and ψ N ( t ) = U N · ψ N (0) we observe thatTr( H f,N ( t )( ψ N ( t ) ⊗ ψ † N ( t ))) = Tr( U N ( t ) H f,N (0) U N ( t ) † U ) N ( t )( ψ N (0) ⊗ ψ N (0) † ) U N ( t ) † )(52) = Tr( H f,N (0)( ψ N (0) ⊗ ψ N (0) † ))(53)Convergence follows from Theorems 6.6 and 6.4. (cid:3) Numerical Experiments
This section describes two numerical experiments. First, a benchmark compu-tation to illustrate the spectral convergence of our method and the conservationproperties in the case of a known solution is considered.8.1.
Benchmark computation.
Consider the vector field ˙ x = − sin(2 x ) for x ∈ S . The flow of this system is given by:Φ tX ( x ) = atan (cid:0) e t tan( x ) (cid:1) . (54)If the initial density is a uniform distribution, ρ = dx , then the the exact solutionof (2) is: ρ ( x ; t ) = (cid:0) e t sin ( x ) + e − t cos ( x ) (cid:1) − | dx | (55)Figure 1 depicts the evolution of ρ ( x ; t ) at t = 1 . UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT19 (a)
Exact (b)
Standard spectral (c)
Algorithm 2
Figure 1.
A benchmark illustration of Algorithm 2 on the exam-ple described in Section 8.1.Here we witness how Algorithm 2 has greater qualitative accuracy than a stan-dard spectral discretization, in the “soft” sense of qualitative accuracy. For exam-ple, standard spectral discretization exhibits negative mass, which is not achievablein the exact system. Moreover, the L -norm is not conserved in standard spectraldiscretization. In contrast, Theorem 7.1 proves that the L -norm is conserved byAlgorithm 2. A plot of the L -norm is given in Figure 2. Finally, a convergence plotis depicted in Figure 3. Note the spectral convergence of Algorithm 2. In terms ofnumerical accuracy, Algorithm 2 appears to have a lower coefficient of convergence. time L n o r m Figure 2.
A plot of the L -norm vs time of a standard spectraldiscretization (solid) and the result of Algorithm 2 (dotted) on theexample described in Section 8.1.In general, Algorithm 3 is very difficult to work with, as it outputs an operatorrather than a classical function. However, Algorithm 3 is of theoretical value, inthat it may inspire new ways of discretization (in particular, if one is only interestedin a few level sets). We do not investigate this potentiality here in the interest offocusing on the qualitative aspects of this discretization. For example, under theinitial conditions g ( x ) = cos( x ) and f = sin( x ) the exact solutions to (14) are: g ( x, t ) = cos( x ) (cid:0) e t sin ( x ) + cos ( x ) (cid:1) − / f ( x, t ) = sin( x ) (cid:0) sin ( x ) + e − t cos ( x ) (cid:1) − /
20 HENRY O. JACOBS & RAM VASUDEVAN -9 -8 -7 -6 -5 -4 -3 -2 -1 Figure 3.
Convergence plot for Algorithm 2 (dotted) and a stan-dard spectral method (solid) in the L -norm.Under the initial condition h = f · g = sin( x ) cos( x ) the exact solution to (14) is: h ( x, t ) = f ( x, t ) g ( x, t ) = cos( x ) sin( x ) (cid:0) cos ( x ) + e t sin ( x ) (cid:1) − . One can compute h by first multiplying the initial conditions and then using Al-gorithm 3 to evolve in time, or we may evolve each initial condition in time first,and multiply the outputs. If one uses Algorithm 3, then both options, as a resultof Theorem 7.2, yield the same result up to time discretization error (which is ob-tained with error tolerance 1 e − L ( M ). As proven by Theorem 7.1, the operator-normis conserved by Algorithm 3. In contrast, the sup-norm drifts over time under astandard discretization. This is depicted in Figure 5 Figure 4.
The discrepancy due to non-preservation of scalarproducts under a standard spectral Galerkin discretization. Thediscrepancy of Algorithm 3 (not plotted) is attributable to ourtime-discretization scheme where we only tolerated error of 10 − in this instance. UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT21 s u p / o p e r a t o r n o r m Figure 5.
A plot of the sup-norm vs time of a standard spectraldiscretization (blue) and the result of Algorithm 3 (red) on theexample described in Section 8.1.8.2.
A modified ABC flow.
Consider the system˙ x = A sin(2 πz ) + C cos(2 πy ) + D cos(2 πx )(56) ˙ y = B sin(2 πz ) + A cos(2 πy ) + D cos(2 πy )(57) ˙ z = A sin(2 πz ) + B cos(2 πy ) + D cos(2 πz )(58)on the three-torus for constants A, B, C, D ∈ R . When D = 0 this system is the wellstudied volume conserving system known as an Arnold-Beltrami-Childress flow [2].When A > B > C > D = 0, and C <<
1, then the solutions to this ODE arechaotic, with a uniform steady state distribution [32]. When D = 0 the operator £ X of (8) is identical to the operator ∂ α ( ρX α ) that appears in (2), and Algorithm1 do not differ from a standard spectral discretization. Therefore we consider thecase where
D >
D > A = 1 . , B = 0 . , C = 0 . , and D = 0 .
5. As an initial condition consider a wrapped Gaussian distribution withanisotropic variance σ = (0 . , . , .
3) centered at (0 , , z -marginal of these densities are illustrated in Figure 6. The toprow depicts the results from using Algorithm 2 using 33 modes along each dimen-sion. The middle row depicts the results from using a Monte-Carlo method with15 = 3375 particles as a benchmark computation. Finally, the bottom row depictsthe results from using a standard Fourier based discretization of (2) using 33 modesalong each dimension. Notice that Algorithm 2 performs well when compared tothe standard discretization approach. This is always the case for a volume conserving system.
Figure 6.
An illustration of the performance of Algorithm 2 (toprow), Monte Carlo (middle row) and a standard spectral Galerkin(bottom row) on the example described in Section 8.2. The domainis the 2-torus. Here we’ve consider an initial probability densitygiven by a wrapped Gaussian. Darker regions represent areas ofhigher-density. 9.
Conclusion
In this paper we constructed a numerical scheme for (1) and (2) that is spec-trally convergent and qualitatively accurate, in the sense that natural invariants arepreserved. The result of obeying such conservation laws is a robustly well-behavednumerical scheme at a variety of resolutions where legacy spectral methods fail.This claim was verified in a series of numerical experiments which directly com-pared our algorithms with standard Fourier spectral algorithms. The importanceof these conservation laws was addressed in a short discussion on the Gelfand Trans-form. We found that conservation laws completely characterize (1) and (2), and thisexplains the benefits of using qualitatively accurate scheme at a more fundamentallevel.9.1.
Acknowledgements.
This paper developed over the course of years fromdiscussions with many people whom we would like to thank: Jaap Eldering, GaryFroyland, Darryl Holm, Peter Koltai, Stephen Marsland, Igor Mezic, Peter Michor,Dmitry Pavlov, Tilak Ratnanather, and Stefan Sommer. This research was madepossible by funding from the University of Michigan.
References [1]
R Abraham, J E Marsden, and T S Ratiu , Manifolds, Tensor Analysis, and Applications ,vol. 75 of Applied Mathematical Sciences, Spinger, 3rd ed., 2009.[2]
V. I. Arnold and B. A. Khesin , Topological Methods in Hydrodynamics , vol. 24 of AppliedMathematical Sciences, Springer Verlag, 1992.[3]
Nusret Balci, Becca Thomases, Michael Renardy, and Charles R Doering , Sym-metric factorization of the conformation tensor in viscoelastic fluid models , Journal of Non-Newtonian Fluid Mechanics, 166 (2011), pp. 546–553.
UALITATIVELY ACCURATE SPECTRAL SCHEMES FOR ADVECTION AND TRANSPORT23 [4]
G. K. Batchelor , An introduction to fluid dynamics , Cambridge Mathematical Library,Cambridge University Press, Cambridge, paperback ed., 1999.[5]
S. Bates and A. Weinstein , Lectures on the geometry of quantization , vol. 8 of BerkeleyMathematics Lecture Notes, American Mathematical Society, Providence, RI, 1997.[6]
Andreas Bittracher, P´eter Koltai, and Oliver Junge , Pseudogenerators of spatialtransfer operators , SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1478–1517.[7]
John P. Boyd , Chebyshev and Fourier spectral methods , Dover Publications, Inc., Mineola,NY, second ed., 2001.[8]
Peter N. Brown, George D. Byrne, and Alan C. Hindmarsh , VODE: a variable-coefficient ODE solver , SIAM J. Sci. Statist. Comput., 10 (1989), pp. 1038–1051.[9]
Marko Budiˇsi´c, Ryan Mohr, and Igor Mezi´c , Applied koopmanism , Chaos: An Interdis-ciplinary Journal of Nonlinear Science, 22 (2012), pp. –.[10]
Mari Calvo, Arieh Iserles, and Antonella Zanna , Numerical solution of isospectralflows , Mathematics of Computation of the American Mathematical Society, 66 (1997),pp. 1461–1486.[11]
Isaac Chavel , Eigenvalues in Riemannian geometry , vol. 115 of Pure and Applied Math-ematics, Academic Press, Inc., Orlando, FL, 1984. Including a chapter by Burton Randol,With an appendix by Jozef Dodziuk.[12]
John B. Conway , A course in functional analysis , vol. 96 of Graduate Texts in Mathematics,Springer-Verlag, New York, second ed., 1990.[13]
Keenan Crane, Ulrich Pinkall, and Peter Schr¨oder , Robust fairing via conformalcurvature flow , ACM Trans. Graph., 32 (2013).[14]
Paul AM Dirac , Lectures on quantum mechanics , Courier Corporation, 2013.[15]
Lawrence C. Evans , Partial differential equations , vol. 19 of Graduate Studies in Mathe-matics, American Mathematical Society, Providence, RI, second ed., 2010.[16]
G. Froyland, O. Junge, and P. Koltai , Estimating long-term behavior of flows without tra-jectory integration: the infinitesimal generator approach , SIAM J. Numer. Anal., 51 (2013),pp. 223–247.[17]
G. Froyland and K. Padberg , Almost-invariant sets and invariant manifolds—connectingprobabilistic and geometric descriptions of coherent structures in flows , Phys. D, 238 (2009),pp. 1507–1523.[18]
I. Gelfand and M. Neumark , On the imbedding of normed rings into the ring of operatorsin Hilbert space , Rec. Math. [Mat. Sbornik] N.S., 12(54) (1943), pp. 197–213.[19]
D. Gottlieb and J.S. Hesthaven , Spectral methods for hyperbolic problems , Journal ofComputational and Applied Mathematics, 128 (2001), pp. 83 – 131. Numerical Analysis2000. Vol. VII: Partial Differential Equations.[20]
David Gottlieb and Steven A Orszag , Numerical analysis of spectral methods: theoryand applications , vol. 26, Siam, 1977.[21]
Jos´e M. Gracia-Bond´ıa, Joseph C. V´arilly, and H´ector Figueroa , Elements of noncom-mutative geometry , Birkh¨auser Advanced Texts: Basler Lehrb¨ucher. [Birkh¨auser AdvancedTexts: Basel Textbooks], Birkh¨auser Boston, Inc., Boston, MA, 2001.[22]
V Guillemin and S Sternberg , Geometric Asymptotics , vol. 14 of Mathematical Surveysand Monographs, American Mathematical Society, 1970.[23]
Emmanuel Hebey , Nonlinear analysis on manifolds: Sobolev spaces and inequalities , vol. 5of Courant Lecture Notes in Mathematics, New York University, Courant Institute of Math-ematical Sciences, New York; American Mathematical Society, Providence, RI, 1999.[24]
D Henrion and M Korda , Convex computation of the region of attraction of polynomialcontrol systems , IEEE Transactions on Automatic Control, 59 (2014), pp. 297–312.[25]
Lars H¨ormander , The analysis of linear partial differential operators. I , Classics in Math-ematics, Springer-Verlag, Berlin, 2003. Distribution theory and Fourier analysis, Reprint ofthe second (1990) edition [Springer, Berlin; MR1065993 (91m:35001a)].[26]
Richard Montgomery (http://mathoverflow.net/users/2906/richard mont-gomery) , Does every compact manifold exhibit an almost global chart . MathOverflow.URL:http://mathoverflow.net/q/177913 (version: 2014-08-06).[27]
R. S. Ismagilov , The unitary representations of the group of diffeomorphisms of the space R n , n ≥
2, Mat. Sb. (N.S.), 98(104) (1975), pp. 55–71, 157–158.[28]
P´eter Koltai , Efficient approximation methods for the global long-term behavior of dynam-ical systems: theory, algorithms and examples , Logos Verlag Berlin GmbH, 2011. [29]
A. Lasota and M. C. Mackey , Chaos, Fractals, and Noise , Applied Mathematical Sciences,Springer Verlag, 1994.[30]
John M Lee , Introduction to smooth manifolds , vol. 218 of Graduate Texts in Mathematics,Springer-Verlag, 2nd ed., 2006.[31]
Randall J. LeVeque , Numerical methods for conservation laws , Lectures in MathematicsETH Z¨urich, Birkh¨auser Verlag, Basel, second ed., 1992.[32]
Andrew J. Majda and Andrea L. Bertozzi , Vorticity and incompressible flow , vol. 27 ofCambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2002.[33]
P. A. Meyer , Quantum probability for probabilists , vol. 1538 of Lecture Notes in Mathemat-ics, Springer-Verlag, Berlin, 1993.[34]
I. Mezi´c , Spectral properties of dynamical systems, model reduction and decompositions ,Nonlinear Dynam., 41 (2005), pp. 309–325.[35]
Bernt Øksendal , Stochastic differential equations , Universitext, Springer-Verlag, Berlin,sixth ed., 2003. An introduction with applications.[36]
Kalyanapuram R Parthasarathy , An introduction to quantum stochastic calculus ,Springer Science & Business Media, 2012.[37]
Joseph E Pasciak , Spectral and pseudospectral methods for advection equations , Mathemat-ics of Computation, 35 (1980), pp. 1081–1092.[38]
William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flan-nery , Numerical recipes , Cambridge University Press, Cambridge, third ed., 2007. The artof scientific computing.[39]
Clarence W Rowley, Igor Mezi´c, Shervin Bagheri, Philipp Schlatter, and Dan SHenningson , Spectral analysis of nonlinear flows , Journal of Fluid Mechanics, 641 (2009),pp. 115–127.[40]
Walter Rudin , Functional analysis. international series in pure and applied mathematics ,1991.[41]
Takashi Sakai , Riemannian geometry , vol. 149 of Translations of Mathematical Monographs,American Mathematical Society, Providence, RI, 1996. Translated from the 1992 Japaneseoriginal by the author.[42]
P J Schmid , Dynamic mode decomposition of numerical and experimental data , Journal ofFluid Mechanics, 656 (2010), pp. 5–28.[43]
Ian Stewart , Galois theory , CRC Press, Boca Raton, FL, fourth ed., 2015.[44]
Michael Taylor , Pseudo differential operators , Lecture Notes in Mathematics, Vol. 416,Springer-Verlag, Berlin-New York, 1974.[45]
C. Truesdell , A First Course in Rational Continuum Mechanics: General Concepts , Aca-demic Press, 1991.[46]
Stanislaw M Ulam and John Von Neumann , On combination of stochastic and determin-istic processes-preliminary report , Bulletin of the American Mathematical Society, 53 (1947),pp. 1120–1120.[47]