Green's function-based time stepping for the Kuramoto-Sivashinsky initial-boundary value problem
GGREEN’S FUNCTION-BASED TIME STEPPING FOR THEKURAMOTO-SIVASHINSKY INITIAL-BOUNDARY VALUEPROBLEM ˚ L. VAN VEEN : Abstract.
Both theoretical and numerical studies of the Kuramoto-Sivashinsky equation havemostly considered periodic boundary conditions. In this setting, the Fourier decomposition of thesolution is central to theoretical ideas, such as renormalization group arguments, as well as to nu-merical solution, allowing for the construction of accurate and efficient time-steppers using standardpseudo-spectral methods. In contrast, fixed boundary conditions induce boundary layers and neces-sitate the use of non-uniform grids, usually generated by orthogonal polynomials. On such bases,numerical differentiation is ill-conditioned and can potentially lead to a catastrophic blow-up ofround-off error. In this paper, we use ideas recently explored by Viswanath (
J. Comput. Phys. (2013), pp. 414-431) to completely eliminate numerical differentiation and linear solving from thetime-stepping algorithm. We use the Green’s function-based method to investigate elements of theKuramoto-Sivashinsky dynamics over a range of five decades of the viscosity.
Key words.
Nonlinear initial-boundary value problem, Kuramoto-Sivashinsky
AMS subject classifications.
1. Introduction.
Boundary conditions form an important part of models ofcontinuous, space and time dependent processes. Prescribing continuous variableslike velocity or concentration, or their derivatives, at material boundaries often leadsto boundary layers, i.e. regions with steep gradients of the dependent variables. Pe-riodic boundary conditions, in contrast, often lead to more spatially homogeneoussolutions. Schemes for numerical simulation must take the peculiar structure of solu-tions, induced by boundary conditions, into account.Periodic boundary conditions arise naturally, and exactly, in some geometries,but are also often used when the actual domain is practically infinite or no more nat-ural choice is available. It is straightforward to list a number of reasons why periodicboundary conditions are pleasant to deal with. Firstly, the natural choice of a grid forspatial discretisation is a regular grid, and the natural choice of a basis to expand thesolution in is one consisting of sines and cosines. The Fast Fourier Transform (FFT)provides an efficient way to switch between grid and spectral representations of asolution. Each basis function satisfies the boundary conditions and is an eigenfunc-tion of any spatial differential operator. Exploiting these facts, we can use standardpseudo-spectral methods to simulate the system. Examples of pseud-spectral methodsfor semilinear partial differential equations, such as the Kuramoto-Sivashinsky (KS)equation, can be found, for instance, in Trefethen’s text book [16]. The essential stepscan be summarised as1. Discretize time using an implicit method for linear terms and an explicitmethod for nonlinear terms. This leads to a periodic Boundary Value Problem(BVP) to solve for each time step.2. Formulate the BVP in Fourier space. The linear differential operator is rep-resented by a diagonal matrix, so that the solution can be written explicitlyin terms of the Fourier transform of the nonlinear terms. ˚ This work was supported by the National Sciences and Engineering Research Council of Canadaand the Japan Society for Promotion of Science : Faculty of Science, University of Ontario Institute of Technology, 2000 Simcoe St. N., Oshawa,ON L1H 7K4, Canada ( [email protected] )1 a r X i v : . [ m a t h . NA ] J un L. van Veen
3. Use the inverse FFT to find the dependent variables and their derivatives ona regular grid, evaluate the nonlinear terms there and use the FFT to findtheir contribution to the BVP.On a fixed number of n grid points in each spatial dimension, this approach yieldsa time step that requires O p n ln n q FLoating point OPperations (FLOPs) for eachspatial dimension and has a bound on the spatial discretization error of the formexp p´ cn q , for come positive constant c , provided that the solution is smooth.With fixed boundary conditions, some complications arise. The Fourier basesand regular grids are no longer optimal, as they leads to spurious oscillations knownas Gibb’s phenomenon. Instead, theory prescribes the use of orthogonal polynomialsand clustered grids generated by their zeros or extrema. These basis functions arenot eigenfunctions of all spatial differential operators and do not usually satisfy theboundary conditions. The most straightforward way to address these issues is toadjust the above scheme as follows:2a Formulate the BVP in spectral space. The linear differential operator isrepresented by a matrix that is structured (e.g. triangular) but not diagonal.2b Delete one linear equation for each boundary condition and replace it by theboundary condition on the polynomial basis.3 Use an appropriate (inverse) fast transform to evaluate the nonlinear termson the clustered grid.4 Solve the resulting linear problem.In step 2b, only a small error is incurred for well-resolved solutions. Furthermore, formany polynomial bases, fast transforms from grid point values to and from expansioncoefficients are know. For Chebyshev bases, for instance, the FFT can be used, asexplained by Trefethen [16], who presents several explicit examples of this approach.The most important issue, however, is the fact that the algorithm now requires asolver for a large linear system which is ill-conditioned for high-order differentiationand fine grids. For differentiation on Chebyshev bases the condition number grows as n p , where p is the order of differentiation. The KS equation is of fourth order andthus the condition number in a naive application of the Chebyshev spectral algorithmdoes not yield a useful bound on the accumulation of errors.These are various ways to improve the conditioning of the linear system. Itappears that the most successful approach combines two elements: preconditioning thelinear system by a spectral integration matrix and expanding the highest derivative ofthe unknown functions, rather than the unknowns themselves, in a polynomial series.Muite [12] compares several subtly different variations of this approach on grids upto 10 ,
000 points. The preconditioned system still has a condition number that growsalgebraically with the number of grid points and is of no use when estimating errors.A detailed examination of the linear systems arising in the various forms of spectralintegration by Viswanath [19] lead to the conclusion that some can yield results nearmachine accuracy in spite of the bad conditioning. This good accuracy hinges on thecancellation of errors and the careful implementation of the boundary conditions.In the current paper, we follow a different, more radical approach, which wasalso investigated by Viswanath [18]. We solve the BVP analytically using Green’sfunction, thereby eliminating the need for numerical differentiation and linear solvingaltogether. Instead, we must construct a proper quadrature, that has an exponentiallysmall error in spite of the finite differentiability of Green’s function.We show, that a combination of barycentric re-interpolation to sub grids andlocal Clenshaw-Curtis quadrature leads to an accurate time-stepping algorithm that reen’s function based time-stepping
2. The KS IBVP.
We consider the following IBVP, originating in the work ofKuramoto [10] and Sivashinsky [15]: u t ` uu x ` u xx ` νu xxxx “ u p´ q “ l, u p q “ r, u xx p´ q “ u xx p q “ u and u x at the boundaries. In that case, the time-stepping algorithm described here does notchange, but Green’s function does, i.e. the analysis presented in Appendix A must beadjusted. In the literature, the KS equation is often presented in the scaling ¯ u “ ? ν u ,¯ x “ p x ` q{? ν and ¯ t “ t { ν , which yields an equation identical to (2.1) but with ν ” r , L s “ r , {? ν s . This scaling is used, for instance, in the numericalwork of Wittenberg and Holmes [20] and the theoretical work of Galaktionov et al. [6]. The former work provides an overview of typical dynamics generated by the KSequation and a rich reference list. In the latter work, the existence of a boundedsolution for any finite time of this IBVP is proven.Another common guise is the integral formulation h t “ p h x q ´ h xx ´ νh xxxx (2.3)where ´ h x “ u . This form is often used in the physics literature when consideringthe KS equation as a model for interface growth, for instance in Refs. [22, 4]. Forour time-stepping algorithm it is more convenient to solve for u on the fixed spatialdomain and consider the viscosity ν as control parameter. We then typically observea transition from equilibrium at large values of ν to time-periodic motion and finallyspatio-temporal chaos as ν is decreased. The spatio-temporal chaos exhibits a formof extensivity, as demonstrated in section 6.
3. Time-stepping based on Green’s function.
In this section, we will firstdescribe the Green’s function-based algorithm on a high level, before specifying andjustifying our choice of grids, interpolation and quadrature rules.It is convenient to consider deviations from a linear profile: v t ` vv x ` R v ` φv x ` R φ ` v xx ` ν v xxxx “ v p´ q “ v p q “ v xx p´ q “ v xx p q “ u “ v ` φ, φ p x q “ l ` R p x ` q , R “ r ´ l L. van Veen ∆ k “ k “ k “ k “ o “ h α k β k o “ h { α k { ´ { β k ´ o “ h { α k { ´ {
11 2 { β k ´ o “ h { α k { ´ {
25 16 { ´ { β k ´ ´ Table 3.1
Table of coefficients of SBDF formula (3.4) for orders up to o “ . The first step is to turn the IBVP intoa linear BVP for every time step. This is done by the use of an implicit-explicit timediscretization. In particular, we use a Semiimplicit Backward Differentiation Formula(SBDF) [2] to obtain L v k ` “ o ÿ s “ α s v k ` ´ s ` ∆ o ÿ s “ β s f p v k ` ´ s q (3.4)where ∆ is a constant that depends on the time step size h and on o , the order of theSBDF formula. The sub scripts denote the approximate solution at different times.We have introduced the linear operator L and the part of the BVP that is treatedexplicitly, f , according to L “ ` ∆ R ` ∆ B xx ` ∆ ν B xxxx (3.5) f p v q “ ´ vv x ´ φ p x q v x ´ Rφ p x q (3.6)The definitions of ∆ and the coefficients α and β are listed in Table 3.1.The time-discretized equation (3.4), together with boundary conditions (3.2), con-stitutes a linear BVP since all quantities on the right-hand side are known explicitly.The solution can be written in terms of Green’s function as v k ` “ o ÿ s “ α s G ˚ v k ` ´ s ` ∆ o ÿ s “ β s G ˚ f p v k ` ´ s q (3.7)where the star denotes the convolution p G ˚ v qp x q “ ż y “ G p x, y q v p y q d y This expression is not suitable for numerical quadrature because the second convolu-tion contains derivatives of the unknown function. We apply integration by parts toobtain v k ` “ o ÿ s “ p α s ` β s ∆ R q G ˚ v k ` ´ s ` ∆ o ÿ s “ β s DG ˚ ˆ v k ` ´ s ` φv k ` ´ s ˙ ´ ∆ RG ˚ φ ” G ˚ I ` DG ˚ I ` J (3.8) reen’s function based time-stepping DG denotes the derivative of Green’s function with respect to the variable ofintegration of the convolution and, for later convenience, we have introduced short-hand notation for the terms that appear in the convolutions with G and DG and theconstant term. Now the only derivative remaining can be computed analytically so theproblem of the bad conditioning of numerical differentiation has been eliminated. Inits place, we face two new challenges. Firstly, we must compute Green’s function andcast it in a form suitable for numerical evaluation. In Appendix A, appropriate ex-plicit expressions are derived. Secondly, we must accurately compute the convolutionsin the knowledge that the integrand is only once or twice continuously differentiable.For this end, we will use polynomial re-interpolation followed by standard quadrature. When using classical, global quadra-ture rules, the finite differentiability of the integrands would lead to a fixed rateof convergence, much like a standard finite difference method for solving the IBVP(2.1,2.2) would have given us. We can solve this issue by using separate quadraturesfor the sub domains r´ , x s and r x, s on which G p x, y q is smooth. This will requirereinterpolation of the integrand onto suitable, non-uniform sub grids.Let us denote the global grid by x i , i “ , . . . , n and the left and right sub grids forthe i th global grid point by x l ,ij , j “ , . . . , n l ,i and x r ,ij , j “ , . . . , n r ,i . In addition,we will denote the approximations of a quantity a p x q on these grids by a , a l ,i and a r ,i .Then the convolutions with G and DG are represented by matrix-vector products b “ G ˚ a Ñ b “ M a M ij “ n l ,i ÿ p,q “ C l ,ip G l ,ipq B l ,iqj ` n r ,i ÿ p,q “ C r ,ip G r ,ipq B r ,iqj b “ DG ˚ a Ñ b “ N a N ij “ n l ,i ÿ p,q “ C l ,ip s G l ,ipq B l ,iqj ` n r ,i ÿ p,q “ C r ,ip s G r ,ipq B r ,iqj (3.9)where ‚ B l ,i is a p n l ,i ` q ˆ p n ` q matrix of interpolation from the global grid tothe i th left sub grid, and likewise for B r ,i , ‚ G l ,i and s G l ,i are the p n l ,i ` q ˆ p n l ,i ` q diagonal matrices with the values G p x i , y q and B y G p x i , y q take on the i th left sub grid, i.e. G l ,ipq “ δ pq G p x i , x l ,ip q ¯ G l ,ipq “ δ pq B G p x, y qB y ˇˇˇˇ p x i ,x l ,ip q (3.10)and likewise for G r ,i and s G r ,i and ‚ C l ,i is a row vector of quadrature coefficients such that C l ,i a l ,i « x i ż x “´ a p x, t q d x and likewise for C r ,i . Here, we will give explicit expressions and assessthe expected degree of accuracy for the combination of barycentric interpolation onclosed Chebyshev grids and Clenshaw-Curtis quadrature.
L. van Veen
We will take both the global grid and the sub grids to consist of Chebyshev pointsof the second kind, i.e. x i “ cos ˆ iπn ˙ , i “ , . . . , nx l ,ij “ p ψ l ,i ˝ cos q ˆ jπn l ,i ˙ , j “ , . . . , n l ,i with ψ l ,i p x q “ ´ ` p x i ` qp x ` q x r ,ij “ p ψ r ,i ˝ cos q ˆ jπn r ,i ˙ , j “ , . . . , n r ,i with ψ r ,i p x q “ p ´ x i qp x ´ q ` B l ,ipq “ ˜ n ÿ k “ w k p x l ,ip ´ x q q w q p x l ,ip ´ x k q ¸ ´ with w q “ $’&’% { q “ p´ q q for 0 ă q ă n p´ q n { q “ n (3.12)and likewise for B r ,i . The extremal points of the local grids coincide with global gridpoints. This leads to a singularity in the formula above, so that we must separatelyspecify that B l ,i q “ δ qi B r ,i q “ δ q B l ,in l ,i q “ δ qn B r ,in r q “ δ qi The Clenshaw-Curtis row vector is a product of a row vector of quadrature weightswith a matrix representing the discrete cosine transform: C l ,ip “ p x i ` q a p n l ,i ÿ k “ c k F kp where a p “ { p “ , n l ,i c k “ $’&’% k “ {p ´ k q if k is even and k ą
00 if k is odd F kp “ n l ,i cos ˆ kpπn l ,i ˙ (3.13)and likewise for C r ,i .The number points in the local grids, n l ,i and n r ,i , should be as least as large asthe number of points in the part of the global grid they span, i.e. n l ,i ą n ` ´ i and n r ,i ą i `
1. Under this condition, the error of interpolation from the global to thelocal grids is negligible as compared the the error of interpolation of the solution ontothe global grid. Importantly, barycentric interpolation is stable under over resolution.If we set n l ,i “ n r ,i “ n `
4. Implementation.
The task of time-stepping the IBVP can now be split upinto two steps. In the first step, we compute the quadrature matrices M and N forgiven ν , ∆ and n and a specific choice of the number of points in the local Chebyshev reen’s function based time-stepping p “ , . . . , P ´
1, and process number p will computethe solution on n p points with indices i p s through i p e . We further assume that P ! n sothat the work load can be almost evenly balanced, i.e. n p « n { P for p “ , . . . , P ´ p n q , bringing the total for this algorithm to O p n q . In the parallel implementationeach processor will thus handle O p n { P q .One more remark about the quadrature algorithm must be made. If the order ofthe global grid is large, it may happen that it coincides with points on local grids upto machine accuracy. When using standard double precision arithmetic, for instance,this start to happen for global grid orders upward from 20,000. This leads to divisionsby zero when evaluating (3.12). This can simply be circumvented by detecting overlapup to finite precision of the grids and replacing each row of B that holds the coefficientsof interpolation onto a overlapping local grid point by a row of zero elements and asingle element equal to unity, just like for the extremal points of the local grid. Itis a remarkable fact, explained in detail by Higham [7], that the evaluation of theelements of B is otherwise stable, in spite of the small denominators. In our test ofthe accuracy of the numerical quadrature in section 5, we replace rows correspondingto local grid points that are closer to global grid points than 10 ´ and the resultingcomputation is stable up to a global grid order of at least 74,000.Algorithm 2 describes the time-stepping. The loop over time steps includes thecomputation of the integrands in equation (3.8), which takes O p n p q FLOPS on proces-sor p , an MPI all-to-all communication of O p n q elements and a matrix vector productwith the elementary FLOP count of O p n p n q . If the MPI routine takes a similaramount of time to complete as O p n q FLOPS, it is obvious that the communicationtime will be negligible as long as n p "
1, as we have assumed.Figure 8 shows the wall time taken for building the quadrature matrices andtime-stepping 4000 times for grid sizes 1000 and 2000. This test was run on a clus-ter computer, on a node with 24 AMD Opteron 2.2GHz processors, 32Gb of RAMmemory, 512Kb of cache and a QDR InfiniBand connection. Clearly, Algorithm 1scales linearly as predicted by the FLOP count since it does not involve any commu-nication. Time-stepping according to Algorithm 2, on the other hand, shows approx-imately linear speed-up up to 9 processors only. This results depends on architecture,on a machine with two quad-core CPUs, for instance, the speed-up saturates at 2processors. There are two possible reasons for this limitation. Firstly, the executionof the all-to-all communication depends on the details of the hardware and the MPIimplementation and is highly sensitive to concurrent use by other processes of cachememory. Secondly, the all-to-all communication is blocking and the execution timefor the partial matrix-vector product may vary over processes.Lastly, we turn to the nontrivial question of initializing the time-stepping. ForSBDF orders two and up, we need the solution at previous time instants. We propose
L. van Veen
ALGORITHM 1:
Computation of the quadrature matrices and constant term
Input : viscosity ν , time-stepping parameter ∆, boundary values l and r and globaland local Chebyshev grid orders n , n L ,i , n R ,i , i “ , . . . , n . Output : p n ` q ˆ p n ` q quadrature matrices M and N and p n ` q vector J ,distributed by rows over P processes. for i “ i p s . . . i p e do % loop over rows stored on processor p Allocate a l j , b l j , j “ , . . . , n l ,i , a r j , b r j , j “ , . . . , n r ,i . Set x l ,ij , j “ , . . . , n l ,i and x r ,ij , j “ , . . . , n r ,i according to (3.11). Compute the Clenshaw-Curtis row vectors according to (3.13) and set a l Ð C l ,i b l Ð C l ,i a r Ð C r ,i b r Ð C r ,i Multiply by Green’s function or its derivative, i.e. set a l Ð n l ,i ÿ p “ a l p G l ,ipq b l Ð n l ,i ÿ p “ b l p s G l ,ipq a r Ð n r ,i ÿ p “ a r p G r ,ipq b r Ð n r ,i ÿ p “ b r p s G r ,ipq Compute the interpolation matrices according to (3.12) and set M ij Ð n l ,i ÿ p “ a l p B l ,ipj ` n r ,i ÿ p “ a r p B r ,ipj N ij Ð n l ,i ÿ p “ b l p B l ,ipj ` n r ,i ÿ p “ b r p B r ,ipj for j “ , . . . , n De-allocate a l j , b l j , a r j , b r j . Compute J i “ ´ ∆ R ř nj “ M ij φ p x j q . return M p , N p and J p , i.e. rows i p s through i p e of the quadrature matrices and theconstant term. ALGORITHM 2:
Green’s function based time-stepping
Input : time step size h , number of inner and outer iterations N i , N o , o initial points v ´ o ` , . . . v , boundary values l and r . M p , N p and J p are stored on processor p . for j “ , . . . , N o do MPI scatter v k ´ o ` i Ñ v k ´ o ` i . . . v p ´ k ´ o ` i for i “ , . . . , o % root to all for k “ , . . . , N i do Compute the integrands I p and I p of (3.8). MPI gather the integrands: I , , . . . , I p ´ , Ñ I , . % all to all Set v pk ` “ M p I ` N p I ` J p MPI gather v k , . . . v p ´ k Ñ v k . % all to root Root: output v k ` φ . four possible solutions:1. Using the first order SBDF formula with a small time step. This is a com-monly used method to seed backward differentiation formulae, and is em-ployed by Ascher et al. [2]. For the Green’s function based time-stepper,however, this method has its limitation. If the time step is taken very small reen’s function based time-stepping w a ll t i m e ( s ) n = n = n = n = √ ν u ( x , t ) -1 -0.5 0 0.5 1 x t / ν -4-2024 Fig. 4.1 . Left: wall time for the computation of the quadrature matrices according to Algorithm1 (red) and for time stepping according to Algorithm 2 (blue) using the first order SBDF. Shown isthe wall time averaged over 20 trial for 4000 time steps with grid orders n “ and n “ .The black line indicates linear scaling. Right: density plot of the solution computed in this test. Theviscosity is set to ν “ ˆ ´ , the step size to h “ ´ and the initial condition is u “ ´ x {? ν . for fixed viscosity, Green’s function approaches a delta distribution as itslength scale, b ´ (see (A.10)) goes to zero. If this scale becomes comparableto the spacing of the local grid near its end points, of order n ´ , an instabilitycan occur.2. Richardson extrapolation from lower order. For instance, we can seed thesecond order SBDF by approximating u p x, h q first by two first order steps ofsize h {
2, giving u h { p x, h q , then by a single step, giving u h p x, h q and finallysetting u p x, h q « u h { p x, h q ´ u h p x, h q . Similar expressions can readily befound for higher order seeding. This method has two disadvantages. Firstly,it relies on the cancellation of error terms , which is likely unstable for thehigher order versions, so that round-off error is introduced. Secondly, for eachstep with a different value of ∆, we need to repeat Algorithm 1, which hasorder O p n q complexity. Therefore, this method is only practical for smallgrids and low SBDF order.3. Using a known exact solution to the KS equation. An exact solution can befound, for instance, in Parkes and Duffy [13]. We use it in section 5 to test theaccuracy of the SBDF formulae. Strictly speaking, this is not a solution to theIBVP, but for small enough viscosity the boundary conditions are satisfiedfar beyond machine accuracy.4. Growing a solution from a small perturbation to the zero solution. If wecompose the perturbation out of eigenmodes of the linear part of (3.1), wecan compute the solution backward in time under the assumption that thenonlinear term in negligible. The disadvantage of this method is that itrequires a long time integration for the perturbations to grow to finite size.This strategy is used in our computation of the finite-size effects in section 6.
5. Convergence and stability tests.
We present two tests to evaluate theaccuracy of the Green’s function based time stepping. The first test demonstrates theexponential convergence of the quadrature computed in each time step. We used the0
L. van Veen − − − − e q n/k k = 6 k = 60 k = 600exp( − αn/k ) Fig. 5.1 . Error of the numerical approximation of the convolution with Green’s function. Testfunction (5.1) is considered with , and oscillations in r´ , s . The solid line denotes thetheoretically expected error of the interpolation and Clenshaw-Curtis quadrature. test function ξ p x, k q “ ` sin p πkx q ´
12 cos p πkx q ´
12 (5.1)which satisfies the homogeneous Dirichlet boundary conditions (3.2) and describes 2 k oscillations on the domain. A similar function was used by Trefethen to demonstratethe convergence of standard spectral methods for functions that can be analyticallycontinued in a neighbourhood of r´ , s in the complex plane [16]. The continuationof our test function has poles at ˘ i ln p ` ? q{p kπ q that determine an upper boundfor the error of Lagrange interpolation on the global grid [3], as follows:max x Pr´ , s | ξ p x, k q ´ p n p x q| ď C exp ´ ´ αk n ¯ where C is a positive constant and p n the interpolant of order n . For k ě
6, we findthat α « .
28. In fact, α “ ln p ` ? q{ π up to corrections of order O p { k q . Thisinterpolation error sets, in turn, an upper bound for the error in the Clenshaw-Curtisquadrature, the difference being a constant factor [17].Let ξ represent the test function evaluated on the global grid, and let ω represent L ξ evaluated on the global grid. Then we measure the quadrature error e q “ } ξ ´ M ω } as a function of the grid order n for fixed ν , h and k . In the first test, we set ν “ ν “ ˆ ´ , h “ ∆ “ ´ and k “ ν “ ν “ ˆ ´ , h “ ∆ “ ´ and k “
60 and in the third ν “ ν “ ˆ ´ , h “ ∆ “ ´ and k “ ? ν forsmall viscosity, as demonstrated in section 6.The three data sets collapse onto a single curve if we plot e q as a function of n { k . Along this curve, the quadrature error decreases approximately as exp p´ αn { k q reen’s function based time-stepping − − − − − − − − k u ( x , T ) − w ( x , T ) k k w ( x , T ) k h/ν h h h h h − o = 1 o = 2 o = 3 o = 4 √ ν u xtν = 0 tν = 8 Fig. 5.2 . Left: relative error in u after time-stepping over t { ν “ for SBDF formulae of orderup to four. The expected error is indicated by a solid line for each order. Right: initial and finalsolution in this test. as predicted. This results indicates that the main error introduced by the numericalevaluation of the convolutions is that of polynomial interpolation. We can make twofurther observations. Firstly, the method is stable to over resolution, as the errordoes not increase beyond the minimum around n { k “ n { k Ç w p x, t q “ c ` ? ? ν „
11 tanh ˆ q ? ν p x ´ ct ´ x q ˙ ´ ˆ q ? ν p x ´ ct ´ x q ˙ (5.2)Here, c is the speed of a solitary wave connecting two constant solutions, x is its initialposition and q “ a {
76 is constant. Of course, w is is not an exact solution to theinitial boundary value problem, but it approaches its left and right limit values at arate of exp p´ qd {? ν q a distance d away from the soliton. We have set ν “ ˆ ´ , c “ x “ ´ . w at the boundaries, and the themagnitude of its second derivative there, are below 10 ´ for 0 ď t { ν ď
8. A similartest was used by Anders [1] et al. , but their viscosity is two orders of magnitude largerand, consequently, they were forced to consider time-dependent boundary conditions.Figure 5(left) shows the error of time-stepping with SBDF formulae up to orderfour. In these tests, the Chebyshev grid order was fixed to n “ k “
6, and consequently we do not expect the interpolation error to play arole. The smallest error achieved is instead determined by the accumulation of round-off error, inversely proportional to the step size h , as indicated with a solid line. Theinitial and final condition in this test are shown in figure 5(right), which shows onlythe centre of the domain.Finally, we consider the stability of the time-stepping method. Figure 5 illustratedthe stability for grid orders n “ n “ L. van Veen − − − − − − − − − h ν ℑ ( p ± ) = stableunstable n = 1000 o = 1 o = 2 o = 3 o = 4 − − − − − − − − − h ν ℑ ( p ± ) = stableunstable n = 2000 o = 1 o = 2 o = 3 o = 4 Fig. 5.3 . Stability of the time-stepping algorithm for SBDF formulae of order up to four. Thealgorithm is stable below the solid lines connecting the symbols for each order and to the right of thedashed line. Left: grid order n “ , right n “ . The solid line at the top of the figures denotesthe transition from imaginary to real-valued roots A.5 of the eigenvalues of the linear operator. Thedashed line indicates the point where ? ν “ { n , i.e. the largest grid spacing is roughly equal to thetypical length scale of the solution. See text for initial condition. randomly chosen amplitudes. A combination of viscosity and time step size waslabeled stable if the integration proceeded up to t { ν “
150 without blow-up. There aretwo boundaries to the stable regime. One lies close to the transition from imaginaryto real roots of the eigenvalues of the linear operator given in A.5. If these rootsare real-valued, Green’s function exhibits global oscillations, meaning that updates tothe solution become dependent over an arbitrarily long distance in a single time step,which renders the step unstable. The other occurs when the typical length scale ofthe solution, expected to scale as ? ν , equals the maximal grid spacing, which is fixedin these experiments.
6. Example computations: finite-size effects.
To demonstrate the power ofthe time-stepping method based on Green’s function and Clenshaw-Curtis quadrature,we generated time series, seeded with random initial conditions of small amplitude,for a range of five orders of magnitude of the viscosity. Each time series extends up to t { ν “ t { ν “ ν “ ´ , we have enlarged one tenth of the domain. The fact that the dynamics lookqualitatively the same as that for ν “ ´ is indicative of scaling behaviour, i.e. forsmall enough viscosity the solutions look similar in the scaled variables ¯ u , ¯ x and ¯ t introduced in section 2.Figure 6 shows the deviation from the time-mean solution near the left boundaryin scaled variables to illustrate the dependence of the boundary layer thickness onviscosity. The fact that the curves approximately overlap indicates that the thick-ness of the boundary layer, in which the boundary conditions strongly influence thedynamics, is about 12 ? ν and is constant in scaled variables on the domain r , L s .
7. Discussion.
We have demonstrated, that the Green’s function based method,in conjunction with the SBDF formulae, barycentric interpolation and Clenshaw-Curtis quadrature, is accurate, stable and reasonably fast for values of the viscosityas small as 10 ´ . To the best of our knowledge, no other time-stepping method forthe KSIBVP has been tested for a viscosity this small or, equivalently, a domain this reen’s function based time-stepping ´ are reliable even for statistical analysis or curve-fitting exercise like thatin the work cited above.In later work, various global and piecewise collocation methods were tested onthe KSIBVP. Often, validation was only done for smooth, viscous solutions, like inKhater & Temsah [9], who used spectral integration on a Chebyshev polynomial basis.Fornberg [5] applied a Chebyshev pseudospectral method, implementing the boundaryconditions in real space, and computed chaotic solutions at ν « ˆ ´ . Piecewisecollocation was tested for a viscosity of order O p ´ q by Xu & Shu [21] and Anders et al. [1]. Based on careful testing and comparison to a priori error estimates, theycould conclude that the spatio temporal chaos they observe numerically is a genuineproperty of the KS equation rather than a numerical artifact, but it remains unclearif their methods are suitable for time-stepping at smaller viscosity.The limiting factor of the Green’s function based method described here is thememory requirement, as it requires storing two p n ` q ˆ p n ` q matrices – albeitdistributes over processors – and the limited scaling of the parallelization. One pos-sible solution is to switch to piecewise Chebyshev grids to avoid excessive clusteringof grid points near the boundaries. This will require a marginally more complicatedprocedure to compute the quadrature matrices and is work in progress. Acknowledgements.
I would like to thank the Faculty of Engineering Science ofOsaka University and the Japan Society for Promotion of Science for making possiblethe sabbatical visit during which much of this work was completed.4
L. van Veen
Appendix A. Computation of Green’s function.
We are looking for Green’s function for the following linear BVP L v “ r ` ∆ R ` ∆ B xx ` ∆ ν B xxxx s v “ r v p´ q “ v p q “ v xx p´ q “ v xx p q “ r will be set by all terms treated explicitly in the time discretization. Withthese boundary conditions, L is symmetric, and has the following spectrum w o k “ sin p kπx q λ o k “ ` ∆ R ´ ∆ π k ` ∆ νπ k (A.2) w e k “ cos ˆ„ k ´ πx ˙ λ e k “ ` ∆ R ´ ∆ π ˆ k ´ ˙ ` ∆ π ν ˆ k ´ ˙ (A.3)where the superscripts denote odd and even. We use the eigenfunction decompositionof Green’s function, given by G p x, y q “ ÿ k “ " w o k p x q w o k p y q λ o k ` w e k p x q w e k p y q λ e k * (A.4)The summations over the odd and even contributions proceeds in a similar fashion,so we will focus on the former. After factorising the denominator as λ o k “ ∆ ν p π k ´ p ` qp π k ´ p ´ q ; p ˘ “ ν ˘ ν a ∆ ´ ν p ` ∆ R q (A.5)we can expand the summation as ÿ k “ w o k p x q w o k p y q λ o k “ ν p p ` ´ p ´ q ÿ k “ " cos p kπ r x ´ y sq π k ´ p ` ´ cos p kπ r x ´ y sq π k ´ p ´ ´ cos p kπ r x ` y sq π k ´ p ` ` cos p kπ r x ` y sq π k ´ p ´ * “ . . . “ (cid:60) ˆ ν p p ` ´ p ´ q ! q f ` p x ´ y q ´ q f ´ p x ´ y q ´ q f ` p x ` y q ` q f ´ p x ` y q )˙ (A.6)where the ellipsis corresponds to some tedious manipulations of the sums to bringthem into the form of the elementary inverse semi discrete Fourier transform q f ˘ p x q “ ˆ π k ´ p ˘ ˙q “ ÿ k “´8 ˆ π k ´ p ˘ ˙ e iπkx “ ´ cos p p ˘ ´ p ˘ | x |q p ˘ sin p p ˘ q (A.7)As can be seen from this expression, Green’s function will exhibit global oscillationsif the discriminant in Eq. A.5 is positive so that at least one of p ˘ is real-valued. Inthat case, the resulting time-stepping scheme will be inaccurate and often unstable,as demonstrated in Sec. 5. We will therefore assume that (cid:61) p p ˘ q ‰
0, which meansthat we impose an upper bound on the time step.Combining the odd and even contributions, we obtain Green’s function in thecompact form G p x, y q “ ν (cid:61) p p ` ´ p ´ q (cid:61) ˆ ´ cos p p ` ´ p ` | x ´ y |q ` cos p p ` p x ` y qq p ` sin p p ` q ˙ (A.8) reen’s function based time-stepping p (cid:61) p p ` qq near the bound-ary, and p ` , in turn, grows as 1 {? ν . This causes large cancellation errors near theboundaries. A more suitable form is G p x, y q “ Q ´ e ´ b | x ´ y | sin p a ´ a | x ´ y | ´ φ q ´ e ´ b ` b | x ´ y | sin p a ´ a | x ´ y | ` φ q´ e ´ b ` b p x ` y q sin p a p x ` y q ´ φ q ` e ´ b ´ b p x ` y q sin p a p x ` y q ` φ q ¯ (A.9)where we have introduced the auxiliary parameters S “ ∆4 ν r ` ∆ R s (0 ă S ă a “ (cid:60) p p ` q “ S ´ { ? ν b ` ? Sb “ (cid:61) p p ` q “ S ´ { ? ν b ´ ? Sθ “ Arg p p ` q “ arctan d ´ ? S ` ? S (0 ă θ ă π { φ “ Arg ˆ p ` sin p p ` q ˙ “ a ´ θ ´ π ` arctan ˆ sin p a q e b ´ cos p a q ˙ Q “ ? ν ∆ d S ? S ´ S a ´ e ´ b cos p a q ` e ´ b (A.10)Where p ` has been chosen to lie in the first quadrant in the complex plane.In this form, Green’s function no longer has exponentially large factors. It isimmediately clear that, if the ratio between ν and ∆ is fixed, then the amplitude of G grows only as 1 {? ν for small viscosity. However, there are still terms of O p q thatcancel near the boundary. To avoid this, we rewrite Green’s function near x “ G p x, y q “ Q sin p a p x ´ qq ” cos p a p y ` q ` φ qp e b p x ´ y ´ q ` e ´ b p x ` y ` q q´ ı cos p a p y ` q ´ φ qp e ´ b p x ´ y q ` e b p x ` y ´ q q ı ´ Q sinh p b p x ´ qq cos p a p x ´ qq ” sin p a p y ` q ´ φ q e b p y ´ q ` sin p a p y ` q ` φ q e ´ b p y ` q ı (A.11)and use the latter form for numerical evaluation if 1 ´ x ă { b .Green’s function satisfies G p x, y q “ G p y, x q G p´ x, ´ y q “ G p x, y q (A.12)the latter property being a consequence of an S symmetry of BVP (A.1), namely p v, x, t q Ñ p´ v, ´ x, t q The BVP is equivariant under this reflection if, and only if, the original boundaryconditions are, i.e. if r “ ´ l in (2.1)-(2.2). As a consequence of (A.12), we only haveto evaluate Green’s function in on the domain 0 ď x ď ´ x ď y ď x . Similarexpressions are readily derived for G y p x, y q near the boundaries x “ y “ ´ L. van VeenREFERENCES[1]
D. Anders, M. Dittmann, and K. Weinberg , A higher-order finite element approach to theKuramoto-Sivashinsky equation , Z. Angew. Math. Mech., 92 (2012), pp. 599–607.[2]
U. M. Ascher, S. J. Ruuth, and T. R. Wetton , Implicit-explicit methods for time-dependentpartial differential equations , SIAM J. Numer. Anal., 32 (1995), pp. 797–823.[3]
J.-P. Berrut and L. N. Trefethen , Barycentric Lagrange interpolation , SIAM Rev., 46(2004), pp. 501–517.[4]
M. C. Cross and P. C. Hohenberg , Pattern formation outside of equilibrium , Rev. Mod.Phys., 65 (1993), pp. 851–1112.[5]
B. Fornberg , A pseudospectral fictitions point method for high order initial-boundary valueproblems , SIAM J. Sci. Comput., 28 (2006), pp. 1716–1729.[6]
V. A. Galaktionov, E. Mitidieri, and S. I Pohozaev , Existence and non-existence of aglobal solution to the Kuramoto–Sivashinsky equation , Doklady Math., 77 (2008), pp. 238–242.[7]
N. J. Higham , The numerical stability of barycentric Lagrange interpolation , IMA J. Numer.Anal., 24 (2004), pp. 547–556.[8]
M. Kardar, G. Parisi, and Y.-C. Zhang , Dynamics scaling of growing interfaces , Phys. Rev.Lett., 56 (1986), pp. 889–892.[9]
A. H. Khater and R. S. Temsah , Numerical solutions of the generalized Kuramoto–Sivashinsky equation by Chebyshev spectral collocation methods , Comput. Math. Appl.,56 (2008), pp. 1465–1472.[10]
Y. Kuramoto and T. Tsuzuki , Persistent propagation of concentration waves in dissipativemedia far from thermal equilibrium , Prog. Theor. Phys., 55 (1976), pp. 356–369.[11]
P. Manneville , Statistical properties of chaotic solutions of a one-dimensional model for phaseturbulence , Phys. Lett., 84A (1981), pp. 129–132.[12]
B. K. Muite , A numerical comparison of Chebyshev methods for solving fourth order semilin-ear initial boundary value problems , J. Comput. Appl. Math., 234 (2010), pp. 317–342.[13]
E. J. Parkes and B. R. Duffy , An automated tanh-function method for finding solitary wavesolutions to nonlinear evolution equations , Comp. Phys. Commun., 98 (1996), pp. 288–300.[14]
H. Sakaguchi , Shock structures in time–averaged patterns for the Kuramoto–Sivashinsky equa-tion , Phys. Rev. E, 62 (2000), pp. 8817–8819.[15]
G. Sivashinsky , Nonlinear analysis of hydrodynamic instability in laminar flames I. derivationof basic equations , Acta Astron., 4 (1977), pp. 1177–1206.[16]
L. N. Trefethen , Spectral methods in Matlab , SIAM, 2000.[17] ,
Is Gauss quadrature better than Clenshaw-Curtis? , SIAM Rev., 50 (2008), pp. 67–87.[18]
D. Viswanath , Navier-Stokes solver using Green’s functions I: Channel flow and plane Couetteflow , J. Comput. Phys., 251 (2013), pp. 414–432.[19] ,
Spectral integration of linear boundary value problems , J. Comput. Appl. Math., 290(2015), pp. 159–173.[20]
R. W. Wittenberg and P. Holmes , Scale and space localization in the Kuramoto-Sivashinskyequation , Chaos, 9 (1999), pp. 452–465.[21]
Y. Xu and C.-W. Shu , Local discontinuous Galerkin methods for the Kuramoto–Sivashinskyequations and the Ito-type coupled KdV equations , Comput. Method. Appl. M., 195 (2006),pp. 3430–3447.[22]
V. Yakhot , Large-scale properties of unstable systems governed by the Kuramoto–Sivashinskyequation , Phys. Rev. A, 24 (1981), pp. 642–644.reen’s function based time-stepping √ ν u ( x , t ) -1 -0.5 0 0.5 1 x t / ν -4-2024 √ ν u ( x , t ) -1 -0.5 0 0.5 1 x t / ν -4-2024 √ ν u ( x , t ) -1 -0.5 0 0.5 1 x t / ν -4-2024 √ ν u ( x , t ) -1 -0.5 0 0.5 1 x t / ν -4-2024 √ ν u ( x , t ) -0.1 -0.05 0 0.05 0.1 x t / ν -4-2024 Fig. 6.1 . Fragment of the time series of u for five orders of magnitude of the viscosity. Fromtop to bottom: ν “ ´ ( n “ ), ν “ ´ ( n “ ), ν “ ´ ( n “ ), ν “ ´ ( n “ ) and ν “ ´ ( n “ ). Shown are the contours of ¯ u “ ? νu for ď t { ν ď .The initial condition in each case is a small random perturbation of the zero solution, the SBDForder is o “ and the time step is h { ν “ ˆ ´ . For ν “ ´ , an enlargement of one tenth ofthe domain is shown so that the qualitative dynamics can be compared to that at ν “ ´ . L. van Veen √ < ¯ u − < ¯ u >> ¯ x ν = 10 − ν = 10 − ν = 10 − ν = 10 − ν = 10 − Fig. 6.2 . Root-mean-square of the departure of the solution from the mean profile near theboundary. The brackets denote averaging from t { ν “ to t { ν “ and the scaled variables are ¯ u “ ? νu and ¯ x “ p x ` q{? νν