A study of an approximate equation for sharp fronts in the generalized Surface Quasi-Geostrophic Equations in the cylinder

(1)

U.U.D.M. Project Report 2019:6

Examensarbete i matematik, 30 hp Handledare: Jordi-Lluís Figueras Examinator: Denis Gaidashev Februari 2019

A study of an approximate equation for sharp fronts in the generalized Surface

Quasi-Geostrophic Equations in the cylinder

Marc Fraile

(2)

(3)

Introduction

1.1 The physical setting

The Surface Quasi-Geostrophic Equation (1.5) is a pseudo-differential equation that originates in the field of geophysics. For an in-depth reference of all the involved concepts see [1]; we will give a quick overview in this section.

In fluid mechanics, focus is given to physical quantities that move along a fluid medium. Let u be the velocity of such a medium: given a physical quantity θ that flows with u, its time evolution is more easily described through the material derivative

Dθ Dt = ∂θ

∂t + u · ∇θ.

When studying atmospherical phenomena dominated by the Earth’s spin, it is sometimes good enough to ignore any vertical motion of the air and describe its winds as a 2-dimensional flow on the surface of a sphere. Using polar coordinates (x, y) ∈ R², where x represents longitude and y latitude, the wind velocity field u can be described as a R² → R² function with components (u, v). This coordinate system is degenerate at the poles, but works well in the mid-latitudes where interesting phenomena like the jet streams happen. Under certain approximations, u can be taken as the perpendicular gradient

∇^⊥ := (−∂y, ∂x) of a stream function ψ:

u = ∇^⊥ψ = (−ψy, ψx).

That is, u flows along the contour lines of ψ. Two interesting cases arise here: First, under the full 3-dimensional setting, the vorticity ω of the flow u is defined as its rotational ω = ∇ × u. When reduced to a 2-dimensional setting, this gets reduced to the z-component:

ω = uy− vx= −ψyy− ψxx= −∆ψ.

In some cases, one can assume ω to be conserved, that is, Dω

Dt = 0, (1.1a)

ω = −∆ψ. (1.1b)

Equations (1.1) are known as the stream-vorticity formulation of the incompressible Euler Equations, as written using Physics notation. Assuming the Laplacian operator can be inverted (that is, the choice of ψ is restricted enough that we can uniquely determine it given a valid ω) and using a more mathematically inclined notation, this becomes:

(7)

Figure 1.1: Jet streams are self-sustaining stripes of high velocity wind in the surface of a planet, caused by the interaction between pressure differences and the Coriolis force. This figure shows the Northern Hemisphere jet stream (pink) on Earth, fueled by the temperature difference between cold polar air (blue) and warmer air (orange). It shows meanders (Rossby Waves) developing (a), (b); then finally detaching a ”drop” of cold air (c). (CC) Wikimedia Commons.

Definition 1.1 (Euler Equations). The Euler Equations are the following set of equations, where ω : R²→ R and u : R²→ R²:

ωt+ u · ∇ω = 0, u = ∇^⊥(−∆)⁻¹ω.

(1.2a) (1.2b) Second, in geophysics a useful quantity is the (quasi-geostrophic) potential vorticity θ, which again can be described as a real-valued function on the sphere. Under a set of approximations which are collectively described as quasi-geostrophic flow, valid for mid-latitudes, θ is conserved (Dθ/Dt = 0), and a peculiar relationship appears between the Fourier transforms of θ and ψ. Formally considering ψ as a plane wave ψ(x) = ˆψ(k)e^ik·xfor some k ∈ R²yields θ as a plane wave θ(x) = ˆθ(k)e^ik·x, with the defining relationship

ψ(k) = −ˆ 1

|k|θ(k).ˆ (1.3)

This is briefly explained in [2] and [3]. It establishes an algebraic relation between the Fourier transforms ψ and ˆˆ θ, which a priori lets us write a conservation law without ψ, simiilar to the substitution we did for (1.2). To describe this law, which we will find to be very similar to (1.2), we need to discuss the fractional Laplacian.

1.2 Fractional Laplacian and Fourier multipliers

1.2.1 General case

Given f ∈ L¹(Rⁿ) ∩ C¹(Rⁿ) and 1 ≤ j ≤ n, assume ∂x_jf ∈ L¹(Rⁿ), and further assume that F [f ] and F [∂x_jf ] are integrable functions. Then

∂x_jf (x) = ∂x_j

Z

Rⁿ

F [f ](k)e^ix·kdk = Z

Rⁿ

F [f ](k)(ikj)e^ix·kdk ∀x ∈ Rⁿ,

∂xjf (x) = Z

Rⁿ

F [∂xjf ](k)e^ix·kdk ∀x ∈ Rⁿ. The uniqueness of the Fourier transform in L¹(Rⁿ) then implies

F [∂x f ] = ikjF [f ].

(8)

Similar operations reveal similar relations. Using multiindices of the form α = (α₁, . . . , α_n) ∈ Nⁿ, and writing derivatives as D^α= (i∂_x₁)^α¹· · · (i∂_x_n)^αⁿ, any linear differential operator of order n with constant coefficients can be written as

P = X

|α|≤n

aαD^α, aα∈ C,

where |a| =Pn

k=1a_k. Using similar reasoning as above, one can then prove that if, all involved functions are regular enough,

F [P f ](k) = p(k)F [f ](k), where p : Rⁿ→ C is the complex polynomial given by

p(k) = X

|α|≤n

aαk^α,

and k^α := k₁^α¹· · · k^α_nⁿ. This means many differential operators can be treated through the Fourier transform as multiplicative functions. As an important particular example, the Laplacian becomes

F [∆f ](k) = −|k|²F [f ](k).

If we ask f to be a zero-mean function, F [f ](0) = 0, we can invert this formula: if g = ∆f and both functions are regular enough, we can uniquely determine f from g through

F [f ](k) = − 1

|k|²F [g](k) ∀k 6= 0, F [f ](0) = 0.

This algebraic relation allows us to define the inverse Laplacian (−∆)⁻¹, and more generally a family of operators which do not correspond to classical differential equations, the fractional Laplacian (−∆)^s, s ∈ R, through the relation

F [(−∆)^sf ](k) = |k|^2sF [f ](k).

To be specific, one needs the right hand side of this equation to be an element of L¹(Rⁿ), so that the formula can be reversed:

Definition 1.2 (Fractional Laplacian). Given s ∈ R, let A be the set of all functions f ∈ L¹(Rⁿ) such that |k|^2sF [f ](k) ∈ L¹(Rⁿ), and let m(k) := |k|^2s. The fractional Laplacian (−∆)^s is defined as the operator A → L^∞(Rⁿ) given by

(−∆)^sf = F⁻¹(mF [f ]) .

This construction can be further generalized to a much wider class of operators. As described in [4] and outlined in Appendix A, the Fourier transform can be defined on functions f : G → C for any locally compact Abelian group (LCA group) G, which is associated to a dual LCA group Γ (unique up to isomorphism), and ˆf is a function from Γ to C. With this in mind,

Definition 1.3 (Fourier multiplier). Given a LCA group G, denote its dual as Γ. Given a function m : Γ → C, define the set A^mas

A_m:= {f ∈ L¹(G) : m · F [f ] ∈ L¹(Γ)}.

The Fourier multiplier T_m: A_m→ L^∞(G) is the operator Tmf = F⁻¹(mF [f ]),

and m is called the symbol of Tm. When no confusion is caused, Tm is just called a multiplier.

(9)

Examples 1.1.

1. Given a ∈ C, the operator f 7→ af is a multiplier with symbol m(k) = a.

2. For G = R or G = T, the derivative operator f 7→ f⁰ is a multiplier with symbol m(k) = ik.

3. For G = Rⁿ and 1 ≤ j ≤ n, the partial derivative operator ∂_x_j is a multiplier with symbol m(k) = ik_j.

4. The Laplacian ∆ is a multiplier with symbol m(k) = −|k|².

5. Given s ∈ R, the fractional Laplacian (−∆)^s is a multiplier with symbol m(k) = |k|^2s. 6. The operator |∂x| := (−∆)^1/2 is a multiplier with symbol m(k) = |k|.

7. The operator log |∂x| is defined as the multiplier over T with symbol m(0) = 0, m(k) = log |k| for k 6= 0.

If m never vanishes, Tmcan be easily inverted:

ˆ

g = m ˆf ⇔ f =ˆ 1

mg.ˆ (1.4)

One must be careful with the domain of definition, since g might not be in L¹(G). To formalize this:

Lemma 1.1 (Fourier multiplier invertibility). Let m : Γ → C^∗. Let Bm = T_m⁻¹(L¹(G) ∩ A_1/m), Cm = Tm(Bm). Then T1

m|C_m is the inverse of Tm|B_m.

Proof. By definition, Cm ⊆ L¹(G), and if g ∈ Cm there is f ∈ Bm such that g = F⁻¹(mF [f ]). Since mF [f ] ∈ L¹(Γ), the Fourier inversion theorem applies, so F [g] = mF [f ]. Also by definition, Cm⊆ A1/m, so (1/m)F [g] = F [f ] is integrable and we can apply the Fourier inversion theorem. Thus, T¹

mg = F⁻¹(F [f ]) = f .

1.2.2 The particular case L

²

(T)

In the case G = T (Γ = Z), restricting Fourier multipliers to T^m: L²(T) → L²(T) lets us use the Hilbert space structure of L²(T) to our advantadge. In particular, it lets us separate neatly the kernel of an operator from its orthogonal complement, and use the relation (1.4) to characterize any multiplier.

Lemma 1.2 (Fourier multipliers in L²(T)). Given m : Z → C, consider the restricted set D_m:= {f ∈ L²(T) : m · F[f ] ∈ L¹(Z)}.

Then Dm⊆ Am, and Tm(Dm) ⊆ L²(T). When dealing with L²(T) we will assume Fourier multipliers to be restricted to Tm: Dm→ L²(T).

Proof. Since T has finite measure, L²(T) ⊂ L¹(T), proving Dm⊆ Am. Since Z is discrete, the opposite relation holds: L¹(Z) ⊂ L²(Z), so if f ∈ D_mthen m · F [f ] ∈ L²(Z). As pointed out in Appendix A, this implies F⁻¹(m · F [f ]) ∈ L²(T).

Lemma 1.3 (Fourier multiplier kernels and inverses). Let m : Z → C. Denote Z := {k ∈ Z : m(k) = 0} ⊆ Z.

Then ker T_m= hφ_ki_k∈Z. Furthermore, T_m|_{(ker T}_m₎⊥ is a bijection into T_m(D_m), with inverse T_m⁻¹g =X

k /∈Z

1

m(k)ˆgkφk, g ∈ Tm(Dm).

Finally, Tm(Dm) ⊆ hφki _∈Z, and {φk} _∈Z ⊂ Tm(Dm).

(10)

Proof.

1. Clearly hφkik∈Z ⊆ ker Tm, since Tmφk = 0 for all k ∈ Z (in particular, hφkik∈Z ⊆ Dm). Consider now f ∈ D_msuch that T_mf = 0. It follows that m(k) ˆf_k = 0 for all k ∈ Z, so ˆf_k can only be nonzero if k ∈ Z, proving ker T_m⊆ hφkik∈Z.

2. Given g ∈ Tm(Dm), it follows that g ∈ L²(T) and F[g] = mF[f ] for some f ∈ D^m, so f ∈ L²(Z) and F [f ] ∈ L²(Z). Consider now the projection of f into (ker T^m)^⊥, f^⊥. It follows that f = f^⊥+ f^k, with f^k∈ ker Tm, so Tmf^⊥= g. Consider now the function ω(k) = 1/m(k) for k /∈ Z, ω(k) = 0 for k ∈ Z. Then ωF [g] = F [f^⊥], and by the Fourier inversion theorem

f^⊥= F⁻¹(ωF [g]) =X

k /∈Z

1 m(k)gˆkφk.

Assume now h ∈ (ker T_m)^⊥ also follows T_mh = g. It follows that f^⊥− h ∈ (ker Tm)^⊥, but also T_m(f^⊥− h) = 0, so f^⊥= h, proving f^⊥ is unique in (ker T_m)^⊥.

3. Since {φk}k /∈Z is a basis of L²(T), T^m(Dm) ⊆ hφkik /∈Z is a direct consequence of ker Tm= hφkik∈Z. Furthermore, for k /∈ Z, Tm(φk/m(k)) = φk.

Note that we cannot say in general that Tm(Dm) = hφki_{k /}_∈Z, since Tm might not be defined in the whole L²(T) space. However, this is still a good rule of thumb when reasoning about such operators.

Introducing the notation R(T ) for the range of an arbitrary operator T , let us explore some important examples.

Examples 1.2.

1. ker ∂x= hφ0i. The inverse operator is

∂_x⁻¹g =X

k6=0

−i

kˆgkφk, g ∈ R(∂x),

which corresponds to the zero-mean integral of f . Unlike ∂_x, ∂_x⁻¹has bounded norm: k∂_x⁻¹φ_kk ≤ 1 for k 6= 0, so k∂_x⁻¹gk ≤ kgk. It can thus be extended by continuity to a bounded linear operator

∂_x⁻¹: hφ0i^⊥→ hφ0i^⊥ with unit norm. Further setting ∂_x⁻¹φ0= 0, it can be taken as bounded linear operator L²(T) → L²(T) with unit norm.

2. ker log |∂_x| = hφ₋₁, φ₀, φ₊₁i. The inverse operator is log |∂x|⁻¹g = X

|k|≥2

(log |k|)⁻¹ˆgkφk, g ∈ R(log |∂x|).

Again, log |∂_x|⁻¹ has bounded norm, so it can be extended by continuity to a bounded linear operator in hφ₋₁, φ₀, φ₊₁i^⊥ with norm (log 2)⁻¹ ' 1.44. Again, it can be further extended to L²(T) → L²(T) while conserving its norm by setting log |∂x|⁻¹φ_k = 0, |k| ≤ 1.

1.3 The Surface Quasi-Geostrophic Equation and a front model

Back to the physical problem at hand, the fractional Laplacion notation gives us a specific meaning of what (−∆)⁻¹ means in Equation (1.2). Furthermore, it lets us rewrite Equation (1.3) to reveal its common form in mathematics:

(11)

Definition 1.4 (Surface Quasi-Geostrophic Equation). The following set of equations is known as the Surface Quasi-Geostrophic Equation (SQG), where θ : R²→ R and u : R²→ R²

θ_t+ u · ∇θ = 0, u = ∇^⊥(−∆)^−1/2θ.

(1.5a) (1.5b) This equation is the focus of [12], where the formation of discontinuous fronts is studied. When written in mathematical terms, Equations (1.2) and (1.5) become particular cases of a one-parameter family of pseudo-differential equations:

Definition 1.5 (Generalized Surface Quasi-Geostrophic equation). The following set of equations is known as the generalized Surface Quasi-Geostrophic equation (gSQG), where θ : R² → R and u : R²→ R²,

θt+ u · ∇θ = 0, u = ∇^⊥(−∆)^−α/2θ,

0 < α ≤ 2.

(1.6a) (1.6b) (1.6c) The case α = 1 corresponds to the SQG Equation and the case α = 2 corresponds to the Euler Equations. An often-studied property of this family of equations is the existence of weak solutions that only take two possible values (see, for example, [13]):

θ(x, t) =

θ₊ x ∈ Ω(t), θ₋ x /∈ Ω(t),

for θ₊, θ₋ ∈ R and Ω(t) a domain whose boundary evolves over time. If the boundary is regular enough, one can model ∂Ω(t) as a function that follows a differential equation derived from the original gSQG Equation. Using vocabulary from Physics, this function is often called a front. In a recent paper [5], John Hunter and Jingyang Shu study fronts which are periodic perturbations of ∂Ω(0) = {y = 0}. In particular, they assume ∂Ω(t) can be parameterized as ∂Ω(t) = (x, ϕ(x, t)), with ϕ(x + 2π, t) = ϕ(x, t).

Using a cubic approximation for small amplitudes, they reach the following dynamic equation for the front:

ϕt+1 2∂x

ϕ²Aϕ − ϕAϕ²+1 3Aϕ³

+ Lϕ = 0, (1.7)

where A and L are operators that depend on the value of α. In the SQG case (α = 1), L = −2 log |∂x| and A = ∂_x²log |∂x|. This brings us to the main matter of interest in this thesis: expanding (1.7) for α = 1 we obtain the equation

ϕt+1 2∂x

ϕ²log |∂x|ϕxx− ϕ log |∂x|(ϕ²)xx+1

3log |∂x|(ϕ³)xx

= 2 log |∂x|ϕx. (1.8)

Hunter and Shu proved in a later paper [6] that this equation is weakly well-posed, but many questions remain open. The goal of this thesis is trying to provide an answer for one of these questions:

Does Equation (1.8) allow for time-periodic solutions?

The interest lies in the fact that the SQG Equation is known to have time-periodic solutions, but it is unclear if this new model captures this property. The rest of this thesis explains our investigations into this subject.

(12)

Chapter 2

Numerical simulation

As a first step to study Equation (1.8), we replicated the numerical results in [5]. Once this proved successful, we explored the use of Newton’s method for Banach spaces to find periodic soliton solutions.

This method showed the existence of possible candidate solutions, and inspired the theoretical approach explored in the next chapter. Before discussing these facts, we need to describe the theory behind numerical simulation using the Fourier transform. The approach towards Fourier transforms in this thesis, as exposed in appendix A, is mostly based on [4]. Other computational Fourier concepts are covered from an engineer’s perspective in [7].

2.1 Fourier pseudo-spectral methods

2.1.1 Theory

For each time point t ∈ R, the functions we deal with are elements of L²(T), so their Fourier transforms are elements of L²(Z) (see appendix A for more details) and they can be described using the Fourier basis {φk}_k∈Z. This is a very convenient setting for theoretical work, but a finite dimensional representation is needed to apply numerical simulations in a computer. The scheme followed by all the tools we used consists in sampling the original function at N equidistant points, so all operations can be taken in Z/N Z instead: given f (x) ∈ L²(T), we can associate it to a function f^N ∈ L²(Z/N Z) through the relation

fN(n) := f 2πn

N

, n ∈ Z/N Z.

We can think of this as the action of an operator:

Definition 2.1 (sampling operator). Let N ∈ N. The sampling operator S^N : L²(T) → L²(Z/N Z) is defined as

(SNf )(n) := fN(n) = f (2πn/N ), ∀f ∈ L²(T), ∀n ∈ Z/NZ.

This opens an important question: how should we represent the actions of the operators ∂x and log |∂x| over fN? The solution is applicable to all Fourier multipliers. Given a function f ∈ L²(T) and a symbol m : Z → C,

Tmf (x) = 1

√ 2π

X

k∈Z

m(k) ˆf (k)e^2πikx=X

k∈Z

m(k) ˆf (k)φk(x),

and we can introduce an equivalent basis for Z/N Z:

(13)

Definition 2.2 (Fourier basis of L²(Z/N Z)). The Fourier basis of L²(Z/N Z) is the collection {φ^NK}_{K∈Z/N Z}⊂ L²(Z/N Z) defined by

φ^N_K(n) = 1

√Ne^2πikn/N, ∀n ∈ Z/NZ, where k ∈ K is any representative of K.

Since e^2πikn/N is N -periodic in k, the choice of representative k ∈ K does not matter here. Further- more, the Fourier inversion theorem in Z/N Z implies that {φ^NK}_{K∈Z/N Z}is indeed a basis for L²(Z/N Z).

Any one choice of representative k ∈ K allows us to apply any L²(T) Fourier multiplier to fN through the formula

(T_mf_N)(n) = X

K∈Z/N Z

m(k) ˆf_N(K)φ^N_K(n), k ∈ K,

but there is an infinite number of choices of representatives and each one yields a different action of Tm over fN. However, there is a common property to the operators we are interested in: kTmφkkL²(T)

is an increasing function of |k|. This means that the L²-norm minimizing choice of representatives is {−N/2, . . . , +N/2} if N is odd, and either {−N/2 + 1, . . . , +N/2} or {−N/2, . . . , +N/2 − 1} if N is even.

This leads to the following definition:

Definition 2.3 (Action of Fourier multipliers over sampled functions). Let N ∈ N. The norm-minimizing choice of representatives K_N : Z/N Z → Z is defined as

KN(K) =

K ∩ {−N/2, . . . , +N/2}, N odd, K ∩ {−N/2, . . . , +N/2 − 1}, N even.

Let m : Z → C, and let T^mbe a multiplier with symbol m. The action Tm: L²(Z/N Z) → L²(Z/N Z) is defined as

(T_mf )(n) = X

K∈Z/N Z

m(K_N(K)) ˆf (K)φ^N_K(n), ∀f ∈ L²(Z/N Z), ∀n ∈ Z/NZ.

2.1.2 Practice

In practical terms, this means most of the computer’s work can be done representing any solution ϕ(x, t) of Equation (1.8) as a collection of timesteps ϕ^t(x) := ϕ(x, t), with each timestep containing N complex samples, or equivalently 2N floating point values, representing the frequencies ˆϕ^t_N. If these complex values are stored in an array indexed k = 0, . . . , N − 1, ∂x can be calculated as

F [∂xϕ^t_N](k) = i

(k mod N ) −N 2

ˆ ϕ^t_N(k),

and similarly for other relevant operators, like log |∂x| and ∂_x². The time evolution can then be calculated through a forward integration method like those in the Runge-Kutta family.

As a general principle, the precision of the simulation improves as N grows, and this particular problem seems to require especially high values – as a reference, Hunter and Shu use N = 2¹⁴and N = 2¹⁵in [5]. This makes it vital to keep the time complexity of any involved algorithms as low as possible. Since we represent the function as its Fourier transform, linear combinations and the application of Fourier multipliers are O(N ) operations. However, function products in the base space become convolutions in the dual space, which have order O(N²). Similarly, a naive implementation of the Fourier transform or its inverse is again an O(N²) operation. These two factors become, a priori, the worse bottlenecks in our simulation.

(14)

Luckily, for highly divisible values of N there are strategies that can be used to reduce the transform and inverse transform times to O(N log N ), collectively known as the Fast Fourier Transform methods (FFT). This also allows base space products (or rather, dual space convolutions) to be realised in O(N log N ) time: first, apply the inverse Fourier transform using FFT, second, perform the function product in base space as an O(N ) operation, third, apply the Fourier transform using FFT again.

Another important computational problem is that of numerical stability. Although the exact spectral form of Equation (1.8) is complicated, a general notion is that the k-th frequency is multiplied at every step by k³log k, ignoring convolutions (every term is differentiated three times and gets log |∂x| applied once). Furthermore, the convolutions expand the quickly growing error from the higher frequencies to the lower frequencies.This means that the evolution equations are stiff, and any numerical errors are expected to grow exponentially as the number of iterations increases. This can be attenuated by choosing smooth initial conditions, since the Fourier transform of a C^∞(T) function is rapidly decaying (see appendix A).

However, the choice of initial conditions is independent of the number of iterations, so it does not scale well with increasing values of N . A filtering step is needed at every step of the iteration process.

Filtering consists of convoluting the original function with a (base space) filter function, or equivalently multiplying the Fourier transform by a (dual space) filter function, so as to reduce unwanted artifacts in the simulation. Two common filters are the sharp spectral filter and the exponential filter :

Definition 2.4 (sharp spectral filter). Let M ∈ N. The sharp spectral filter with cutoff M is the Fourier multiplier with symbol

m(k) =

1, |k| ≤ M, 0, |k| > M.

Definition 2.5 (exponential filter). Let α, β > 0. The exponential filter with parameters α, β is the Fourier multiplier with symbol

m(k) = exp −α|k|^β .

These filters act on L²(Z/N Z) as described in Definition 2.3. A general rule of thumb is to discard the higher third of frequencies. For the sharp spectral filter, this means choosing M = N/3. For the exponential filter, [14] suggests using

m(k) = exp −36 2|k|

N

19! .

Another reason why filtering is needed is avoiding aliasing phenomena: the choice of representatives KN effectively embeds L²(Z/N Z) into L²(T) as the trigonometric polynomials of degree at most N/2.

However, the product in L²(T) of two degree-N/2 trigonometric polynomials is a degree-N polynomial, whereas the product in L²(Z/N Z) of two functions stays in L²(Z/N Z) and gets embedded as an at most degree-(N/2) polynomial. Hence, the high frequency information of an L²(T) product becomes mixed with the lower frequencies when this product is approximated in L²(Z/N Z). There are various ad-hoc ways to attenuate this problem; the theoretically simplest is to use a sharp spectral filter with cutoff N/2 to avoid any frequency overlap.

2.2 Replicating the results in Hunter-Shu

In [5], Hunter and Shu discuss the numerical simulation of Equation (1.8) with two sets of initial conditions:

(15)

Figure 2.1: Our own simulation with N = 2¹⁰for the initial condition (2.1).

ϕ₀(x) = cos(x + π) +1

2cos[2(x + π + 2π²)], (2.1)

ϕ0(x) = sech² 5(x − π) 2

. (2.2)

When we set out to replicate the numerical simulations in [5], our first approach consisted in programming the simulation directly in C. For this purpose, we used a commonplace FFT library, called FFTW [8], that implements O(N log N ) Fourier transform and inverse Fourier transform algorithms for highly divisible N , as discussed previously. Following Hunter and Shu, we restricted ourselves to powers of two to obtain high performance. The time iteration was realised using the standard 4th-order Runge-Kutta method (RK4) with a fixed time step. This involves rearranging Equation (1.8) as

ˆ

ϕt= F [T ϕ], T ϕ := 2 log |∂x|ϕx−1 2∂x

ϕ²log |∂x|ϕxx− ϕ log |∂x|(ϕ²)xx+1

3log |∂x|(ϕ³)xx

, (2.3) and calculating F [T ϕ] for various values of ˆϕ. In the customary Runge-Kutta notation of solving ˙y = f (y), this corresponds to y = ˆϕ, f (y) = F [T ϕ].

We tested the initial conditions (2.1) and (2.2), but also a simple sine f (x) = sin(x), the discontinuous sawtooth function f (x) = x, and low-degree trigonometric polynomials. Similarly, the code was tested not only against Equation (2.3), but also against two commonplace equations: the Transport Equation ϕ_t= ϕ_xand the Heat Equation ϕ_t= ϕ_xx. The simpler cases highlighted how unstable spectral numerical simulations are, and the need for filtering. As a consequence, we tried three options: no filtering, applying an M = N/3 sharp spectral filter at the end of the F [T ϕ] calculations, and applying the previously mentioned exponential filter at the end of the F [T ϕ] calculations.

(16)

Figure 2.2: Our own simulation with N = 2¹⁰for the initial condition (2.2).

Using RK4 with filtering and smooth initial conditions proved stable enough to correctly simulate the Transport and Heat Equations, but insufficient for Equation (2.3). This led to a second approach: using MatLab and its built-in routines. The achievable performance in MatLab is necessarily lower than pure C programming, but it allows for easier and faster iteration of the code. It also includes adaptive time step integration methods whose manual implementation would be beyond the scope of this thesis, and easy to use FFT methods (that internally use the FFTW library, according to MatLab’s documentation). We implemented integration through Matlab’s ode45 method, which uses the Runge–Kutta–Fehlberg method for adaptive timestep integration (see [9]).

The new approach was met with mixed success: the adaptive time step integration allowed for more stable results in the Heat and Transport Equations with lower values of N , but was still insufficient for Equation (2.3). This led us to contact Dr. Hunter, who suggested two modifications: using a viscosity approximation method, and applying an M = N/2 sharp cutoff filter before every single function product to avoid aliasing artifacts. It was only after these modifications that satisfactory results were achieved.

Figures 2.1 and 2.2 show the results we obtained for the initial conditions (2.1) and (2.2), respectively. In them, oscillatory phenomena similar to those described in the original paper can be observed, although considering the high numerical instability we cannot determine if this is a real property of the underlying equation or a property of the simulation.

2.3 A Newton method

After we succeeded in replicating the numerical results in [5], we moved on to explore numerically the possibility of periodic solutions in (1.8). This process can be simplified by considering soliton solutions, that is, constant shape solutions. This allows us to reduce the problem to functions f ∈ L²(T), represent them by their Fourier transforms, and apply Newton’s method.

(17)

2.3.1 Soliton solutions

For the context of this thesis, we define soliton solutions as follows:

Definition 2.6 (soliton). A map ϕ : T × R → R is called a soliton if there exist f : T → R, a ∈ R such that

ϕ(x, t) = f (x + at).

Assuming the solution is a soliton, Equation (1.8) can be written as af⁰+1

2∂x

f²log |∂x|f⁰⁰− f log |∂x|(f²)⁰⁰+1

3log |∂x|(f³)⁰⁰

= 2 log |∂x|f⁰, f : T → R, a ∈ R.

Similarly to Equation (2.3), this can be written in operator form:

Faf = 0, Faf := af⁰− 2 log |∂x|f⁰+1

2∂x

f²log |∂x|f⁰⁰− f log |∂x|(f²)⁰⁰+1

3log |∂x|(f³)⁰⁰

.

(2.4a) (2.4b)

Definition 2.7 (Hunter-Shu soliton operator). Let A ⊆ L²(T) be the set of functions f ∈ L²(T) such that

∃g ∈ L²(T), g = 1 2∂_x

f²log |∂_x|f⁰⁰− f log |∂x|(f²)⁰⁰+1

3log |∂_x|(f³)⁰⁰

,

where all terms in the equation must be valid expressions (in particular, f must be in the domain of log |∂x|∂x). The operator F : R × A → L²(T) is defined as F (a, f ) := F^af , as defined in Equation (2.4b).

Lemma 2.1 (Properties of the Hunter-Shu soliton operator).

1. Constant functions are a solution to Equation (2.4).

2. For each a ∈ R, the operator F^a:= F (a, ·) is unbounded.

3. The domain of log |∂_x| includes C¹(T), and k log |∂x|f k ≤ k∂_xf k for all f ∈ C¹(T).

4. C^∞(T) ⊆ A.

5. A is dense in L²(T).

Proof.

1. This is a simple check on Equation (2.4b).

2. Consider the sequence kFaφkk_L2, for k ∈ N:

F_aφ_k = ik(a − 2 log |k|)φ_k+1 2∂_x

φ²_klog |∂_x|(φk)⁰⁰− φklog |∂_x|(φ²_k)⁰⁰+1

3log |∂_x|(φ³_k)⁰⁰

=

= ik(a − 2 log |k|)φ_k+1

2∂_x 4 log 4 − 3 log 3 2π k²φ_3k

=

= ik(a − 2 log |k|)φk+ i3(4 log 4 − 3 log 3) 4π k³φ3k,

and since φk and φ3k are orthogonal with unit norm, this means kFaφkk ≥ |k|³C, where C :=

3(4 log 4 − 3 log 3)/4π.

(18)

3. Given n ∈ N, log n < n, so

k log |∂x|φkk = log |k| < |k| = k∂xφkk, |k| ≥ 1,

and log |∂x|φ0 = ∂xφ0 = 0. Therefore, k log |∂x|f k ≤ k∂xf k whenever ∂xf ∈ L²(T). Furthermore, C¹ functions are in the local integrable sets L^p_loc for all 1 ≤ p ≤ ∞, and T has finite measure, so C¹(T) ⊂ L^p(T). In particular, C¹(T) ⊂ L²(T).

4. This is a direct consequence of the previous point.

5. This is a direct consequence of the previous point and the fact that C^∞(T) is dense in L²(T).

F_af can further be subdivided into a linear part L_af⁰ and a nonlinear part N f : F_af = L_af⁰+ N f,

Laf := (a − 2 log |∂x|)f, N f := 1

2∂x

f²log |∂x|f⁰⁰− f log |∂x|(f²)⁰⁰+1

3log |∂x|(f³)⁰⁰

. Lemma 2.2 (Properties of La and N ).

1. For each a ∈ R, La is an unbounded linear operator.

2. Given n ∈ N, ker L2 log n= hφ_−n, φ_+ni. For other values of a ∈ R, ker La= {0}.

3. N is cubic, in the sense that N (λf ) = λ³N f . 4. N is unbounded.

5. Both La and N are defined in the dense subset A ⊂ L²(T).

Proof.

1. This follows from the fact that log |∂x| is an unbounded linear operator.

2. This can be obtained by direct calculation:

Laf = 0 ⇔ af = 2 log |∂x|f ⇔ a ˆfk = 2 log |k| ˆfk ∀|k| ≥ 1,

so that either ˆfk= 0 or a = 2 log |k|. In the case k = 0 this means ˆf0= 0, independently of a ∈ R.

3. This follows from the definition of N .

4. The same proof used for the unboundedness of Fa applies in this case.

5. The definition of A implies both these facts.

(19)

2.3.2 Newton’s method: theory

Newton’s method was originally designed for univariate functions and is usually taught for multivariate functions, but it can be generalized to Banach spaces. However, this requires a notion of differentiation for functions acting on those spaces – in our case, in the space of operators from L²(T) to itself. It turns out the notion of a derivative in Rⁿ extends naturally to this domain:

Definition 2.8 (Fr´echet derivative). Let A and B be Banach spaces. Let U ⊆ A be an open set, and let F : U → B. F is said to be Fr´echet differentiable at x ∈ U if there exists a bounded linear operator DF (x) : A → B such that

lim

h→0

kF (x + h) − F (x) − DF (x)hk_B khkA

= 0. (2.5)

DF (x) is called the Fréchet derivative of F at x. If F is Fréchet differentiable at every point in U , the map DF : U → L(A, B) mapping x to DF (x) is called the Fréchet derivative of F .

Lemma 2.3 (Properties of the Fr´echet derivative).

1. If it exists, DF (x) is unique. This extends to DF . 2. The map F 7→ DF is linear.

3. In the case A = Rⁿ, B = R^m, this corresponds with the usual notion of a multivariate derivative.

We shall skip the proof, since it is exactly the same covered in any multivariate course. The important point to keep in mind when working with spaces of infinite dimension is that DF (x) must be bounded.

As described in [10], this allows us to define Newton’s method in such spaces:

Definition 2.9 (Newton’s method). Let A, B be Banach spaces, U ⊆ A an open set and F : U → B a Fr´echet differentiable map with derivative DF . Newton’s method consists in calculating the sequence {Xn}^∞_n=0⊆ U , for arbitrary X0∈ U , and Xn+1implicitly defined by

DF (X_n)(X_n+1− X_n) = −F (X_n). (2.6)

If DF (x) is invertible for each x ∈ U , this can be rewritten in explicit form as

Xn+1= Xn− DF (Xn)⁻¹F (Xn). (2.7)

If DF (x) is not generally invertible, a strategy must be provided to find a valid Xn+1 at each step. A possible strategy is finding the minimum-norm least-squares solution, that is, finding the set of approximate solutions that minimize kDF (X_n)(X_n+1− Xn) + F (X_n)k² and then choosing the minimum-norm solution among those.

In general, Newton’s method is more effective if the initial value X₀ is close to a root. Since we are interested in nontrivial solutions, we may consider bifurcations from a known constant solution.

Considering the use of the operators ∂x and log |∂x| in Fa, we cannot expect DF (x) to be invertible in general. However, as pointed out in section 1.2.2, we might have a better chance inverting those operators if we restrict our solution to the space of zero-mean functions. This strongly suggests investigating bifurcations from f0= 0. Let us start from this point and try to calculate the Fr´echet derivative of Fa. It is worthwhile to introduce the notation

[F, G] := F G − GF, T f := log |∂_x|f⁰⁰.

We follow the same scheme we would use in multivariate calculus: given f0, h ∈ L²(T), we calculate the increment ∆Fa(f0)h = Fa(f0+ h) − Faf0, and then separate its linear and nonlinear parts. In this case,

∆F_a(f₀)h = ah⁰− 2 log |∂_x|h⁰+1

∂_xf₀²T h + 2f₀[h, T f₀] − [h, T f₀²] + O(h²).

(20)

Thus, the linear part is simply

DF_a(f₀)h = ah⁰− 2 log |∂x|h⁰+1

2∂_xf₀²T h + 2f₀[h, T f₀] − [h, T f₀²] , (2.8) and Equation (2.5) is mantained. However, this operator is not bounded for arbitrary f₀, independently of the value of a ∈ R. Consider the case f0= 0:

DFa(0)h = ah⁰− 2 log |∂x|h⁰= Lah⁰,

which we know to be unbounded. This means Fa is not Fr´echet differentiable in f0 = 0. However, the derivative is invertible in hφ₋₁, φ0, φ+1i^⊥, and its inverse is bounded. Therefore, Equation (2.7) is well-defined, and we can explore the results of Newton’s method.

Of course, we cannot choose X0 = 0, since then {Xn}^∞_n=0 would stay fixed in the trivial solution.

By the implicit function theorem, bifurcations can only happen for values of a ∈ R such that DF^a(0) is singular. Since we consider zero-mean functions, ∂x is invertible, and this means La must be singular.

As discussed in Lemma 2.2, this can only happen if a = 2 log n for some n ∈ N. Let us choose a fixed n ∈ N and consider the possible bifurcation from a0:= 2 log n as a one-parameter family of functions f_ε for ε > 0, where

F_a₀_+εf_ε= 0.

Consider then a series expansion of f_εby the parameter ε, that is,

f_ε= 0 + εf₁+ ε²f₂+ · · · =

∞

X

n=1

f_nεⁿ, f_n: T → R ∀n ∈ N.

This suggests a first-order approximation in terms of ε > 0 of the form fε ' εf1. Furthermore, since Fa₀+ε= Fa₀+ ε∂x,

∂ε(Fa₀+εfε) = Fa₀∂εfε+ ∂xfε+ ε∂ε∂xfε= Fa₀(f1+ 2εf2+ · · · ) + f_ε⁰+ ε(f₁⁰+ 2εf₂⁰ + · · · ), so

∂ε(Fa0+εfε)|ε=0= Fa0f1+ f₀⁰ = Fa0f1,

but since F_a₀_+εf_ε= 0, this means F_a₀f₁= 0. For f₁ real this implies f₁ = αφ_+n+ αφ_−n, α ∈ C. This suggests using the Newton method with the set of initial values

fˆ_k = 0, |k| 6= n, fˆ_+n= εα, fˆ_−n= εα, for a suitable choice of α ∈ C, ε > 0.

2.3.3 Newton’s method: practice

With all of this in mind, we wrote a script in MatLab that simulates the infinite-dimensional New- ton’s method in a finite-dimensional representation. As we did before, we represented f ∈ L²(T) as the L²(Z/N Z) element F[SNf ], with N a suitably large power of two. We calculated DF_a(X_n) as an N × N complex matrix following Equation (2.8), and solved Equation (2.6) using MatLab’s lsqminnorm least-squares norm minimizing algorithm.

There are several parameters we had to choose for this: the value n ∈ N that determines a0:= 2 log n and f1∈ hφ_−n, φ+ni, the value ε > 0, and the value α ∈ C. We expressed α in polar form as α = Ae^iθ.

(21)

Figure 2.3: l^∞ norm of successive iterations of ˆf , according to Newton’s method. Each line corresponds to a different set of initial values. Horizontal axis: number of iterations. Vertical axis: log k ˆf k_l^∞.

We set the script to repeat the simulation for each possible combination of values, with each parameter chosen from the following possibilities:

n ∈ {2, 3, 4, 5, 8, 11},

ε ∈ {10⁻³, 10⁻², 10⁻¹, 10⁰, 10¹, 10²}, A ∈ {10⁻¹, 10⁰, 10¹},

θ ∈ {0, π/2, π, 3π/2}.

The script runs a maximum of 10³iterations for each combination of values and uses the maximum absolute value of the stored frequencies (k ˆfNk_L^∞_{(Z/N Z)}, in representation of k ˆf kl^∞) to discard sequences that are considered divergent (the calculated value surpasses 10²⁰) or zero (the calculated value falls under 10⁻⁷). It records as successes the results of any combination for which the maximum absolute value stored for Faf , kF [Faf ]k_L^∞_{(Z/N Z)}, becomes smaller than 10⁻¹² before any other end conditions are met.

For the first version of this script, the results obtained were negative. However, the theoretical tools described in this chapter suggested a theoretical approach that is discussed in the next chapter. This theoretical approach has been so far inconclusive, but does shed clarity on the subject.

After later revision of the code, some implementation problems were found. This prompted a revision of the script, and yielded the final version described in this section. Figure 2.3 shows the maximum norm of X_n for different sets of initial values. As can be seen, most combinations either quickly converge to the trivial solution X_n = 0 or diverge, but a few solutions stay stable. Upon further investigation, some of these functions seem meaningless, but others suggest real solutions – and are remarkably similar in structure to what is described in the next chapter. Figures 2.4, 2.6 and 2.9 showcase some cases that were detected as successful, both promising ones and false positives.

(22)

-150 -100 -50 0 50 100 150 k

-0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04

F[f](k)

Figure 2.4: Real (blue) and imaginary (orange) parts of ˆfkfor the Newton method run with values n = 5, ε = 10⁻¹, A = 10, θ = 0, which was detected as successful. This is a purely real-valued successful case, showing symmetry with respect to k.

0 1 2 3 4 5 6 7

x -8

-6 -4 -2 0 2 4 6 8

f(x)

10^-5

Figure 2.5: Spatial form f (x) for the Newton method run with values n = 5, ε = 10⁻¹, A = 10, θ = 0.

(23)

-150 -100 -50 0 50 100 150 k

-0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04

F[f](k)

Figure 2.6: Real (blue) and imaginary (orange) parts of ˆfkfor the Newton method run with values n = 5, ε = 10⁻¹, A = 10, θ = π/2, which was detected as successful. This is a purely imaginary-valued successful case, showing anti-symmetry with respect to k.

0 1 2 3 4 5 6 7

x -8

-6 -4 -2 0 2 4 6 8

f(x)

10^-5

Figure 2.7: Spatial form f (x) for the Newton method run with values n = 5, ε = 10⁻¹, A = 10, θ = π/2.

(24)

-150 -100 -50 0 50 100 150 k

-0.005 0 0.005

0.01 0.015 0.02 0.025

F[f](k)

Figure 2.8: Real (blue) and imaginary (orange) parts of ˆfk for the Newton method run with values n = 11, ε = 10⁰, A = 10, θ = 0, which was detected as successful. This is a purely real-valued dubious case. Notice the scale of the peaks is less than 1/40.

0 1 2 3 4 5 6 7

x -2.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

f(x)

10^-5

Figure 2.9: Spatial form f (x) for the Newton method run with values n = 11, ε = 10⁰, A = 10, θ = 0.

(25)

Chapter 3

A theoretical approach

This chapter contains a theoretical approach to the existence of periodic solutions for the Soliton Equation (2.4). It is heavily inspired by section 2.3 and reuses all the theoretical background that we have exposed so far. Section 3.1 describes how a bifurcation fεfrom the trivial solution f0= 0, parameterized by ε > 0, can be expanded into a formal series. The terms in this series are defined recursively and provide a formal solution. Section 3.2 explores the convergence of this series using a majorant method for a simplified version of the equation (3.2.1) and for the full case (3.2.2); it is shown that this approach fails.

3.1 An assimptotic expansion

Following the steps of section 2.3, we want to find a one-parameter family of solutions fε, with ε > 0 small enough. Let us start from the most evident case, f0 = 0 and a0 = 2 log 2, and this time consider a_εas a function of ε too. Thus, we are looking for a mapping

ε 7→ (a_ε, f_ε) ∈ R × L²(T), (a₀, f₀) = (2 log 2, 0), F_a_εf_ε= 0.

Assuming a_εand f_εare regular enough, they can be simultaneously expanded as series on ε:

f_ε = 0 + εf₁ + ε²f₂ + · · · = P∞

n=1f_nεⁿ, f_n: T → R, aε = a0 + εa1 + ε²a2 + · · · = P∞

n=0anεⁿ, an∈ R. (3.1) We are interested in finding f_n that are regular enough and have zero mean, so the operators ∂_x, log |∂_x| and their combinations are well defined and are invertible (this would imply f_n ∈ L²(T) and f_n ∈ hφ₋₁, φ₀, φ₊₁i^⊥). In this scenario, Equation (3.1) lets us expand the linear part of F_a_ε, L_a_εf_ε⁰, as a series on ε, {a}^∞_n=0and {f_n⁰}^∞_n=1:

La_εf_ε⁰ = (aε− log |∂x|)∂xfε=

∞

X

n=0

anεⁿ− log |∂x|

!

∂x

∞

X

n=1

fnεⁿ

!

=

∞

X

n=1 n

X

m=1

an−mf_m⁰ − log |∂x|f_n⁰

! εⁿ=

∞

X

n=1 n−1

X

m=1

an−mf_m⁰ + La0f_n⁰

! εⁿ.

A similar calculation can be done on N f_ε, but it is more involved. For now, let us simply write

N f_ε=

∞

XεⁿN_nf_ε,

(26)

and note that, since N is cubic, N_nf_εcan be described by {f₁, . . . , f_n−2}. In particular, N₁f_ε= N₂f_ε= 0.

Putting it all together, we can write

Fa_εfε=

∞

X

n=1 n−1

X

m=1

an−mf_m⁰ + La₀f_n⁰ + Nnfε

! εⁿ, and since we expect Fa_εfε= 0 for all values of ε, this implies

0 =

n−1

X

m=1

a_n−mf_m⁰ + L_a₀f_n⁰ + N_nf_ε, ∀n ∈ N. (3.2)

This set of equations lets us give a formal recursive definition of an−1 and f_n⁰, assuming the values {a0, . . . , an−2} and {f1, . . . , fn−1} are known. To deduce these relations, we need to consider the kernel Z := ker La0 = hφ−2, φ+2i, and the projectors PZ : L²(T) → Z and PZ^⊥ : L²(T) → Z^⊥:

PZf = X

|k|=2

fˆkφk= ˆf−2φ−2+ ˆf+2φ+2,

PZ^⊥f = X

|k|6=2

fˆkφk.

Since ∂x, log |∂x|, PZ, PZ^⊥ and their linear combinations (including La0) are diagonal operators in the Fourier basis, they commute with each other whenever their joint action is well-defined. With this in mind, the cases n = 1, 2 provide initial conditions:

n = 1 ⇒ La0f₁⁰ = 0 ⇒ f₁⁰ ∈ Z,

n = 2 ⇒ La₀f₂⁰ = −a1f₁⁰ ⇒ PZ^⊥La₀f₂⁰ = 0.

Using commutation, La₀PZ^⊥f₂⁰ = 0. Separating f₂⁰,

L_a₀f₂⁰ = L_a₀(PZf₂⁰+ PZ^⊥f₂⁰) = 0 + L_a₀PZ^⊥f₂⁰ = 0,

so f₂⁰ ∈ Z and a1 = 0. A priori there are no restrictions on the precise values of f₁⁰ and f₂⁰ within Z, but we will treat both cases differently: let us choose f₁⁰ freely as any real-valued function, so that f₁⁰ = αφ₋₂+ αφ₊₂ for some nonzero α ∈ C^∗, and let us take f₂⁰ = 0. Assuming f_n to have zero mean, this uniquely determines f1and f2as

f1= ∂⁻¹_x f₁⁰ = i

2αφ−2− i

2αφ+2, f2= ∂_x⁻¹f₂⁰ = 0.

The fact that P^Zf_n⁰ has been undetermined so far turns out to be a recurring pattern which will allow us to obtain formulas for f_n⁰ and an. For this, consider the operator (La₀|_Z⊥)⁻¹ : R(La₀) → Z^⊥. This operator has norm (2[log 3 − log 2])⁻¹' 1.23, so we can consider its extension by continuity, the bounded linear operator G : Z^⊥ → Z^⊥. With this tool, we can provide the following result:

Lemma 3.1 (formulas for f_n⁰, an). Let α ∈ C^∗. There are a unique sequence {an}^∞_n=0⊂ C and a unique sequence {fn}^∞_n=1⊂ L²(T) ∩ C^∞(T) of zero mean functions that follow the Equation Set (3.2), and such that

a₀= 2 log 2, f₁⁰ = αφ₋₂+ αφ₊₂, f₂= 0, f_n∈ Z^⊥ ∀n ≥ 2.

Furthermore, each f_n⁰ is determined by

f_n⁰ = −G

n−1

X

m=2

a_n−mf_m⁰ + PZ^⊥N_nf_ε

!

, ∀n ≥ 3, (3.3)

A study of an approximate equation for sharp fronts in the generalized Surface Quasi-Geostrophic Equations in the cylinder

U.U.D.M. Project Report 2019:6

Examensarbete i matematik, 30 hp Handledare: Jordi-Lluís Figueras Examinator: Denis Gaidashev Februari 2019