Extended analysis of a pseudo-spectral approach to the vortex patch problem

(1)

UPTEC F18 024

Examensarbete 30 hp Juni 2018

Extended analysis of a pseudo-spectral approach to the vortex patch problem

Mattias Bertolino

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:

Box 536 751 21 Uppsala Telefon:

018 – 471 30 03 Telefax:

018 – 471 30 00 Hemsida:

http://www.teknat.uu.se/student

Abstract

Extended analysis of a pseudo-spectral approach to the vortex patch problem

Mattias Bertolino

A prestudy indicated superior accuracy and convergence properties of a pseudo-spectral method compared to a spline-based method implemented by Còrdoba et al. in 2005 when solving the α-patches problem. In this thesis we further investigate the numerical properties of the pseudo-spectral method and make it more robust by implementing the Nonequispaced Fast Fourier Transform. We present a more detailed overview and analysis of the pseudo-spectral method and the

α-patches problem in general and conclude that the pseudo-spectral method is superior in regards to accuracy in periodic settings.

ISSN: 1401-5757, UPTEC F18 024 Examinator: Tomas Nyberg Ämnesgranskare: Gunilla Kreiss Handledare: Jordi-Lluis Figueras

(3)

Acknowledgement

I wish to thank my supervisor Jordi-Lluis Figueras for his patient guidance, discussions and support during the thesis. I would also like to thank my subject reader, Professor Gunilla Kreiss, for her interest in what was achieved and her valuable feedback.

Special thanks to the authors of Nonequispaced fast fourier transform for the support pro- vided by email and Annette Ekblom for her excellence as a courier. I wish to thank Adam Waks for discussion regarding Fourier transforms and nectarine. All of your inputs have contributed considerably to this thesis.

Infine vorrei ringraziare Graziano Bertolino per il suo aiuto con i disegni, och resten av min familj för deras stöd och glädje.

Keywords

α-patches problem, discontinuous vorticity equations, 2D Euler equations, Surface Quasi-Geo- strophic equations, self-similar solutions, Pseudo-spectral methods, Nonequispaced fast fourier transform.

(4)

Populärvetenskaplig sammanfattning

En av matematikens olösta gåtor är huruvida Navier-Stokes ekvationer är välställda. Det är ett av Clay Mathematics Institutes Millenniumproblem som formulerades för att lösa några välkända och olösta problem. Att ekvationerna är välställda betyder att rimliga värden in i ekvationerna ger rimliga värden ut. Det är av stor vikt eftersom ekvationerna används för att modellera en mängd fenomen i naturen som kan beskrivas med fluider så som orkaner och vattenflöden men även luftflöden efter jetmotorer och kring flygplansvingar.

För att få insikt i hur Navier-Stokes ekvationerna beter sig kan det vara av nytta att först under- söka förenklade, men relaterade, problem. Ett sådant är Euler-ekvationerna som också beskriver fluider, men där viskositeten, ungefär fluidens inre friktion, har bortsetts ifrån. Vad gäller Euler- ekvationerna så är det bevisat att de i två dimensioner är välställda men liksom för Navier-Stokes ekvationer saknas ett sådant bevis för Euler-ekvationerna i tre dimensioner. Matematiker har dessutom funnit att ekvationer vid namn Surface Quasi-Geostrophic Equations i två dimensioner har stora matematiska likheter med Euler-ekvationerna i tre dimensioner. Det vore alltså intressant att undersöka hur en blandning av Euler-ekvationerna i två dimensioner och Surface Quasi Geostrophic ekvationerna beter sig.

Båda ekvationerna beskriver hur ett område med konstant rotation rör sig när det är om- slutet av ett område utan rotation, vilket kan liknas vid orkaner. Ett första steg är att använda beräkningsmetoder för att undersöka hur ett rotationsområde förändras i tiden och det här arbetet jämför i första hand olika beräkningsmetoder för att beskriva förändringen av rotationsområdet i tid. Viktigt vid utveckling av beräkningsmetoder är att de ska uppfylla kriterier på noggrannhet, robusthet och tidskomplexitet.

I en förstudie utvecklades en beräkningsmetod som både var mer noggrann och snabbare än tidigare utvecklade metoder. Däremot har det visat sig att den utvecklade metoden bara var noggrann för tydligt tillrättalagda problem och därför oanvändbar i praktiken. I den här studien söker vi lösning på problemet och presenterar en metod som på nästan lika lång tid beräknar utvecklin- gen av rotationsområden noggrant i allmänhet och med möjlighet att placera ut beräkningspunkter där det är mer intressant att studera; till exempel där området böjs kraftigt.

Förutom att bättre verktyg ökar möjligheterna att i detalj studera dessa rotationsområden så ökar möjligheterna att hitta nya infallsvinklar att studera Navier-Stokes ekvationer samt ökar förståelsen för hur fluider beter sig.

(5)

1 Introduction

Discontinuous vorticity problems have been studied by both mathematicians and physicists and with the increased development of numerical methods and scientific computing, the possibilities to study these problems from a different perspective have increased. One commonly studied discontinuous vorticity problem is the vortex patch problem in which there exists a constant non-zero local vorticity in a two-dimensional bounded region D(t). Vortex patches can be used as prototyp- ical models of hurricanes where the actual region exhibits a strong vorticity inside the hurricane but where the surroundings appear completely irrotational.

The main reason to study the vortex patch problem is to further gain knowledge whether or not the Navier-Stokes equations develop singularities in finite time. They are used to model flows of viscous, homogeneous, fluids in Rⁿ, n = 2 and 3, and appear in a range of physical phenomena since they couple how the velocity, pressure, temperature and density of a moving fluid are related.

Since the formalization of the Millennium Prize problems [1], the answer to possible formation of singularities in Navier-Stokes in three dimensions has further gained industrial and academic interest. The many components of the Navier-Stokes equations make it possible to, in different ways, relax some of the properties in order to deal with a simpler, but possibly still informative, problem. If, for example, the viscosity of the moving fluid is neglected, the Navier-Stokes equations reduce to the Euler equations.

The vortex patches are related to Navier-Stokes and the Euler equations since they are weak solutions to the 2D Euler Equations (2DEE) in the space L¹(R²)∩L^∞(R²) as proved by Yudovich [30]

and extended by Majda [23]. However, the proofs do not suffice as proofs of uniqueness. As will be clarified in the following sections, the relation between the stream function and vorticity in 2DEE can be altered in order to yield the Surface Quasi-Geostrophic Equation (SQGE) which has several similarities with the 3D Euler Equations (3DEE), despite being a 2D model. The similarities to 3DEE as well as the interpretation and computational advantages of SQGE being a two-dimensional model makes it a popular choice to study (see e.g. [9,10,27,17,8]).

The link between 3DEE and SQGE is the key in this thesis as well as in the above mentioned studies. If one would present a rigorous proof of blow-up (or the absence of it) in SQGE, it could possibly be extended to 3DEE and Navier-Stokes and perhaps one way to gain deeper knowledge of what would cause a blow-up is to study the α-patches problem. By interpolating the linear operator of 2DEE and SQGE by introducing a parameter α, it is possible to study a problem related to both 2DEE and SQGE. This is of interest since 2DEE has been rigorously proven to not develop singularities whereas SQGE has yet not.

In a prestudy to this thesis, [3], the numerical velocity field evaluation of vortex patches of the type α-patches was studied using a Fourier spectral method with collocation points, a pseudo-spectral method, as compared to a cubic interpolation-based method first introduced by Dritschel [12]. The pseudo-spectral method proved superior in terms of both accuracy and speed although the use of a pseudo-spectral method introduces two sources of error, that is truncation error and aliasing, which need to be taken care of. The main source of error comes from the placement of the nodes, which by Fast Fourier Transform (FFT) is assumed to be uniform with respect to a parametrization of the patch’s contour. If using an adaptive node-handling scheme to increase the node density at high-curvature segments, the number of frequencies to account for in FFT increases. Therefore, the comparisons performed in the prestudy are only applicable in the special case of uniform sampling of the curve, which is only obtained in the initial setup of the patches.

In this thesis, we extend the comparison between Dritschel’s cubic interpolation-method and pseudo-spectral methods for velocity field evaluation to include more general samplings of the curve, with varying node density along the curve segments. We are interested in investigating the pseudo-spectral method in evaluating the divergent integral with numerical quadratures. We do so by first comparing two submethods of the pseudo-spectral method in the frequency expansion step using FFT. The first submethod method uses a standard FFT approach as in the prestudy to retrieve samples of the integrand of the contour dynamics equation. The second method uses the Nonequispaced Fast Fourier Transform (NFFT), in which an oversampled FFT is applied to a deconvolution of the spatial signal with a window function. This is done in a pretest in order to explain the accuracy results when evaluating the whole contour integral.

In addition to the introduction of NFFT applied in the frequency expansion step of the pseudo- spectral method, we introduce a proof-of-concept of a radial formulation of the contour dynamics

(7)

equation. This formulation can be used with any pseudo-spectral method and relies on an extended polar formulation of the evolution of the patch by using Frenet-Serret formulas to describe the radial evolution of a patch from an initial patch. This directly solves the normal projected contour integrals and transform the 2D problem to a 1D problem, hence cutting the computational costs by a factor two. Due to time-constraints, this has not been implemented numerically and is merely left with the theoretical motivations.

Apart from studying the performance of the pseudo-spectral approach in the usual periodic contour dynamics equation setting, we theoretically study how adaptable pseudo-spectral methods are when integrating along a contour which is not periodic. Córdoba et al. [11] found that a scaling law exists near the formation of singularities in the α-patches setting which led them to propose a system of rescaled variables which admits a stationary profile when singularities are formed. The numerical part of the rescaled problem in this thesis simulates a perturbation of a fixed point to examine whether it is attracting nearby solutions or not. An attractive fixed point would be of significant importance in the search for rigorous proofs of singularities in the α-patches problem which possibly could be extended to the more involved 3DEE.

In the rescaled variables, contrary to the periodic setting, the area of the patches approaches infinity as the singularities are formed. This means that the fixed point we look for can be expressed as a function on R. The infinite domain poses a numerical problem, however, since we are in need of a finite set of nodes. We thus investigate the Clenshaw-Curtis quadrature (CC) for its simplicity and convergence properties on infinite domains which is based on Chebychev-polynomial expansion of the integrand and the Trapezoidal rule quadrature. Both the Riemann sum, used in the periodic problem and CC quadratures are weighted quadratures in which the integral is approximated as an inner product of the samples of the integrand and their corresponding weights. In the Riemann sum, the weights are identical 1/N , assuming equidistant spacing of the samples, and in the CC, the weights are given in terms of the Chebychev-polynomials.

There are mainly two reasons for investigating the CC quadrature. The first reason is because of the possibility to formulate the Chebychev-polynomial weights as a type-I discrete cosine transform implying that the weights can be efficiently retrieved with a computational complexity of O(N log N ). Secondly, it includes a transformation to an even, periodic integrand which is useful on aperiodic domains or over an infinite interval. This makes it possible for us to study the vicinity of the fixed point using numerical tools without simplifying the model by using large, but finite, areas of the patches.

The drawback of CC is that it can only be used on decaying integrands. This is a concern as it hinders the use of it on the perturbation of the fixed point proposed by Mancho [24]. This leads us to propose a conjecture that the rescaled patches can be decomposed into a decaying part solved pseudo-spectrally with CC and a non decaying part solved with another scheme. Although this conjecture has not yet been verified, we present pretests on solving known integrals of decaying functions to show the usability of CC.

Further, we accompany the numerical tests by formalizing the pseudo-spectral method for standard FFT as frequency expansion submethod and further investigate the numerical properties of a pseudo-spectral approach to evaluate the velocity field in the α-patches problem by bounding the truncation error and the aliasing using appropriate norms.

This thesis is partly a broad introduction to using pseudo-spectral methods for solving discontinuous vorticity problems. The goal is not only to present accuracy results compared to Dritschel’s method but also to lay the groundwork in terms of theoretical error estimations of the frequency expansions and to complement the accuracy results with novel ideas on how to use the pseudo- spectral method also on infinite domains.

Discontinuous vorticity problems have several facets to them and from a numerical perspective it needs to be treated with care. In a crude sense, their continuous formulations are reduced to a system of functional differential equations that can be solved by series expansion and a time- integrator of choice. However, the evaluation of the velocity field required in the evolution entails several difficulties including the numerical quadrature of divergent integrals, which in the rescaled problem is over infinite intervals, and high-curvature which requires a compromised global and local node handling scheme.

The outline of this thesis is as follows. In Section2 we give an introduction to the α-patches problem by showing analogies of SQGE and 3DEE both in regards to conservation of physical properties and conditions of the forming of singularities. We further provide more extensive theoretical description of the physical system studied in [3] and go deeper in the analogies between

(8)

the α-patches problem and 3DEE than was done in the prestudy. We also include a derivation of the α-patches problem’s contour dynamics form from its strong form formulation. The described α-patches problem in this section is identical to the one in the pre-study as it deals with area-preserving patches over a finite domain. We further present the α-patches problem in the extended polar form and in the rescaled, self-similar, variables to broaden the theoretical aspects of the problem.

The formalization of the pseudo-spectral method is given in Section 3. Here, we state the error bounds of the pseudo-spectral method when using bare FFT as the submethod for frequency expansion. We also give a theoretical description of the Nonequispaced Fast Fourier Transform and describe the notion of the Riemann sum used for periodic settings and Clenshaw-Curtis used for infinite domain settings.

In Section4 we give the implementation of NFFT in terms of used window functions, planner flags, etc. We further state how CC can be used as the chosen quadrature for a pseudo-spectral method on the self-similar problem, given that the conjecture that splitting the integrand in a decaying and non-decaying function holds.

In Section5we give the setup for the numerical measurements performed in this thesis, alongside the assumptions that motivates the measurements. We further present both the software and hardware used for simulations and as testing platforms.

Further, the obtained numerical results are presented in Section6 and discussed in Section7.

The conclusions drawn from the discussions and the thesis as a whole are given in Section8.

(9)

2 The vortex patch problem

In this section, we present the α-patches problem and its relation to other discontinuous vorticity problems through its analogies to 3DEE. The α-patches problem has essentially two components in the form of 2DEE and SQGE in which the α-patches problem can be seen as an interpolation between the linear operators relating the stream and vorticity of the two. Thus, in order to understand the importance of studying the α-patches problem, the analogies between 3DEE and SQGE need to be made clear. As an introduction to the α-patches problem, we therefore present SQGE and 3DEE in vorticity form and explicitly remarking both the similarities equation-wise, but also in terms of the conserved quantities. We further state two theorems with conditions for singularities to form for both SQGE and 3DEE. These three types of similarities are the core motivation for studying SQGE and hence the α-patches problem.

The background of SQGE is from geophysics and was introduced by Charney et al. [6] to model atmospheric turbulence. However, after finding that SQGE and 3DEE have several properties in common, Constantin, Majda and Tabak [8] suggested that SQGE could be used as a 2D model for 3DEE in which the potential temperature of SQGE plays the role of the vorticity in 3DEE.

The α-patches problem, which is numerically studied in this thesis, is stated in contour dynamics equation form and to reach that, we first present the concept of a vortex patch in the context of 2DEE. We further present the α-patches problem in extended polar form and in the rescaled, self- similar, form. For the polar form, we provide a proof that there exists local coordinates close to the reference curve used to represent the patch, which ensures us that the polar form is equivalent to the normal-projected velocity field in cartesian formulation.

2.1 3D Euler and Surface Quasi-Geostrophic

SQGE is given by

∂tθ(x, t) + u(x, t) · ∇θ(x, t) = 0, (x, t) ∈ R²× R+, (1) The vorticity formulation for the incompressible 3DEE, on the other hand, is given by

∂tω(x, t) + u(x, t) · ∇ω(x, t) = ∇u(x, t) · ω(x, t), (x, t) ∈ R³× R+. (2) In both SQGE and 3DEE, the gradient of the velocity u is given by singular integral operators of convolutional type (Riesz transforms in Euclidean spaces) of the vorticity θ and ω. For SQGE that is

u(x) = − Z

R²

∇^⊥θ(x + y)

|y| dy, (3)

and likewise for 3DEE,

u(x) = 1 4π

Z

R³

y × ω(x + y)

|y|³ dy.

By taking the normal derivative, given by ∇^⊥= (−∂x₂, ∂x₁), of Equation (1), we obtain

∂t(∇^⊥θ) + u · ∇(∇^⊥θ) = (∇u) · (∇^⊥θ),

similar to Equation (2). The quantity ∇^⊥θ is an active scalar, meaning that it modifies the physical properties like density, contrary to a passive scalar which does not. The analogies of Equation (2) and Equation (1) are clear by letting the active scalar ∇^⊥θ in SQGE play the role of vorticity in 3DEE. In both systems the velocity, as well as the vorticity ω and the active scalar ∇^⊥θ are divergence free in two and three dimensions respectively.

Further, it follows from Equation (1) that ∇^⊥θ is tangent to the level sets in SQGE, θ is constant, and that the level sets move with the fluid flow since it is a transport equation. The level sets are hence analogs of the vortex lines in 3DEE where the vorticity is tangent to the vortex lines.

Regarding the conserved quantities in both equations, they have several properties in common.

For example, the L²-energy is conserved in time for both 3DEE and SQGE, that is

||u||_L²(T )= ||u||_L2(0),

for any fixed time T . Further, Beale, Kato and Majda [2] showed a sufficient condition for a singularity to form at time T for 3DEE, presented in Theorem1.

(10)

Theorem 1 (BKM Necessary and sufficient singularity condition). Let u be a solution to 3DEE in Equation (2) in the class u ∈ C ([0, T ]; H^s) ∩ C¹ [0, T ]; H^s−1. Assume further that there exists a time T_∗ such that the solution is excluded from the class. Then the following holds.

Z T∗

0

||ω(t)||_L^∞dt → +∞.

Theorem1assures us that if we control the L^∞norm of the vorticity, then there is no singularity formed. Based on this result Constantin, Fefferman and Majda [7] found that if you control the direction field of the vorticity, then no singularities form; as presented in Theorem2.

Theorem 2 (CFM Singularity condition). If the direction field of the vorticity, ξ = ω

|ω|,

is smooth in and near the regions of high |ω|, then singularities cannot form.

Theorems 1- 2 also holds when interchanging ω ↔ ∇^⊥θ, namely by moving from 3DEE to SQGE.

2.2 Vortex patches

A vortex patch is a bounded, simply connected, region D(t) ⊂ R²of constant vorticity,

θ(x, t) =

(θ0, x ∈ D(t), ϑ(x, t), x ∈ R²\D(t).

In our case, we deal with the special case where the vorticity outside the patch is constant zero, that is ϑ(x, t) ≡ 0. Note that we use θ to denote vorticity, rather than the more common ω. In SQGE setting, which will be of more of interest to study than 2DEE in the following sections, θ is the well-adopted symbol. In both these equations, the vorticity and the stream function are related by a linear operator, L(ψ) = θ. If in the case of SQGE, the linear operator is given by L(ψ) = (−∆)^1/2ψ. If, on the other hand, in the case of 2DEE we have the relation between the stream function and the vorticity given by Laplace’s equation,

L(ψ) = ∆ψ. (4)

Vortex patches are weak solutions to the 2DEE and it is possible to invert the Laplace’s equation to obtain a contour dynamics equation. In the following section, we derive this form for the α-patches problem to be able to describe the evolution of a patch completely determined by the evolution of the boundary of the patch.

2.3 Contour dynamics equation

All the analysis in this thesis, both numerical and theoretical, deals with the equations transformed to the contour dynamics equation form. This form is an equivalent formulation of the equations which is useful since it transforms the problem to a system of functional differential equations which can be solved for the velocity of the curve and hence time integrated for the evolution of a curve. In essence, the transformation is done by inverting the linear operator that relates the stream and vorticity with the use of the associated Green’s function followed by the use of Green’s formula to move from a surface integral to a contour integral.

In the case of 2DEE, the linear operator relating the vorticity and the stream is given by Laplace’s equation in Equation (4), and in order to invert Laplace’s equation, we need Green’s function associated with Laplace’s equation. In the following numerical studies, we do not directly work with the Laplacian as in 2DEE but rather with a fractional Laplacian in the α-patches problem. The derivation of the contour dynamics equation for an α-patch is given in Section2.3.1.

However, the derivation of the contour dynamics equation for 2DEE serves as an aperitif to the more involved derivation of the case for the fractional Laplacian, and hence the contour dynamics equation for 2DEE is derived in AppendixA.

(11)

2.3.1 Contour dynamics equation of the α-patch

The α-patches problem can be seen as an interpolation between the two linear operators in 2DEE and SQGE by introducing a parameter α ∈ [0, 1] which yield the relation between the vorticity and the stream in the fractional Laplace’s equation

L(ψ) = (−∆)^1−α/2ψ = θ, (5)

where α = 0 is equivalent to 2DEE and α = 1 is equivalent to SQGE. As stated in Section2.3, we need to find the Green’s function associated with the linear operator, but instead of Green’s function associated with Laplace’s equation, we need to find Green’s function associated with the fractional Laplacian in Equation (5), that is to find G(x, x⁰) such that

(−∆)^1−α/2G2(x, x⁰) = δ(x − x⁰),

where δ(x − x⁰) denotes the Dirac delta function centered at x⁰. In two dimensions, this is given by

G2(x, x⁰) = c_2,1−α/2 1

|x − x⁰|^1+α, (6)

in accordance with Kwaśnicki [22]. The fraction-dependent constant cd,s in d dimensions and fraction s is given by

cd,s= Γ(^d−s₂ ) 2^sπ^d/2Γ(^s₂).

Similar to the 2DEE case in AppendixA, the stream function is retrieved with the Riesz transformation with the kernel as in Equation (6)

ψ = θ0

2π Z Z

R²

dx⁰

|x − x⁰|^1+α/2, and the velocity is given by

u(x, t) = θ0

2π Z Z

R²

∇^⊥ dx⁰

|x − x⁰|^1+α/2,

where we have absorbed the constant arising in Green’s function in the vorticity θ₀ = ˜θc_2,1−α/2. By the use of Green’s formula, we get the evolution of an α-patch expressed in terms of a contour integral as given by

∂x

∂t(γ, t) = θ0

2π Z

∂D(t)

∂x

∂γ(γ⁰, t)

|x(γ, t) − x(γ⁰, t)|^αdγ⁰. (7) The integrand in Equation (7) is singular at x(γ, ·) = x(γ⁰, ·), but its Cauchy Principal value exists. For computational reasons, the tangential component of Equation (7) is removed since only the normal component of the velocity field is able to deform the curve, in accordance with the divergence theorem. This way, the singularity in the integrand is removed without the loss of information. We end up with

∂x

∂t(γ, t) _n

= θ0

2π Z

∂D(t)

∂x

∂γ(γ⁰, t) −^∂x_∂γ(γ, t)

|x(γ, t) − x(γ⁰, t)|^αdγ⁰, (8) which is what has been used in the implementation.

To summarize the connections and similarities of the above discussed vorticity problems, and how the α-patches problem is related to them, we present three flowcharts in Figure1which mark out how the parameter α can be chosen between 0 and 1 to move between SQGE and 2DEE. The red color of the entities Navier-Stokes, 3DEE and SQGE indicates that absence of blow-up is yet to be proven, the blue color of 2DEE indicates that absence of blow-up has been proven and the purple color of the α-control indicates that the dynamics of the α-patches problem tend more to either SQGE or 2DEE depending on the chosen value.

(12)

α Navier-

Stokes

3D Euler

Surface Quasi- Geostrophic

Neglect viscosity

2D Euler

(a)

α Navier-

Stokes

3D Euler

Neglect viscosity

2D Euler

(b)

α Navier-

Stokes

3D Euler

Neglect viscosity

2D Euler

(c)

Figure 1: The relation of different vorticity problems and the tuning of the α-patches problem to inhibits the dynamics of (a) 2DEE (b) 2DEE and SQGE (c) SQGE.

2.4 Polar formulation of the contour dynamics equation

An alternative to the cartesian formulation of the contour dynamics equation is in extended polar coordinates. Assuming a patch can be expressed in polar form, the normal-projected evolution of the patch is equivalent to the evolution of the radial component according to the divergence theorem. Starting from Equation (8), describing the evolution of the normal component of the curve, we make the ansatz

x(γ, t) = R(γ, t)ψ(γ) = R(γ, t)e^iγ. (9)

In these coordinates, Equation (8) becomes

∂x

∂t(γ, t) _n

=∂R

∂t(γ, t)ψ(γ) = θ₀ 2π

Z

∂D(t)

∂R

∂γ(γ⁰, t)eîγ⁰ + iR(γ⁰, t)eîγ⁰−^∂R_∂γ(γ, t)eîγ− iR(γ, t)eîγ

|R(γ⁰, t)e^iγ⁰− R(γ, t)e^iγ|^α dγ⁰. The evolution of the radial component is thus given by

∂R

∂t(γ, t) = θ0

2π Z

∂D(t)

∂R

∂γ(γ⁰, t)e^i(γ⁰^−γ)+ iR(γ⁰, t)e^i(γ⁰^−γ)−^∂R_∂γ(γ, t) − iR(γ, t) (R(γ⁰, t)²+ R(γ, t)²− 2R(γ⁰, t)R(γ, t) cos(γ⁰− γ))^α/2 dγ⁰. By the use of Euler’s formula, the right hand side is simplified to

∂R

∂t(γ, t) = θ0

2π Z

∂D(t)

∂R

∂γ(γ⁰, t) cos(γ⁰− γ) − R(γ⁰, t) sin(γ⁰− γ) − ^∂R_∂γ(γ, t) (R(γ⁰, t)²+ R(γ, t)²− 2R(γ⁰, t)R(γ, t) cos(γ⁰− γ))^α/2dγ⁰, where the imaginary terms are equal to zero, since the patch is real. The computational reasons to formulate the problem in polar form is to only treat the radial component. In this formulation, the degrees of freedom is reduced from two to one, cutting the number of operations to a half. Further, having a curve described is polar coordinates directly yields a parametrization of the curve.

Using standard polar coordinates, as in Equation (9), is not a general ansatz however as not all curves can be expressed in polar form. We therefore propose the use of extended polar coordinates in accordance with Theorem3to treat simple (i.e. non-overlapping) curves.

(13)

Theorem 3 (Extended polar coordinates). Given a reference curve C(γ) such that C : [0, 1] → R²,

we can represent any nearby curve x(γ) as

x : [0, 1] → R², by writing

x(γ) = C(γ) + p(γ)n(γ),

where p(γ) is the scalar distance from the reference curve along its normal n(γ) to the curve x(γ).

The theorem essentially states that for any parametrizable curve, we can find a nearby curve represented with the same parametrization. To prove this, we prove that the curve has local coordinates for each bounded segment of the curve, namely that the characteristics of the differential equation do not overlap at any nearby point to the curve. In the following proof, we prove Theorem3by proving C : [−ε, ε] → R².

Proof. The proof is given by contradiction. Suppose that we have

C(˜γ) − C(γ) = pn(γ) − ˜pn(˜γ), (10)

for ˜γ 6= γ, ˜p 6= p. Locally, ˜γ = γ + ε, and ˜p = p + δ. By Taylor expansion around γ we get C(γ + ε) − C(γ) = C⁰(γ)ε + O(ε²)

n(γ + ε) = n(γ) + n⁰(γ)ε + O(ε²).

By inserting this in Equation (10), we obtain

C⁰(γ)ε = n(γ)δ + O(ε², εδ). (11)

Since C⁰(γ) and n(γ) are orthogonal, they do not cancel out. Likewise, the higher order terms do not cancel out the first order terms. Hence Equation (11) does not hold for ε, δ > 0, which concludes the proof.

In Figure 2 we illustrate the idea of local coordinates in each segment of the reference curve.

Locally, the curvature is negligible.

C(γ)

Figure 2: Segment of a reference curve C(γ) with local coordinates.

With the extended polar coordinates as in Theorem3 we can write Equation (7) in terms of these. We have

∂x

∂t(γ, t) = θ 2π

Z

∂D(t)

∂C

∂γ(γ⁰, t) +^∂p_∂γ(γ⁰, t)n(γ⁰, t) + p(γ⁰, t)^∂n_∂γ(γ⁰, t)

|C(γ, t) + p(γ, t)n(γ, t) − C(γ⁰, t) − p(γ⁰, t)n(γ⁰, t)|^αdγ⁰, which yields the evolution of distance p(γ, t) by following the same steps as above.

To illustrate the representation of a curve using extended polar coordinates we present two contours described by two different reference curves, namely a circle and an arbitrary curve in Figure3(b).

(14)

C(0) C(γ

¹

)

x(γ

¹

,t) p(γ

¹

,t)

(a)

C(γ

¹

)

C(0) x(γ

¹

,t)

p(γ

¹

,t)

(b)

Figure 3: Illustrations of extended polar coordinates using arbitrary curves represented by a reference curve.

2.5 Self-similar problem in rescaled variables

While the previous section highlighted the importance of studying α-patches by remarking its analogies to 3DEE and derived its contour dynamics equation form, this section poses the numerically studied problem in rescaled self-similar space variables. If a singularity is formed at time t_∗ in the coordinates x_∗(t_∗), we are able to rescale the space variable in Equation (8) in accordance with Mancho [24]. The rescaled variable y(η, t) is given by

y(η, t) − y∗= x(η, t) − x_∗(t) (t_∗− t)^δ ,

where the parameter δ is restricted to be equal to 1/α for physical solutions, which follows from Theorem2. At the collapse point, the velocity field is given by

dx_∗(t) dt =

N_k

X

k=1

θk

2π Z

∂D(t)

∂x_k

∂η (η⁰, t)

|x_∗(t) − x_k(η⁰, t)|^αdη⁰.

(15)

In the rescaled variables, the evolution of the system is thus given by

(t_∗− t)∂x

∂t − δ(y − y_∗) =

N_k

X

k=1

θk

2π Z

Υk(t)

∂y_k

∂η (η⁰, t)

|y(η, t) − y_k(η⁰, t)|^α −

∂y_k

∂η (η⁰, t)

|y_∗− y_k(η⁰, t)|^α

!

dη⁰, (12) where Υ_k represents the kth rescaled contour. By introducing the pseudo time

τ = − log(t_∗− t), Equation (12) projected on the normal of the curve is given by

∂y

∂τ(η, τ ) _n

= δ(y − y∗) _n + θ₀

2π Z

Υ₁(τ )

∂y₁

∂η(η⁰, τ ) −^∂y_∂η(η, τ )

|y(η, τ ) − y(η⁰, τ )|^αdη⁰ (13a)

− θ0

2π Z

Υ₂(τ )

∂y₂

∂η(η⁰, τ ) −^∂y_∂η(η, τ )

|y(η, τ ) − y(η⁰, τ )|^αdη⁰ (13b)

− θ0

2π Z

Υ1(τ )

∂y₁

∂η(η⁰, τ ) −^∂y_∂η(η, τ )

|y_∗− y(η⁰, τ )|^α dη⁰ (13c) + θ0

2π Z

Υ2(τ )

∂y₂

∂η(η⁰, τ ) −^∂y_∂η(η, τ )

|y_∗− y(η⁰, τ )|^α dη⁰, (13d) for a system of two patches with constant vorticity θ0 and where η is the parametrization of the rescaled curves. The integral in Equation (12) is divergent for α = 1, and also computationally unfavourable when using a pseudo-spectral approach, and hence we project the velocity field on the normal of the curve as done in Equation (13).

By rescaling the space variable, the area of the patches, which is related to the energy of the initial data, are no longer conserved and the system is no longer symplectic. Instead, the area within the patches approaches infinity when singularities are formed. Also, in the original variables, the collapse occurs during a short time whereas in the rescaled variables the time where singularities form is mapped to an infinite time interval which enhances the possibilities to monitor the formation of singularities.

In the search of a singularity in Equation (13), we seek a fixed point, ^∂y_∂τ(η, τ ) = 0. An attractive fixed point would be of importance as it would point in direction towards a rigorous proof of blowup. The infinite area in the rescaled contours made Mancho suggest that there is a stationary profile of the form Υ = (x, y(x)), namely that we can describe the contours as functions on R. For the fixed point to occur, we need the integrals to cancel each other out. More precisely, we need the first and the third integrals to cancel out as well as the second and the fourth. This gives a restriction of the rescaled space variable, and without loss of generality we look at the absolute value of the integral,

Z +∞

M

y⁰(η, ·)

|y∗− y(η, ·)|^α − y⁰(η, ·)

|y(η, ·)|^αdη

≤ Z +∞

M

y⁰(η, ·)

|y∗− y(η, ·)|^α− y⁰(η, ·)

|y(η, ·)|^α

dη,

where we have omitted the subindex, k, and denoted the derivative with respect to its parametrization with y⁰(η, ·). By the inequality |y_∗− y(η, ·)| ≥ |y(η, ·)| − |y_∗|, the integrand, I, has an upper bound as

I ≤ y⁰(η, ·)

(|y(η, ·)| − |y_∗|)^α − y⁰(η, ·)

|y(η, ·)|^α

= |y⁰(η, ·)|

|y(η, ·)|^α

1

(1 − |y_∗|/|y(η, ·)|)^α− 1

. By Taylor expansion in the neighborhood of y∗ we get

I ≤ |y⁰(η, ·)|

|y(η, ·)|^α C|y_∗|

|y(η, ·)|

= C|y_∗| |y⁰(η, ·)|

|y(η, ·)|^1+α.

(16)

Therefore, the criterion for the rescaled variable yk(η, ·) in Equation (13) to attain a fixed point is that yk(η, ·) is such that

η→∞lim

|y_k⁰(η, ·)|

|yk(η, ·)|^1+α = 0.

If a fixed point attracts nearby solutions it would give an indication of local development of singularities of the α-patches problem which possibly could be translated to SQGE and extended to 3DEE.

(17)

3 The pseudo-spectral method

In this section, we present how the discretization of the model is done and on which premises the pseudo-spectral method relies. We also investigate the numerical properties of the method in regards to accuracy by defining the used operations and by presenting a set of suitable spaces.

The error bounds derived and given in this section treat the pseudo-spectral method when using the standard FFT to retrieve the coefficients, hence the results may up to an uncertain degree be useful when determining error bounds for the pseudo-spectral method using the Nonequispaced Fast Fourier Transform.

Further, we give derivations of the frequency expansion step with NFFT in Section3.3 as an option to standard FFT. In Section3.4we introduce the two quadratures used in this thesis, that is Riemann sum and Clenshaw-Curtis quadrature. The pseudo-spectral method assumes that the patch, sampled at a set of collocation nodes, can be expressed in terms of a convergent Fourier series. The choice of Fourier series expansion instead of Chebychev polynomials expansion, for instance, is based on the periodicity of a vortex patch with respect to its parametrization. For a non-periodic problem however, Chebychev polynomials could be used as a basis instead of the trigonometric polynomials in a Fourier series.

With the chosen basis, the velocity field acting on each discrete node is calculated by the use of a quadrature. This is essentially what makes the method pseudo-spectral, rather than Galerkin spectral. The underlying parametrization of a patch is unknown and may not be unique. As the two main components of the pseudo-spectral method, namely basis expansion and quadrature, may rely on different assumptions of the underlying parametrization the final discrete system may vary. For simplicity, we assume in Section3.1that the patches’ parametrization is performed on a uniform grid and that a Riemann sum is chosen as quadrature to evaluate the integral.

3.1 Discretization of the system

Let the parametrization of the patch be over the domain

Ω = (0, 2π], (14)

and let an equidistant mesh over the domain be γi∈ Ω, for i = 0, ..., N . The semi-discrete function values are sampled at collocation points as

x(γ_i, t).

Further, the contour integral over a patch is numerically evaluated with a quadrature of the Fourier coefficients using the nonlinear quadrature weights of the denominator of Equation (8). The evolution of a particle x(γ_j, t) is thus approximately given by

∂x

∂t(γ_j, t) ≈ θ0

2π

N

X

i=1

i6=j

∂x

∂γ(γ_i, t) −^∂x_∂γ(γ_j, t)

|x(γj, t) − x(γ_i, t)|^α,

where the singularity in the integrand at i = j is effectively removed. This results in a 2N × 2N system of partial differential equations which grows with the number of nodes, given by





 x⁽¹⁾_t x⁽²⁾_t ... x^{(N )}_t







= θ1

2π







−PN i=2

1

|x⁽¹⁾−x⁽ⁱ⁾|^α

1

|x⁽¹⁾−x⁽²⁾|^α · · · _|x₍₁₎_−x¹_{(N )}_|_α

1

|x⁽²⁾−x⁽¹⁾|^α −PN

i=1

i6=2 1

|x⁽²⁾−x⁽ⁱ⁾|^α · · · _|x₍₂₎_−x¹_{(N )}_|_α

... ... . .. ...

1

|x^{(N )}−x⁽¹⁾|^α

1

|x^{(N )}−x⁽²⁾|^α · · · −PN −1 i=1

1

|x^{(N )}−x⁽ⁱ⁾|^α











 x⁽¹⁾γ

x⁽²⁾γ

... x^{(N )}γ





 ,

which in short notation is given by

x_t= θ1

2πV (x)x_γ, (15)

where x_tand x_γ represent the derivative of the nodes with respect to time and parametrization.

(18)

3.2 Retrieving the coefficients

With the position of the patch given by x(γ, t), in accordance with Equation (7), we search a basis expansion of x(γ, ·) to use the discrete model given in Equation (15). We also need to define a set of operations, that is projection and interpolation, in suitable spaces of trigonometric polynomials in order to analyze the numerical method’s properties.

For any function f (γ) ∈ L²(Ω), with the domain Ω as in Equation (14), we have the Fourier expansion

f (γ) =

∞

X

k=−∞

fˆ_ke^ikγ, (16)

where the coefficients ˆfk are given by fˆk =

Z

Ω

f (γ, ·)e^−ikγdγ.

We next introduce the space Γ^N which is the space of trigonometric polynomials in γ of degree N or less. There are mainly two sources of errors which are introduced in the pseudo-spectral method, that is the truncation error and aliasing. In order to analyze the consistency and accuracy of the pseudo-spectral method applied on the system in Equation (15) we need to define the L²-projection onto Γ^N, as well as an interpolation between the collocation points. The truncation error will be due to the projection onto an insufficiently large space ΓN and the aliasing error will be explained in terms of the interpolation.

Definition 3.1 (L²-Projection). The L²-projection PN of a function f (γ) ∈ L²(Ω) onto the space ΓN of trigonometric polynomials is given by

Z

Ω

(f − P_Nf )vdγ = 0, ∀v ∈ ΓN.

Since the L²-projection belongs to the space ΓN, we can write it as a linear combination of the trigonometric polynomials spanning ΓN, that is

P_Nf =

N

X

k=−N

fˆ_ke^ikγ.

Further, the L²-projection is the best approximate estimate of the function f , as stated in Theo- rem4.

Theorem 4 (Best approximation result). The L²-projection as defined in Definition3.1 satisfies the best approximation estimate

||f − PNf ||_L2(Ω)≤ ||f − v||_L2(Ω), ∀v ∈ ΓN

Proof. We use the definition of the L²-norm while rewriting f − PNf = f − v + v − PNf , where v ∈ ΓN to obtain

||f − PNf ||²_L2(Ω)= Z

Ω

(f − PNf )(f − v + v − PN)dγ

= Z

Ω

(f − PNf )(f − v)dγ + Z

Ω

(f − PNf )(v − PN)dγ

= Z

Ω

(f − PNf )(f − v)dγ

≤ ||f − PNf ||_L2(Ω)||f − v||_L2(Ω), which concludes the proof.

To show consistency of the pseudo-spectral method applied on the system in Equation (15), we need a bound on how large the projection error is in terms of the exact solution when evaluating the velocity field, i.e. the right hand side of Equation (15). The following theorem is useful for

(19)

bounding the projection error. In order to present the theorem, we introduce proper spaces and norms.

Let the Sobolev norms and seminorms on the set of infinitely differentiable periodic functions on Rⁿ, C_p^∞(Ω), for the positive integers µ and σ be given by

||f ||σ= X

|α|≤σ

||Dαf ||²_L2(Ω)

^1/2 ,

and

|f |µ= X

|α|=µ

||Dα||²_L2(Ω)

^1/2

respectively, where

Dα= ∂^α¹

∂x^α₁¹ . . . ∂^αⁿ

∂x^αnⁿ

.

The Hilbert space H_p^σ(Ω) is the completion of C_p^∞(Ω) under the norm || · ||σ in accordance with Pasciak [26]. For µ, σ ≥ 0, the space H_p^σ(Ω) is given be complex interpolation as stated by Krein and Petunin [19].

Theorem 5 (Projection error bound). For any µ, σ ∈ R such that 0 ≤ µ ≤ σ, the error of projecting a function f ∈ H_p^σ(Ω) on ΓN can be bounded by

||f − PNf ||_µ≤ CN^µ−σ|f |σ.

For a proof, see Canuto and Quarteroni [5]. It follows from the differentiability of f and the fact that the projection and differentiation commute,

d

dγPNf = PN

df dγ,

that the derivative with respect to the parametrization γ of the projection can be bound according to

||∂γ(f − PNf ) ||µ≤ CN^µ−σ|f |σ, f ∈ H_(p)^σ+1(Ω).

Analogously, this applies to the time-derivative as

||∂t(f − PNf ) ||µ ≤ CN^µ−σ|∂tf |σ, f ∈ H_(p)^σ (Ω),

As stated in Theorem 4, the projection onto ΓN gives the best possible estimate of all the functions in ΓN. However, given a set of nodes, or collocation points, that we can use in the estimate, we want to bound the error which results from choosing a sub-optimal set of coefficients.

To estimate coefficients given a certain set of nodes, we define the interpolation between these in Definition3.2.

Definition 3.2 (Fourier interpolation). Let a function f (γ, t) ∈ H_(p)^σ (Ω) be sampled along its parametrization γ. The Fourier interpolation is given by

(I_Nf ) (γ) =

N

X

k=−N

( ˆf_c^N)_ke^2πikγ^L .

The coefficients in the interpolation are those that are retrieved using FFT. Note that the coefficients in the projection are not equal to the coefficients in the interpolation, i.e. ˆfk6= ( ˆf_c^N)k. The difference between the projection and the interpolation is the aliasing error, for which a bound is given in Theorem6 followed by a proof of the special case of one dimension with an L²-bound.

For a proof of a H^r-bound in d dimensions, see Gottlieb [15].

Theorem 6 (Interpolation error bound). Let a function f (γ, t) ∈ Γ^pN for p, N ∈ N+ in d dimensions be interpolated as in Definition 3.2, the error introduced by the interpolation is then bounded by

||(INf )||²_Hr ≤ p^d||f ||²_Hr.

(20)

Proof. For any f ∈ Γ^pN in one dimension, the Fourier series expansion is given by

f (γ) =

pN

X

k=−pN

fˆ_ke^2πikγ^L .

Further, by interpolation of a set of collocation nodes we have

(I_Nf )(γ) =

N

X

k=−N

( ˆf_c^N)_ke^2πikγ^L ,

with the coefficients, showing sign of aliasing,

( ˆf_c^N)_k =







fˆk+ ˆfk+2N +1+ ... + ˆfk+pN +1, −N ≤ k ≤ −1,

fˆk, k = 0,

fˆk+ ˆfk−2N −1+ ... + ˆfk−pN −1, 1 ≤ k ≤ N.

(17)

With the use of Parseval’s inequality, we get a bound for f as in

||f ||²_L2 =

pN

X

k=−pN

| ˆfk|²,

and for the interpolation as

||(INf )||²_L2 =

pN

X

k=−pN

|( ˆf_c^N)k|².

The collocation coefficients in Equation (17) can be bounded by the triangle inequality

|( ˆf_c^N)_k|²=







| ˆfk+ ˆfk+2N +1+ ... + ˆfk+pN +1|²≤ p

| ˆfk|²+ | ˆfk+2N +1|²+ ... + | ˆfk+pN +1|²

, −N ≤ k ≤ −1.

| ˆfk+ ˆfk−2N −1+ ... + ˆfk−pN −1|²≤ p

| ˆfk|²+ | ˆfk−2N −1|²+ ... + | ˆfk−pN −1|²

, 1 ≤ k ≤ N, from which we get

||INf ||²_L2 ≤ p||f ||²_L2, which concludes the proof.

With the numerical scheme presented in Section3.1, and the projection and interpolation errors defined and bounded, we can draw the conclusion that the method is of geometrical convergence for periodic domains and that the model is consistent. The best possible approximation, as stated in Theorem 4 is retrieved by solving a minimization problem to retrieve the coefficients of the projection. However, with the interpolation between the nodes as collocation points, we can obtain a bounded estimate of the velocity field.

3.3 Nonequispaced Fast Fourier Transform

The use of FFT in the pseudo-spectral method, or in any spectral method for that matter, relies on that the signal is sampled uniformly in space or time. In the expansion of Fourier coefficients to a given curve it is assumed that the curve is sampled equidistantly with respect to its parametrization, more precisely that the samples of the coordinates correspond to a uniform grid of their parameter.

However, in the case of evaluating a velocity field of a given curve, a uniform underlying grid might be suboptimal in terms of computational efficiency, as one typically would make use of an adapting grid which puts more nodes near parts of high curvature and closely spaced curves. It is also inevitable that the sampling of the curve can not be done uniformly as the curve evolves in time since nodes move differently depending on their positions in the velocity field. This situation, where obtaining a equispaced grid may be impractical and/or hardly achieved motivates the use of a method which is capable of nonequidistant sampling. The Nonequispaced Fast Fourier Transform (NFFT) is designed to handle the nonequidistant sampling, at the same time retaining the spectral accuracy of FFT, and can incorporate adapting node handling methods for better computational properties.

(21)

In this section, we give a theoretical description of NFFT which has been used when evaluating the velocity field acting on a vortex patch, both in the rescaled and initial variables. The concepts that follow in this section are not particularly new in the sense that the discrete Fourier transform is still evaluated using a divide and conquer approach of complexity O(N log N ), but it includes a series of convolution and deconvolution steps with a window function to cope with the nonequispaced samples. Since we are no longer constrained to have a square system matrix, this yields a computational complexity of order O(N log N + P ), for performing NFFT on a grid of P samples.

The computation of the inverse transform on the other hand becomes in this setting equivocal, because the absence of an inverse of the system matrix may occur and therefore, needs to be treated in a non-direct approach as explained in Section3.3.2. In AppendixCwe provide a tutorial on how to use NFFT in C. The tutorial is written to be understood independently of this thesis, which motivates some degree of overlapping theory.

3.3.1 Forward Nonequispaced Fast Fourier Transform

To evaluate a given function f (γ) ∈ L²(Ω) sampled nonequidistantly in γ as a sum of trigonometric polynomials, we first give a definition of the window function that will be used to deconvolve the Fourier coefficients ˆf into ˆg ← ˆf / ˆϕ. A formal definition of the window function is given by Kunis [20] and restated in Definition3.3.

Definition 3.3 (Window function). Let a window function ϕ : R → C such that it is satisfies

|ϕ(γ)| ≤ A|γ|^−1−α, for some constants α, A ∈ R⁺, be given by the Fourier transform ϕ(γ) :=

Z

R

ω(z)e^−2πiγzdz, (18)

where ω is an even and continuous function satisfying ω : R → R⁺, |ω(z)| ≤ B|z|^−1−β, with

|z| ≤^N₂ for some constants β, B ∈ R⁺ and N ∈ N.

With the window function as in Equation (18), we seek its 1-periodic version

˜

ϕ(γ) =X

r∈Z

ϕ(γ + r),

with an absolute convergent Fourier series,

˜

ϕ(γ) =X

k∈Z

c_k( ˜ϕ)e^−2πikγ,

with the Fourier coefficients

c_k( ˜ϕ) = Z 1/2

−1/2

˜

ϕ(γ)e^−2πikγdγ.

Several window functions can be used as long as these conditions are fulfilled. The window function used in this thesis is the Kaiser-Bessel window as given in Equation (19), and it maximizes the energy in the main lobe. With the definition of the window function in place, we seek a function g such that f can be approximated as a deconvolution of g with the window function, that is

f (γ) ≈ s(γ) = X

l∈I_M

glϕ˜

γ − l

M

,

where M := σN is the oversampled number of frequencies, and σ ∈ N⁺is the oversampling factor.

The Fourier transform of s1gives

s1(γ) =X

k∈Z

ˆ

gkck( ˜ϕ)e^−2πikγ,

where the summation can be split at the frequency support IM to obtain

s₁(γ) = X

k∈IM

ˆ

g_kc_k( ˜ϕ)e^−2πikγ+ X

r∈Z\{0}

M −1

X

k=−M

ˆ

g_kc_{k+M r}( ˜ϕ)e−2πi(k+M r)γ.

Extended analysis of a pseudo-spectral approach to the vortex patch problem

Examensarbete 30 hp Juni 2018

Extended analysis of a pseudo-spectral approach to the vortex patch problem

Mattias Bertolino

Abstract

Extended analysis of a pseudo-spectral approach to the vortex patch problem

Mattias Bertolino

Populärvetenskaplig sammanfattning

Contents

1 Introduction

2 The vortex patch problem

2.1 3D Euler and Surface Quasi-Geostrophic

2.2 Vortex patches

2.3 Contour dynamics equation

2.4 Polar formulation of the contour dynamics equation

C(γ)

C(0) C(γ

)

x(γ

,t) p(γ

,t)

C(γ

)

C(0) x(γ

,t)

p(γ

,t)

2.5 Self-similar problem in rescaled variables

3 The pseudo-spectral method

3.1 Discretization of the system

3.2 Retrieving the coefficients

3.3 Nonequispaced Fast Fourier Transform