Wave Resolution Properties and Weighted Essentially Non-Oscillatory Limiter for Discontinuous Galerkin Methods by Xinghui Zhong B.S., University of Science and Technology of China, Hefei, China, 2007 M.S., Brown University, RI, 2009 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2012 c Copyright 2012 by Xinghui Zhong This dissertation by Xinghui Zhong is accepted in its present form by Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Chi-Wang Shu, Ph.D., Advisor Recommended to the Graduate Council Date Johnny Guzm´an, Ph.D., Reader Date Jan Hesthaven, Ph.D., Reader Approved by the Graduate Council Date Sheila Bonde, Dean of the Graduate School iii Vitae Biographical Information Born on August 7th, 1984 in Puyang, China Education Brown University, Providence, RI, USA Ph.D. in Applied Mathematics, May 2012 (expected) Advisor: Chi-Wang Shu Brown University, Providence, RI, USA M.S. in Applied Mathematics, May 2009 University of Science and Technology of China, Hefei, China B.S. with distinction (top 3%), in Mathematics and Applied Mathematics, 2007 Research Experience Brown University, Providence, RI, USA, 2009-present Research Assistant, Dissertation: Wave resolution properties and weighted ENO limiters for discontinu- ous Galerkin methods. Institute for Mathematics and its Applications (IMA), Minneapolis, MN, USA, August, 2009 Participated a 10-day workshop on Mathematical Modelling in Industry XIII. Conducted research on fast computational methods for reservoir flow models under the direction of Dr. Robert Shuttleworth from ExxonMobil Upstream Research. University of Science and Technology of China, Hefei, China, July-October 2005 Participated the “Undergraduate Research Project” program. Conducted research on boundary evaluation of interval vector operations under the direction of Prof. Jiansong Deng. iv Publications/Preprints • X. Zhong and C.-W. Shu, Numerical resolution of discontinuous Galerkin meth- ods for time dependent wave equations, Computer Methods in Applied Mechan- ics and Engineering, v200 (2011), pp. 2814-2827. • X. Zhong and C.-W. Shu, A simple weighted essentially nonoscillatory limiter for Runge-Kutta discontinuous Galerkin methods, submitted to Journal of Computational Physics. • J. Zhu, X. Zhong, C.-W. Shu and J.X. Qiu, Runge-Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes, submitted to SIAM Journal on Scientific Computing. • W. Guo, X. Zhong and J.-M. Qiu, Superconvergence of discontinuous Galerkin method: eigen-structures analysis based on Fourier approach, in preparation. Teaching Experience Brown University Teaching Assistant Methods of Applied Mathematics - II, Spring 2009 Methods of Applied Mathematics - I, Fall 2008 Math Resource Center Tutor, 2007-2011 Certificate III, Sheridan center for Teaching and Learning, 2011-2012 Certificate I, Sheridan Center for Teaching and Learning, 2009- 2010 Awards Brown University Dunmu Ji Award, 2012 AMS Graduate Student Travel Grant, January 2012 NSF Travel Support, June 2011 Graduate Student Fellowships in Fall 2007, Spring 2008, Fall 2011 University of Science and Technology of China Outstanding Undergraduate Research Project Award, 2005 Cyrus Tang Scholarship, 2004-2006 Outstanding Freshman Scholarship, 2003 v Acknowledgements First and foremost, I owe my deepest gratitude to my advisor Prof. Chi-Wang Shu for his patience, inspiration and enthusiasm. His wide knowlege has helped me in all the time of research and writing of this thesis whilst allowing me the room to work in my own way. His guidance has helped me to see life and science in their full depth. His detailed and constructive comments have had a remarkable influence on my entire study. One simply could not wish for a better or friendlier supervisor. I also wish to express my warm and sincere thanks to his wife Din-Sui Loh for he ever smiling face and encouraging words. It is an honor for me to thank the members of my committee, Prof. Jan Hes- thaven and Prof. Johnny Guzman for their detailed review, constructive criticism and excellent advice during the preparation of this thesis. I am grateful to the Division of Applied Mathematics and graduate school of Brown University for generous financial support in the past five years. I am indebted to many of my close friends with whom I have been able to share the joys during the five long years of my PhD studies. I cannot thank enough my best friend in high school (Xiaoping Cheng, Suhong Yang) and my best friends as an vi undergraduate (Liquan Huang, Han Lv) for all the emotional support, comraderie, entertainment, and caring they provided. Many friends at Brown University have helped me a lot, too. It is a pleasure to thank my former roommates and friends (Yanan Liu, Xiaomin Ma) and my present roommate Luan Lin, whose words and encouragements have always been very special to me. I am thankful to all my previ- ous and current group members for methodically looking over my work and helping me at every instance. Last but definitely not the least, I owe my loving thanks to my parents for giving birth to me at the first place and supporting me spiritually throughout my life. My special gratitude is due to my brothers for their unflagging love. I could never finish this dissertation without them. I finish with a final silence of gratitude for my life. Abstract of “Wave Resolution Properties and Weighted Essentially Non-Oscillatory Limiter for Discontinuous Galerkin Methods” by Xinghui Zhong, Ph.D., Brown Uni- versity, May 2012 This dissertation presents wave resolution properties and weighted essentially non- oscillatory limiter for discontinuous Galerkin methods solving hyperbolic conserva- tion laws. In this dissertation, using Fourier analysis, we provide a quantitative error anal- ysis for the semi-discrete DG method applied to time dependent linear convection equations with periodic boundary conditions. We apply the same technique to show that the error is of order k + 2 superconvergent at Radau points on each element and of order 2k + 1 superconvergent at the downwind point of each element, when using piecewise polynomials of degree k. An analysis of the fully discretized approxima- tion is also provided. We compute the number of points per wavelength required to obtain a fixed error for several fully discrete schemes. We also investigate a simple limiter using weighted essentially non-oscillatory (WENO) methodology for the Runge-Kutta discontinuous Galerkin (RKDG) meth- ods solving conservation laws, with the goal of obtaining a robust and high order limiting procedure to simultaneously achieve uniform high order accuracy and sharp, non-oscillatory shock transitions. The idea of this limiter is to reconstruct the en- tire polynomial, instead of reconstructing point values or moments in the classical WENO reconstructions. That is, the reconstruction polynomial on the target cell is a convex combination of polynomials on this cell and its neighboring cells and the nonlinear weights of the convex combination follow the classical WENO procedure. The main advantage of this limiter is its simplicity in implementation, especially for multi-dimensional meshes. Contents Vitae iv Acknowledgments vi 1 Introduction 1 2 Review of Discontinuous Galerkin Methods for Conservation Laws 4 2.1 Algorithm Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Time Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Error Estimates for Smooth Solutions . . . . . . . . . . . . . . . . . . 8 2.4 Minmod Limiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3 Numerical Resolution of Discontinuous Galerkin Methods for Time viii Dependent Wave Equations 15 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3 Implementation of DG Algorithm . . . . . . . . . . . . . . . . . . . . 21 3.4 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4.1 The Case of P 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4.2 The Case of P 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.4.3 The Case of P 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.5 Superconvergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.5.1 Superconvergence at Radau points . . . . . . . . . . . . . . . 37 3.5.2 Superconvergence at the Downwind Point . . . . . . . . . . . 40 3.6 Fully Discrete Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.6.1 The Case of P 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.6.2 The Case of P 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.6.3 The Case of P 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.6.4 A Measure of Work per Wavelength . . . . . . . . . . . . . . . 48 ix 3.7 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.8 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4 A Simple Weighted Essentially Non-Oscillatory Limiter for Runge- Kutta Discontinuous Galerkin Methods 62 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 DG Algorithm Formulation . . . . . . . . . . . . . . . . . . . . . . . 65 4.3 A New WENO Limiter . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.1 Identify the Troubled Cells . . . . . . . . . . . . . . . . . . . . 66 4.3.2 Reconstruction of the New Polynomials in the Troubled Cells Using a WENO Limiter: Scalar Case . . . . . . . . . . . . . . 68 4.3.3 WENO Limiting Procedure for Systems . . . . . . . . . . . . 77 4.4 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.4.1 Scalar Conservation Laws . . . . . . . . . . . . . . . . . . . . 82 4.4.2 Euler System in One Dimension . . . . . . . . . . . . . . . . . 88 4.4.3 Euler System in Two Dimension . . . . . . . . . . . . . . . . . 95 5 Conclusion 103 x A 106 A.1 Cardano’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 A.2 Ferrari’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 B 108 B.1 Lagrange Interpolation polynomial . . . . . . . . . . . . . . . . . . . 108 B.1.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 B.2 Lagrange interpolation polynomial . . . . . . . . . . . . . . . . . . . 109 B.2.1 Properties of Lagrange interpolation polynomial and Lagrange basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 B.2.2 Existence and uniqueness of Lagrange interpolation polynomials110 B.2.3 Error in Lagrange interpolation . . . . . . . . . . . . . . . . . 111 xi List of Tables 2.1 The CFL numbers for polynomials of degree k . . . . . . . . . . . . 13 2.2 ut + ux = 0 with u(x, 0) = sin x at t = 2π with RKDG method . . . . 14 3.1 Leading term of the lower bound of Mk as a function of q . . . . . . . 36 3.2 {ζk,l}kl=0 : the roots of Pk+1 (x) − Pk (x) . . . . . . . . . . . . . . . . . . 38 3.3 Lower bound of Mpf as a function of q . . . . . . . . . . . . . . . . . . 41 3.4 Leading term of the lower bound of Mk as a function of q . . . . . . . 48 3.5 When using P k polynomials on a uniform mesh of N cells, t = 1 . . . 51 3.6 When using P k polynomials on a uniform mesh of N cells, t = 100 . . 51 3.7 When using P k polynomials on a uniform mesh of N cells, t = 1000 . 52 3.8 Error at Radau points using P 1 DG and RK(9,9) on a uniform mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 xii 3.9 Error at Radau points using P 2 DG and RK(9,9) on a uniform mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.10 Error at Radau points using P 3 DG and RK(9,9) on a uniform mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.11 Error at Radau points using P 1 DG and RK(9,9) on a random mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.12 Error at Radau points using P 2 DG and RK(9,9) on a random mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.13 Error at Radau points using P 3 DG and RK(9,9) on a random mesh of N cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.14 Necessary number of points per wavelength for ek ≤ ε when time period q = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.15 Necessary number of points per wavelength for ek ≤ ε when time period q = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.16 When time period q = 1 . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.17 When time period q = 10 . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1 1D Transport equation. t = 2π. Uniform mesh with N cells . . . . . . 82 4.2 2D Transport equation. t = π. Uniform mesh with N × N cells. . . . 83 xiii 4.3 1D Burgers equation. t = 0.5. Uniform mesh with N cells. . . . . . . 84 4.4 1D Burgers equation. t = 0.5. Nonuniform mesh with N cells. . . . . 85 4.5 2D Burgers equation at t = 0.25. Uniform mesh with N × N cells. . . 86 4.6 2D Burgers equation at t = 0.25. Nonuniform mesh with N × N cells 86 4.7 1D Euler system with initial condition ρ(x, 0) = 1+0.2 sin x, v(x, 0) = 1, p(x, 0) = 1 at t = 2π. Uniform mesh with N cells. . . . . . . . . . . 89 4.8 1D Euler system with initial condition ρ(x, 0) = 1+0.2 sin x, v(x, 0) = 1, p(x, 0) = 1 at t = 2π. Nonuniform mesh with N cells. . . . . . . . . 90 4.9 2D Euler equation with initial condition ρ(x, y, 0) = 1 + 0.2 sin(x + y), u(x, y, 0) = 0.7, v(x, y0) = 0.3, p(x, y, 0) = 1 at t = 2π. Uniform mesh with N × N cells. . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.10 2D Euler equation with initial condition ρ(x, y, 0) = 1 + 0.2 sin(x + y), u(x, y, 0) = 0.7, v(x, y0) = 0.3, p(x, y, 0) = 1 at t = 2π. Nonuni- form mesh with N × N cells. . . . . . . . . . . . . . . . . . . . . . . . 96 4.11 2D Euler System of Smooth Vortex Evolution at t = 2. Uniform mesh with N × N cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.12 2D Euler System of Smooth Vortex Evolution at t = 2. Nonuniform mesh with N × N cells. . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.13 Percentage of troubled cells in the double Mach problem . . . . . . . 99 xiv List of Figures 3.1 Local basis functions of P 1 case. . . . . . . . . . . . . . . . . . . . . . 24 3.2 Fourier analysis procedure for P 1 case. . . . . . . . . . . . . . . . . . 27 3.3 Mesh points related to the local basis functions of P 2 case. . . . . . . 31 3.4 Mesh points related to the local basis functions of P 3 case. . . . . . . 34 3.5 Local basis function of P 1 based on Radau points. . . . . . . . . . . . 37 3.6 Time evolution of L∞ error (left) and e1 error (right) for P 1 DG solution with N = 40, cfl = 0.3. . . . . . . . . . . . . . . . . . . . . . 53 3.7 Time evolution of L∞ error (left) and e2 error (right) for P 2 DG solution with N = 80, cfl = 0.2. . . . . . . . . . . . . . . . . . . . . . 53 3.8 Time evolution of L∞ error (left) and e3 error (right) for P 3 DG solution with N = 160, cfl = 0.1. . . . . . . . . . . . . . . . . . . . . . 53 3.9 Necessary points per wave for kek∞ ≤ ε . . . . . . . . . . . . . . . . 59 xv 3.10 Comparison of time evolution of the errors with two different initial interpolations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1 DG polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.2 Linear weights for two dimensional rectangular mesh. . . . . . . . . . 77 4.3 Burgers equation at t = 1.5 with N = 80 cells. . . . . . . . . . . . . . 84 4.4 2D Burgers solution at t = 0.75 with 80 × 80 cells. . . . . . . . . . . . 87 4.5 2D Burgers solution that cuts along the diagonal at t = 0.75 with 80 × 80 cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.6 Buckley-Leverett problem at t = 0.4 with 80 cells. . . . . . . . . . . 88 4.7 Sod problem. t = 2. M = 0.01. N = 100. Density. . . . . . . . . . . . 91 4.8 Sod problem. t = 2. M = 30. N = 100. Density. . . . . . . . . . . . . 91 4.9 Lax problem. t = 1.3. M = 0.01. N = 100. Density. . . . . . . . . . . 91 4.10 Lax problem. t = 2. N = 100. Density. . . . . . . . . . . . . . . . . . 92 4.11 The shock density wave interaction problem. t = 1.8. M = 0.01. N = 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.12 The shock density wave interaction problem. t = 1.8. M = 300. N = 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 xvi 4.13 The blast wave problem. t = 0.038. M = 0.01. N = 400. . . . . . . . 94 4.14 The blast wave problem. t = 0.038. M = 200. N = 400. . . . . . . . 94 4.15 Double Mach reflection problem. M = 0.01. 960 × 240 cells. Thirty two equally spaced density contours from 1.35 to 23. Top: k = 1. Bottom: k = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.16 Double Mach reflection problem. M = 100. 960 × 240 cells. Thirty two equally spaced density contours from 1.35 to 23. Top: k = 1. Bottom: k = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.17 Double Mach reflection problem. 480 × 120 cells. M = 0.01 (top), M = 100 (middle) and M = 200 (bottom). Thirty two equally spaced density contours from 1.35 to 23. Left: k = 1. Right: k = 2. . . . . . 101 4.18 Double Mach reflection problem. 960 × 240 cells. M = 0.01 (top), M = 100 (middle) and M = 200 (bottom). Thirty two equally spaced density contours from 1.35 to 23. Left: k = 1. Right: k = 2. . . . . . 102 xvii Chapter 1 Introduction 1 2 Discontinuous Galerkin (DG) methods are a class of finite element methods us- ing completely discontinuous basis functions, which are usually chosen as piecewise polynomials. The first discontinuous Galerkin method was introduced by Reed and Hill [53] in 1973, to solve the neutron transport equation, i.e. a time independent linear hyperbolic equation. The type of DG methods we will discuss in this dissertation is the Runge-Kutta discontinuous Galerkin (RKDG) method [24, 23, 21, 20, 25], which is using explicit and nonlinearly stable high order Runge-Kutta method for time discretization and the DG method for space discretization. This method has several advantages, such as local conservation, the allowance of arbitrary triangu- lation, excellent parallel efficiency, the capability in h-p adaptivity and the ability to sharply capture discontinuities (especially contact discontinuities) in the solution. Furthermore, it has certain superconvergence properties. DG methods are known to provide good wave resolution properties, especially for long time simulation. The critical issue for wave phenomena is how many points are needed to resolve a wave, namely to make sure that the error is smaller than a given tolerance. In 1972, Kreiss and Oliger [42] introduced the error analysis of the wave problem using the finite difference method and developed the estimates of the neces- sary points per wavelength. 1974, Swartz and Wendroff [65] extended these ideas to the finite element method. In this dissertation, we combine and extend these ideas to DG methods. Using Fourier analysis, we develop a quantitative error estimate for DG methods and an estimate of the number of points per wavelength required to obtain a fixed error. The quantitative analysis provided in this dissertation can be used in the application to guide the choice of mesh size and polynomial degree in the DG method. Moreover, by investigating the quantitative error at Radau points using the same technique, we show that the errors at Radau points are of order k + 2 3 with the exception of the Radau point at the downwind end of the element, which is of order 2k + 1, when using the DG method with piecewise polynomials of degree k. DG methods have a lot of advantages, one of which is that they can compute so- lutions to hyperbolic conservation laws, which are either smooth or have weak shocks and other discontinuities, without further modification. However, if the discontinu- ities are strong, the polynomial representation will generate significant oscillations and even nonlinear instability. To avoid such difficulties, a slope limiter is bor- rowed from the finite volume methodology to control the numerical solution. Many such limiters exist in the literature, all with advantages and disadvantages. It is no exaggeration to state that the design of good, robust limiters remains one of the bottlenecks to the development of DG methods for solving conservation laws. In this dissertation, we use the WENO methodology to design a new and simpler limiter for the RKDG methods. The idea of this simple WENO limiter is that the reconstruction polynomial on the target cell is a convex combination of the polynomials on this cell and its immediate neighboring cells, with the nonlinear weights of the linear combination following the classical WENO procedure. Since this procedure uses information only from immediate neighbors and simple positive linear weights, it is significantly simpler to implement than previous WENO type limiters for RKDG methods. This dissertation is organized as follows: we first present an overview of DG methods in Chapter 2. Then we investigate the wave resolution property and su- perconvergence property of DG methods in Chapter 3. In Chapter 4, we design a new and simple WENO limiter for RKDG methods. Concluding remarks are given in Chapter 5. Chapter 2 Review of Discontinuous Galerkin Methods for Conservation Laws 4 5 The discontinuous Galerkin method was first designed as an effective numerical methods for solving hyperbolic conservation laws, which may have discontinuous solutions. In this chapter, we will review the algorithm formulation, error estimates and minmod type limiters for the DG method solving hyperbolic conservation laws. 2.1 Algorithm Formulation In this section, we give an overview of the algorithm formulation of the RKDG method solving one-dimensional conservation law ut + f (u)x = 0 (2.1) To define the DG method for (2.1), we consider a partition of the computational domain [a, b] in N cells as follows a = x 1 < x 3 < · · · < xN + 1 = b. 2 2 2 We denote 1  Ij = [xj− 1 , xj+ 1 ], xj = x 1 + xj− 1 2 2 2 j+ 2 2 as the cell and cell center, respectively. We again denote ∆xj = xj+ 1 − xj− 1 ; h = max ∆xj . 2 2 1≤j≤N We assume that the mesh is regular, namely there is a constant c > 0 independent of h such that ∆xj ≥ ch, 1 ≤ j ≤ N. 6 Define the approximation space as  Vhk = v : v|Ij ∈ P k (Ij ); 1 ≤ j ≤ N where P k (Ij ) denotes the set of polynomials of degree up to k defined on the cell Ij . The semi-discrete DG method for solving (2.1) is defined as follows: find the unique function uh = uh (t) ∈ Vhk such that, for j = 1, · · · , N, Z Z (uh )t vh dx − f (uh )(vh )x dx + fˆj+ 1 vh (x− j+ 1 ) − fˆj− 1 vh (x+ j− 1 )=0 (2.2) 2 2 2 2 Ij Ij holds for all test functions vh ∈ Vhk . fˆj+ 1 is the numerical flux, which is a single 2 valued function defined at the cell interface and in general depends on the values of the numerical solution uh from both sides of the interface   fˆj+ 1 = fˆ uh (x− j+ 1 , t), uh (x + j+ 1 , t) . 2 2 2 We use monotone numerical fluxes from finite difference and finite volume schemes for solving conservation laws, which satisfy the following conditions: • Consistency: fˆ(u, u) = f (u). • Continuity: fˆ(u− , u+ ) is at least Lipschitz continuous with respect to both arguments u+ and u− . • Monotonicity: fˆ(u− , u+ ) is a non-decreasing function of its first argument u− and a non-increasing function of its second argument u+ . Symbolically fˆ(↑, ↓). Here and below u− , u+ denote the left and right limits of the function u at the cell interface, respectively. 7 There is a class of monotone fluxes called upwind fluxes, which satisfies    f (u− ), if f ′ (u) ≥ 0, ∀u ∈ [min(u− , u+ ), max(u− , u+ )], ˆ − + f(u , u ) =   fˆ(u+ ), if f ′ (u) < 0, ∀u ∈ [min(u− , u+ ), max(u− , u+ )]. Well know monotone fluxes include the Lax-Friedrichs flux  1    fˆLF u− , u+ = f u− + f u+ − α u+ − u− , α = max |f ′ (u)| ; 2 u the Godunov flux    minu− ≤u≤u+ f (u),  if u− < u+ , fˆGod − + u ,u =   maxu+ ≤u≤u− f (u), if u− ≥ u+ ; and the Engquist-Osher flux Z u− Z u+ fˆEO = ′ max (f (u), 0) du + min (f ′ (u), 0) du + f (0). 0 0 Clearly, fˆGod and fˆEO are upwind fluxes. For more details about monotone fluxes, see, e.g., [45]. 2.2 Time Discretization The semi-discrete scheme (2.2) can be written as ut = L(u), 8 where L(u) is the spatial discretization operator. To discretize the time variable, we often use a class of high order nonlinearly stable Runge-Kutta methods, which are convex combinations of first order forward Euler steps. Therefore they maintain strong stability properties of the forward Euler step in any semi-norm. These methods were first developed in [62, 58], and later generalized in [31, 32]. The most popular scheme in this class is the following total variation diminishing (TVD) third order Runge-Kutta method [62]: u(1) = un + ∆tL(un ), 3 1 1 u(2) = un + u(1) + ∆tL(u(1) ), (2.3) 4 4 4 1 2 2 un+1 = un + u(2) + ∆tL(u(2) ). 3 3 3 Other TVD, or strong stability preserving (SSP) time discretizations can of course also be used. For details, see the survey paper [60] and the review paper [32]. 2.3 Error Estimates for Smooth Solutions If we assume the exact solution of (2.1) is smooth, we can obtain optimal L2 error estimates. The first error estimate of the DG method of Reed and Hill [53] was done by Lesaint and Raviart [44] in 1974. They proved that this method is convergent  with optimal order of accuracy, namely O hk+1 , in L2 norm, when piecewise tensor product polynomials of degree k are used as basis functions. Later in 1986, Johnson 1 and Pitk¨aranta [39] proved that L2 error estimate of this method is O hk+ 2 where k is again the polynomial degree and h is the mesh size, when the solution is sufficiently smooth for arbitrary meshes. 9 Such error estimates can be obtained for the general conservation law (2.1): Proposition 2.3.1. The solution uh of the DG scheme (2.2) for (2.1) with a suffi- ciently smooth solution u satisfies the following error estimate ku − uh k ≤ Chk+1 , if fˆ is an upwind flux; (2.4) 1 ku − uh k ≤ Chk+ 2 , if fˆ is general monotone flux, (2.5) where C depends on u and its derivatives but is independent of h. For the proof of Proposition 2.3.1, see [71, 72]. In actual numerical computations, one almost always observe the optimal O(hk+1) accuracy. 2.4 Minmod Limiters In practice, especially for problems containing strong discontinuities, the DG scheme is significant oscillatory near discontinuities and even nonlinear instable. We often need to apply nonlinear limiters to control these oscillations. In this section, we discuss the minmod type limiter [24, 23, 21, 20, 25]. For simplicity, we consider the forward Euler time discretization of the semi- discrete scheme (2.2). Starting from a solution unh ∈ Vhk at time level n (for the initial condition, u0h is taken as the L2 projection of the analytical initial condition into Vhk ). We would like to “limit” it to obtain a new function un,new h before advancing it to next time level. That is: to find un+1 h ∈ Vhk , such that, for j = 1, · · · , N, Z Z un+1 − un,new h h vh dx − f (un,new h n,new )(vh )x dx + fˆj+ 1 vh (x− j+ 1 n,new ) − fˆj− 1 vh (x+ j− 1 )=0 Ij ∆t Ij 2 2 2 2 10 holds for all test functions vh ∈ Vhk . The limiting procedure to go from un to un,new should satisfy the following con- ditions: • The cell averages of un,new h and unh are the same. That is to ensure the conser- vation. • In the smooth regions this limiter does not change the solution, un,new h (x) = unh (x). There are many limiters discussed in the literature. We provide here an overview of the minmod limiter [23] in the one-dimensional scalar case. Denote the cell average of the solution u as Z 1 u¯j = uh dx (2.6) ∆xj Ij and further denote u˜j = u− j+ 1 − u¯j , u˜˜j = u¯j − u+ j− 1 . (2.7) 2 2 u˜j and u˜˜j are modified by the minmod limiter [33], (mod) (mod)  u˜j = m (˜ uj , ∆+ u¯j , ∆− u¯j ) , u˜˜j = m u˜˜j , ∆+ u¯j , ∆− u¯j , (2.8) where ∆+ u¯j = u¯j+1 − u¯j , ∆− u¯j = u¯j − u¯j−1 , 11 and the minmod function m is defined by    s min1≤i≤l |ai | if s = sign(a1 ) = · · · = sign(al ), m(a1 , · · · , al ) = (2.9)   0, otherwise. Then the new point values at the cell boundary are     (mod) (mod) (mod) (mod) uh − xj+ 1 = u¯j + u˜j , uh xj− 1 = u¯j − u˜˜j + (2.10) 2 2 (mod) by definition (2.7). The limited solution uh is then reconstructed to maintain the old cell average (2.6) and the new point values (2.10). Clearly, this reconstruction (mod) is unique for k ≤ 2. For k > 2, we have extra freedom in obtaining uh . It has been proved that the DG scheme with this “pre-processing” limiter is TVD by using Harten’s Lemma [33]. We now need to verify that this limiter does not affect accuracy in smooth regions. After a simply Taylor expansion, we have 1 1 u˜j = ux (xj )∆xj + O(h2 ), u˜˜j = ux (xj )∆xj + O(h2), 2 2 1 1 ∆+ u¯j = ux (xj )(∆xj + ∆xj+1 ) + O(h2 ), ∆− u¯j = ux (xj )(∆xj + ∆xj−1 ) + O(h2 ). 2 2 Clearly, if we are in a smooth and monotone domain, i.e. ux (xj ) is away from zero, the first argument in the minmod function (2.8) is of the same sigh as the second and third argument and is smaller in magnitude, when h is small. Therefore, the (mod) (mod) modified values u˜j and u˜˜j in (2.8) will take the unmodified values u˜j and u˜˜j , respectively. That is, this limiter does not affect accuracy in smooth monotone regions. However, from numerical results, this TVD limiter does kill accuracy at smooth extrema. Therefore, in practice, we often use a total variation bounded 12 (TVB) modified minmod limiter [57]    a1 if |a1 | ≤ Mh2 , m(a ˜ 1 , · · · , al ) = (2.11)   m(a1 , · · · , al ), otherwise, where the TVB parameter M depends on the solution of the problem. For scalar problems it is possible to estimate M by the initial condition (M is proportional to the second derivative of the initial condition at smooth extrema); however it is more difficult to estimate M for the system case. If M is chosen too small, accuracy may degenerate at smooth extrema of the solution; however, if M is chosen too large, oscillations will appear. For more details, see [23]. 2.5 Numerical Examples In this section, we demonstrate the performance of the DG method, using the linear version of (2.1): ut + ux = 0, 0 ≤ x ≤ 2π (2.12) with the initial condition u(x, 0) = sin x and periodic boundary conditions. The exact solution of (2.12) is u(x, t) = sin(x − t). ˆ − , u+ ) = u− . We use The monotone flux is taken as the simple upwind flux f(u the 3rd order RK method discussed in Section 2.1 for time discretization and DG methods with polynomials of degree k for space discretization. Table 2.1 lists the CFL number, which are chosen according to the CFL table in [26]. For the P k case 13 with k ≥ 3, ∆t is further reduced in the accuracy test. Table 2.2 lists the errors and Table 2.1: The CFL numbers for polynomials of degree k k 0 1 2 3 4 5 6 RK3 1.20 0.40 0.20 0.14 0.08 0.06 0.05 orders of the accuracy for the RKDG method with P k polynomials, for 0 ≤ k ≤ 6. We observe that the error is of order k + 1. Note here we collect L2 , L1 , L∞ errors. For the error tables in the following chapters, we only collect L1 and L∞ errors, since Z 2π 1 kf kL1 = |f (x)|dx 2π 0 s s Z 2π Z 2π 1 2 ≤ |f (x)| dx 1dx 2π 0 0 s Z 2π 1 = |f (x)|2 dx 2π 0 = kf kL2 ≤ kf kL∞ . For more numerical examples about DG methods and minmod limiters, see, e.g. [24, 23, 21, 20, 25]. 14 Table 2.2: ut + ux = 0 with u(x, 0) = sin x at t = 2π with RKDG method N L2 error order L1 error order L∞ error order 20 4.50E-01 4.05E-01 6.43E-01 40 2.78E-01 0.70 2.50E-01 0.70 3.98E-01 0.69 P 0 80 1.55E-01 0.84 1.40E-01 0.84 2.22E-01 0.84 160 8.25E-02 0.91 7.42E-02 0.91 1.18E-01 0.92 320 4.25E-02 0.96 3.83E-02 0.96 6.05E-02 0.96 10 2.41E-02 1.99E-02 4.44E-02 20 4.56E-03 2.40 3.76E-03 2.40 9.84E-03 2.17 P 1 40 1.02E-03 2.16 8.20E-04 2.20 2.90E-03 1.76 80 2.46E-04 2.05 1.91E-04 2.10 7.79E-04 1.90 160 6.07E-05 2.02 4.61E-05 2.05 2.01E-04 1.95 10 8.63E-04 7.28E-04 2.33E-03 20 1.05E-04 3.04 8.70E-05 3.06 3.05E-04 2.94 P 2 40 1.30E-05 3.01 1.07E-05 3.02 3.86E-05 2.98 80 1.62E-06 3.00 1.33E-06 3.01 4.87E-06 2.99 160 2.02E-07 3.00 1.66E-07 3.00 6.10E-07 3.00 10 8.12E-05 6.90E-05 1.71E-04 20 5.08E-06 4.00 4.31E-06 4.00 1.08E-05 3.98 3 P 40 3.17E-07 4.00 2.69E-07 4.00 6.82E-07 3.99 80 1.98E-08 4.00 1.68E-08 4.00 4.27E-08 4.00 160 1.24E-09 4.00 1.05E-09 4.00 2.67E-09 4.00 10 9.27E-06 8.33E-06 1.32E-05 20 2.91E-07 5.00 2.61E-07 4.99 4.16E-07 4.99 P4 40 9.09E-09 5.00 8.19E-09 5.00 1.29E-08 5.01 80 2.85E-10 5.00 2.56E-10 5.00 4.05E-10 4.99 10 2.46E-06 2.21E-06 3.51E-06 20 3.85E-08 6.00 3.46E-08 6.00 5.55E-08 5.98 P 5 30 3.38E-09 6.00 3.04E-09 6.00 4.86E-09 6.01 40 6.01E-10 6.00 5.41E-10 6.00 8.68E-10 5.99 50 1.58E-10 5.99 1.42E-10 5.99 2.28E-10 6.00 10 8.93E-07 8.04E-07 1.26E-06 20 6.99E-09 7.00 6.29E-09 7.00 9.88E-09 7.00 P6 30 4.09E-10 7.00 3.68E-10 7.00 5.79E-10 7.00 40 5.50E-11 6.98 4.95E-11 6.98 7.77E-11 6.98 Chapter 3 Numerical Resolution of Discontinuous Galerkin Methods for Time Dependent Wave Equations 15 16 3.1 Introduction In this chapter, we consider the following time dependent linear wave problem ut + aux = 0, x ∈ [0, 2π], t > 0 (3.1) iωx u(x, 0) = e , x ∈ [0, 2π] with periodic boundary conditions, where a is the phase speed (for simplicity, assume √ that a > 0), i = −1 and ω is the wave number (for convenience, assume that ω > 0). We develop a quantitative error analysis for the discontinuous Galerkin (DG) solution via Fourier analysis and study the superconvergence property of the DG solution. Because the equation (3.1) as well as the DG scheme are linear, a general L2 initial condition can be written as a sum of the simple waves eiωx (Fourier series) and the numerical solution for such general initial condition is just a superposition of the numerical solution of the equation (3.1) with the simple wave initial condition. The results of this chapter is therefore applicable to such general initial conditions when we have a clear understanding of the wave numbers we would like to resolve. The quantitative analysis provided in this chapter can be used to guide the choice of mesh size and polynomial degree in the DG method, when we are solving a linear wave equation, and would like to resolve the first k waves, to a given threshold of error, up to a certain time t. The results of this chapter is also valid for one dimensional linear hyperbolic system, and can provide useful guidelines when the DG method is used to solve more complicated linear and nonlinear wave equations. The superconvergence behavior of the classical finite element method has been analyzed for many years. For example, in [27], Douglas and Dupont showed that for a general class of two-point boundary value problems, the rate of convergence at the 17 mesh points is of order 2k when polynomials of degree k are used. In [10], Bakker proved that the C 0 Galerkin solution of a two-point boundary value problem using piecewise P k polynomials, has superconvergence of order 2k at the knots and of order k + 2 at the Lobatto points of each segment, and the gradient is superconvergent of order k + 1 at the zeros of a Legendre polynomial shifted to the elements of the partition. In [11], the same author proved that for two classes of Galerkin methods: the Ritz-Galerkin method for 2m-th order self-adjoint boundary value problems and the collocation method for arbitrary m-th order boundary value problems, the solu- tion exhibits superconvergence of order k + 2 at the zeros of the Jacobi polynomial Pnm,m (σ) shifted to the elements of the partition, and the derivative of the solution m−1,m−1 is superconvergent of order k + 1 at the zeros of the Jacobi polynomial Pn+1 (σ) shifted to the elements, where n = k +1−2m and k is the degree of the finite element space. Based on the results of Adjerid et al. [1] for singularly-perturbed parabolic sys- tems, Biswas et al. [13] used the assumption that the DG solutions of hyperbolic con- servation laws using k-degree polynomial approximation exhibit superconvergence at the roots of Radau polynomial of degree k + 1 (Radau points), to construct a poste- riori estimate of spatial discretization error. Adjerid et al. [2] proved that the DG solution of the ordinary differential equa- tion (ODE) u′ − f (u) = 0 is superconvergent of order 2k + 1 at the downwind end of each element while maintaining an order of k + 2 at the remaining Radau points. Numerical examples for the partial differential equations (PDEs) were also shown without analysis. Later, these results were extended to two-dimensional problems on rectangular meshes [4] and nonlinear hyperbolic problems [5]. Castillo [15] inves- tigated the existence of superconvergent points for DG methods applied to elliptic problems and showed that on each element the k-degree local discontinuous Galerkin 18 (LDG) solution gradient is superconvergent of order k + 1 at the shifted roots of the k-degree Legendre polynomial. For these cases, the model problems are independent of time. In [3], Adjerid et al. showed that the LDG solutions of convection-dominated problems are superconvergent of order k + 2 at the shifted Radau points on each element. For diffusion-dominated problems, the derivative of the LDG solution is superconvergent of order k + 2 at Radau points. Later, Adjerid et al. proved that the DG solution is superconvergent of order k + 2 at Radau points for linear sym- metric hyperbolic systems [7] and for linear symmetrizable hyperbolic systems [6], when t = O(1). In this chapter, by investigating the quantitative error at Radau points using Fourier analysis, we show that, even for t greater than O(1), the errors at Radau points are of order k + 2 with the exception of the Radau point at the downwind end of the element, which is of order 2k + 1, when using DG method with piecewise polynomials of degree k. See, for example [40, 16, 17, 18, 19, 22] for more superconvergence results of the DG method. 3.2 Motivation This work was motivated by [42], where the error analysis is given using the finite difference method. We first provide an overview of this method solving the linear wave problem (3.1). The exact solution of (3.1) is U(x, t) = eiω(x−at) . (3.2) 2π Assume a uniform grid with ∆x = N . Using the following second order central 19 finite difference method to approximate the spatial derivative in (3.1) uj+1 − uj−1 D0 u = , 2∆x yields a semi-discrete version of (3.1) with a system of differential-difference equations given by duj (t) uj+1(t) − uj−1 (t) +a = 0, dt 2∆x (3.3) iωxj uj (0) = e . It is easy to compute the solution of (3.3), which is uj (t) = eiω(xj −ct) with the numerical phase speed a sin ξ c= , ξ = ω∆x. ξ By Taylor expansion, the leading term in the error between the exact solution U(x, t) and the approximate solution u(x, t) is ωatξ 2 e = kU − uk∞ := max |U(xj , t) − uj (t)| ≃ . (3.4) 0≤j≤N 6 It is evident that the error is a function of time. Under the assumption of periodic- ity, the important quantity is not the time elapsed, but rather the number of time periods, which is denoted by q, that is, t q= , t = ωat. (3.5) 2π 20 Now the critical issue is: how many grid points are needed to resolve a wave, to make sure that the error is smaller than a given tolerance ε? We introduce the number of points per wavelength, M, which is given by N 2π M= = . (3.6) ω ξ Note that M has a theoretical minimum of 2, since it takes a minimum of two points per wavelength to uniquely specify a wave. Rewriting the error (3.4) in terms of M and q, yields  2 πq 2π e≃ . (3.7) 3 M Then, the lower bound on M  π 1/2  q 1/2 M ≥ 2π , 3 ε required to ensure a specific error ε can be obtained from (3.7). According to the analysis above, assume that, for a given problem, we know the maximum wave number ω needed to adequately resolve the solution. Then we can choose the number of points to ensure that the wave is well resolved. This can reduce the cost of computation while preserving the accuracy. If more detailed information is known about the spectral distribution of the Fourier coefficients, then this can be used to obtain sharper comparisons of efficiency by weighting the error function appropriately. In 1974, Swartz and Wendroff [65] extended the ideas of [42] to the finite element method, using smooth splines as basis functions and presented an analysis of the 21 fully discretized approximation. Our approach is to combine and extend the ideas of [42] and [65] to the DG method, using the techniques introduced in [69] (see also [68, 70]) where the explicit formula of the DG solution is given by Fourier analysis. See, for example [55, 36, 37, 8, 9] for recent development and studies of the dispersion and dissipation errors of DG methods using Fourier analysis. 3.3 Implementation of DG Algorithm In this section, we present the implementation of the DG method solving the model problem (3.8): ut + aux = 0, x ∈ [0, 2π], t > 0 (3.8) u(x, 0) = sin ωx, x ∈ [0, 2π] with periodic boundary conditions. The exact solution to (3.8) is U(x, t) = sin(ω(x − at)). (3.9) As discussed in Section 2.1, to define the DG method for the model problem, we first consider a uniform partition of the computational domain [0, 2π] in N cells of size ∆x = 2πN . Denote the cell by Ij = [xj− 1 , xj+ 1 ] and the cell center by xj = 2 2   1 2 xj+ 1 + xj− 1 , j = 1, · · · , N, where 2 2 0 = x 1 < x 3 < · · · < xN + 1 = 2π 2 2 2 22 and denote the approximation space as  Vhk = v : v|Ij ∈ P k (Ij ); 1 ≤ j ≤ N (3.10) where P k (Ij ) denotes the set of polynomials of degree up to k defined on the cell Ij . The semi-discrete DG method using the upwind flux for solving (3.8) is defined as follows: find the unique function u = u(t) ∈ Vhk such that, for j = 1, · · · , N, Z Z ut v dx − a uvx dx + au− j+ 1 v(x− j+ 1 ) − au− j− 1 v(x+ j− 1 )=0 (3.11) Ij Ij 2 2 2 2 holds for all test functions v ∈ Vhk . Here and below u+ , u− denote the left and right limits of the function u at the cell interface, respectively. We now look at the implementation of the scheme (3.11). If a local basis of P k (Ij ) is chosen and denoted as φlj (x) for l = 1, 2, · · · , k + 1, then the numerical solution can be represented as k+1 X u(x) = ulj φlj (x), x ∈ Ij . (3.12) l=1 After substituting (3.12) into (3.11), we have k+1 Z ! X φlj (x)φm l j (x)dx (uj )t l=1 Ij | {z } Mm,l         !    Z      X  k+1  l m l l m l l m l = a φj (x)(φj (x))x dx uj + φj φj x− uj − φj−1φj x+ uj−1 ,   Ij 1 j+ 2 1 j− 2   l=1  | {z } | {z }    l   A u l 1m,l j B u 1m,l j−1 23 which can be written as the following matrix form duj T M = aA1 uj + aB1 uj−1 uj = u1j , · · · , uk+1 j . dt After inverting a local (k + 1) × (k + 1) mass matrix M , the DG scheme (3.11) can be written as duj a = (Auj + Buj−1 ) , (3.13) dt ∆x T where uj = u1j , · · · , uk+1 j , A and B are (k + 1) × (k + 1) constant matrices. As in [69], the local basis of P k (Ij ) is chosen to be the Lagrangian polynomials, based on the following k + 1 equally spaced points   2l − k xj+ 2l−k = xj + ∆x, l = 0, · · · , k. 2(k+1) 2(k + 1) In this way, uj , the coefficients of the solution u inside the cell Ij is a vector of length k + 1 containing the values of the solution at these points. The DG finite element scheme (3.13) becomes a finite difference scheme on a globally uniform mesh (with a ∆x mesh size k+1 ). However, it is not a standard finite difference scheme because each point in the group of k + 1 points belonging to the cell Ij obeys a different form. 3.4 Error Analysis In this section, we present the details of the DG algorithm formulation with the basis functions discussed in Section 3.3 and derive the error estimation of the semi-discrete scheme using piecewise P k polynomials, with k = 1, 2, 3. At the end of this section, 24 We provide an estimate of the necessary number of points per wavelength required to obtain a fixed error. 3.4.1 The Case of P 1 In this subsection, we consider the piecewise linear case, i.e. k = 1. We present the details of the error analysis for this case. φj− 41 φj+ 14 @ @ @ @ c c c @c c c xj− 1 xj+@1 4 4 xj− 23 xj− 21 xj+ 21 xj+ 23 Figure 3.1: Local basis functions of P 1 case. In this case, the local basis functions inside cell Ij are φj− 1 (x), φj+ 1 (x), which 4 4 are Lagrangian polynomials based on the points xj− 1 , xj+ 1 (Figure 3.1). With these 4 4 basis functions, the solution inside the cell Ij is then represented by u(x) = uj− 1 φj− 1 (x) + uj+ 1 φj+ 1 (x). (3.14) 4 4 4 4 After substituting (3.14) into (3.11), we get the finite difference representation of the DG method a   u′j− 1 = −5uj− 5 + 15uj− 3 − 7uj− 1 + 3uj+ 1 4 4∆x 4 4 4 4 (3.15) a   u′j+ 1 = uj− 5 − 3uj− 3 + 11uj− 1 − 9uj+ 1 . 4 4∆x 4 4 4 4 25 The scheme (3.15) can be rewritten into the matrix form duj a = (Auj + Buj−1 ) , (3.16) dt ∆x where       uj− 1 1 −7 −3 1 −5 15  4    uj =  , A=  , B=  . (3.17) 4 4 uj+ 1 11 −9 1 −3 4 To solve (3.16), the standard Fourier analysis is used here. This analysis depends heavily on the assumption of uniform mesh and periodic boundary conditions. As- sume      uj− 41 (t)   uˆ− 41 (t)  iωx  =  e j, (3.18)     uj+ (t) 1 uˆ (t) 1 4 4 and after substituting (3.18) into the DG scheme (3.16), the coefficient vector satisfies the following ODE system      uˆ′− 1 (t)  uˆ− 14 (t)   4  = G , (3.19)     uˆ′1 (t) uˆ 1 (t) 4 4 where G is the amplification matrix, given by a G= (A + Be−iξ ), ξ = ω∆x. (3.20) ∆x The two eigenvalues of G are a  −iξ p  λ1,2 = −e − 2 ∓ e−2iξ + 10e−iξ − 2 , (3.21) ∆x 26 and the corresponding eigenvectors are  √  iξ iξ 2iξ −1 + e ∓ 4 1 + 10e − 2e   V1,2 = . (3.22) iξ 1 + 11e Then the general solution of the ODE system (3.19) is    uˆ− 41 (t)    = C11 eλ1 t V1 + C12 eλ2 t V2 . (3.23)   uˆ 1 (t) 4 The coefficients C11 and C12 in (3.23) can be determined by the initial condition     − i4ξ  uˆ− 41 (0)   e   =     iξ  . (3.24) uˆ 1 (0) e4 4 We thus have the explicit solution of the DG scheme with piecewise linear polyno- mials for solving (3.8). Comparing this DG approximate solution with the exact solution (3.9) will give us the quantitative error. The error analysis procedure can be summarized in Figure 3.2. Now let us first introduce some notations: ke− 1 k∞ := max U(xj− 1 , t) − uj− 1 (t) , 4 1≤j≤N 4 4 (3.25) ke+ 1 k∞ := max U(xj+ 1 , t) − uj+ 1 (t) . 4 1≤j≤N 4 4 27 duj a = Auj + Buj−1 dt ∆x ˆ iωxj uj = ue ? ˆ ′ = Gu u ˆ λ, V ? ˆ = C1 eλ1 t V1 + C2 eλ2 t V2 u - uj 6 U(x) = sin(ω(x − at)) I.C. ? ? C1 , C2 error estimate Figure 3.2: Fourier analysis procedure for P 1 case. 28 From the error analysis procedure, we have ke+ 1 k∞ = max U(xj+ 1 , t) − uj+ 1 (t) 4 1≤j≤N 4 4   n o = max sin iω(xj+ 1 − at) − Im eiωxj uˆ+ 1 (t) 1≤j≤N 4 4 | {z } Θ             ∆x  = max sin Θ − Im eiΘ exp iω at − uˆ+ 1 (t) 1≤j≤N   4 4    | {z }   Taylor expansion    1 2 9i + 8ωt 3 = max sin Θ − Im eiΘ 1− ξ − ξ +··· 1≤j≤N 24 576           1 2 9i + 8ωt 3  iΘ = max sin Θ − Im e 1− ξ − ξ +··· 1≤j≤N   | 24 576     {z }  1+C3 +iC4 = max sin Θ − ((1 + C3 ) sin Θ + C4 cos Θ) 1≤j≤N = max |−C3 sin Θ + C4 cos Θ| 1≤j≤N q = C32 + C42 . Notice that (Taylor expansion part):    q ∆x exp iω at − uˆ+ 1 (t) − 1 = |1 + C3 + iC4 − 1| = C32 + C42 . 4 4 Therefore    ∆x ke+ 1 k∞ = exp iω at − uˆ+ 1 (t) − 1 4 4 4 29 By a simple Taylor expansion, we get 1 2 t 211 4 1861 t 5 ke+ 1 k∞ = ξ + ξ3 − ξ − ξ 4 24 72 27648 414720 2 ! 2808495 + 1364224 t + ξ6 1592524800 3 ! 11(−2854295 t + 965888 t ) − ξ7 33443020800 ! (3.26) 2 4 −6386427605 + 581426944 t + 1359970304 t + ξ8 12842119987200 3 5 ! 48576117895 t − 7643182848 t + 6799851520 t − ξ9 192631799808000 + O(ξ 10) where t = ωat. 1 2 From (3.26), for short time t, the first term 24 ξ , which is independent of t,   dominates. As t increases to O 1ξ , the second term 72t ξ 3 begins to dominate. If t continues increasing, another term with higher degree of t as coefficient may domi- nate. This means, during different time intervals, the dominant terms are different.   At t = O 1ξ , the error can be approximated by the first two terms, i.e. 1 2 1 ke− 1 k∞ ≃ ξ − tξ , 3 4 24 72 (3.27) 1 1 3 ke+ 1 k∞ ≃ ξ2 + tξ . 4 24 72  Denote e1 = max ke− 1 k∞ , ke+ 1 k∞ , where the subscript 1 is used here to indicate 4 4 30 the case k = 1. Clearly, the error e1 satisfies  1 1 3 e1 = max ke− 1 k∞ , ke+ 1 k∞ ≃ ξ 2 + tξ . (3.28) 4 4 24 72 Similar to the finite difference case, we would like to rewrite the error in terms of t the number of points per wavelength M1 and the number of time periods q = 2π . Notice that two points {xj− 1 , xj+ 1 } are used for each cell, therefore, 4 4 2N 4π M1 = = . (3.29) ω ξ Substituting (3.29) and (3.5) into (3.28), we get  2  3 1 4π πq 4π e1 ≃ + . (3.30) 24 M1 36 M1 Then the number of points per wavelength required to guarantee the error e1 ≤ ε satisfies the following inequality  2  3 1 4π πq 4π + ≤ ε. (3.31) 24 M1 36 M1 3.4.2 The Case of P 2 In this subsection, we use the same procedure as we did in the piecewise linear case to present the details of the error estimation of the DG method using piecewise P 2 polynomials.  In this case, the local basis functions φj− 1 (x), φj (x), φj+ 1 (x) inside cell Ij are 3 3 the Lagrangian polynomials based on the points xj− 1 , xj , xj+ 1 (Figure 3.3). The 3 3 31 c c c xj− 1 xj xj+ 1 3 3 xj− 21 xj+ 21 Figure 3.3: Mesh points related to the local basis functions of P 2 case. solution inside the cell Ij is then represented by u(x) = uj− 1 φj− 1 (x) + uj φj (x) + uj+ 1 φj+ 1 (x). 3 3 3 3 Denote  T uj = uj− 1 uj uj+ 1 , 3 3 then the DG scheme can be written into the matrix form (3.13) with  43 29 1   23 115 115  − − −  16 24 16   16 24 16               69 15 15   9 15 45  A= − − , B= − − . (3.32)  16 8 16   16 8 16              19 139 71 1 5 5 − − − − 16 24 16 16 24 16 The standard Fourier analysis is again used here. First we make an ansatz of the form     uj− 1 (t) uˆ− 1 (t)  3   3          iωx  u (t)  =  uˆ (t)  e j. (3.33)  j   0          uj+ 1 (t) uˆ 1 (t) 3 3 After substituting (3.33) into the DG scheme (3.13) with (3.32), the coefficient vector 32 satisfies     uˆ′− 1 (t) uˆ 1 (t)  3   −3           uˆ′ (t)  = G  uˆ (t) , (3.34)  0   0          uˆ′1 (t) uˆ 1 (t) 3 3 where G is given by (3.20) with A and B defined in (3.32). Then with three eigen- values of G and the corresponding eigenvectors computed by Mathematica (Here we use Cardano’s method to obtain the eigenvalues, see Appendix A.1), the general solution of the ODE system (3.34) is   uˆ 1 (t)  −3       uˆ (t)  = C21 eλ1 t V1 + C22 eλ2 t V2 + C23 eλ3 t V3 . (3.35)  0      uˆ 1 (t) 3 The coefficients C21 , C22 and C23 in (3.35) can be determined by the initial condition     − iξ uˆ 1 (0) e 3  −3             uˆ (t) = 1 . (3.36)  0            iξ uˆ 1 (0) e3 3 We thus have the explicit solution of the DG scheme with P 2 polynomials for solving (3.8). Comparing this with the exact solution (3.9) will give us the quantitative error 33 estimates. ke− 1 k∞ := max U (xj− 1 , t) − uj− 1 (t) 3 1≤j≤N 3 3 ! 2 ξ3 3223ξ 5 241 tξ 6 110175241 t = + − + − + ξ 7 + O(ξ 8 ), 1296 3110400 1008000 134369280000 80000 ke0 k∞ := max |U (xj , t) − uj (t)| 1≤j≤N ! 2 ξ3 4979ξ 5 11 tξ 6 28245223 t = − − + + ξ 7 + O(ξ 8 ), 240 15552000 1814400 2015539200000 432000 ke+ 1 k∞ := max U (xj+ 1 , t) − uj+ 1 (t) 3 1≤j≤N 3 3 ! 2 23ξ 3 67567ξ 5 263 tξ 6 451421052697 t = + + − − ξ 7 + O(ξ 8 ). 6480 357696000 4636800 8174355148800000 368000 Clearly, the coefficients of high order terms depend on t. Similarly, during different n time intervals, the dominant terms are different. Denote e2 := max ke− 1 k∞ , ke0 k∞ , 3 o   ke+ 1 k∞ . At t = O 1ξ , the error e2 satisfies 3 ξ3 e2 ≃ . (3.37) 240   1 At t = O ξ2 , the error e2 satisfies 2 ! 3 5 6 ξ 4979ξ 11 tξ 28245223 t e2 ≃ − − + + ξ 7. (3.38) 240 15552000 1814400 2015539200000 432000 In this case, the number of points per wavelength M2 satisfies 3N 6π M2 = = , (3.39) ω ξ since three points are used in each cell. Similarly, by rewriting (3.37) and (3.38) in t terms of M2 and q = 2π , we can obtain the number of points per wavelength M2 34 required to ensure that e2 ≤ ε. 3.4.3 The Case of P 3 In this subsection, we apply the same procedure discussed in the previous subsec- tions to the DG method using piecewise P 3 polynomials. The local basis functions c c c c xi− 3 xi− 1 xi+ 1 xi+ 3 8 8 8 8 xi− 21 xi+ 21 Figure 3.4: Mesh points related to the local basis functions of P 3 case. n o φj− 3 (x), φj− 1 (x), φj+ 1 (x), φj+ 3 (x) inside cell Ij are the Lagrangian polynomials 8 8 8 8 based on the points xj− 3 , xj− 1 , xj+ 1 , xj+ 3 (Figure 3.4). With these basis functions, 8 8 8 8 the DG scheme can be written in the matrix form (3.13) with   15109 4521 255 403  − 6144 − 2048 2048 6144         43577 7699 1115 959   − − −   6144 2048 2048 6144    A=     11761 10747 5629 7097   − − −   6144 2048 2048 6144         4723 7983 21993 43211  − − − 6144 2048 2048 6144 35 and   2865 12033 20055 20055  − 2048 2048 − 2048 2048         1685 7077 11795 11795     2048 − 2048 2048 − 2048    B=  .   2555   − 365 1533 − 2555   2048 2048 2048 2048         615 2583 4305 4305  − − 2048 2048 2048 2048 The standard Fourier analysis is once again used here to obtain the explicit solution of the DG scheme (Here we use Ferrari’s method to solve the eigenvalues of the ampli- fication matrix, see Appendix A.2). The error is also approximated by the dominant n o terms of Taylor series. Denote e3 := max ke− 3 k∞ , ke− 1 k∞ , ke+ 1 k∞ , ke+ 3 k∞ . At 8 8 8 8   1 t = O ξ2 , e3 ≃ 0.00036369202628968253968253968254ξ 4. (3.40)   At t = O ξ13 , e3 ≃0.00036369202628968253968253968254ξ 4 − 7.25319954423243464358961570099621 × 10−6 ξ 6 (3.41) + 7.08616780045351473922902494331066 × 10−7 tξ 7 . In this case, the number of points per wavelength M3 satisfies 4N 8π M3 = = , (3.42) ω ξ 36 since four points are used in each cell. Similarly, the number of points per wavelength M3 required to ensure that e3 ≤ ε can be obtained by rewriting (3.40) and (3.41) in terms of M3 and q. Table 3.1 shows the lower bound of Mk to ensure the specific error ε, which is a function of q. we only present the leading term here. The leading term is computed using asymptotic analysis. Take (3.31) and ε = 0.1 as an example:  2  3 1 4π πq 4π + ≤ 0.1. (3.43) 24 M1 36 M1 We look for a solution of the form M1 = Cq α . Using this expansion in the equation, we have 2π 2 α 16π 4 −0.1C 3 q 3α + Cq + q = 0. (3.44) 3 9 Then equating coefficients of q with highest degree to zero, we get 16π 4 (1) 3α = α, then α = 0. However, 9 6= 0. (2) α = 1, then −0.1C 3 should be 0. C = 0. 16π 4 (3) 3α = 1, then α = 31 . And −0.1C 3 = 9 ⇒ C ≈ 12. 1 Therefore, the leading term of M1 is 12q 3 . Table 3.1: Leading term of the lower bound of Mk as a function of q M 1   M2    M3   1 1 1 t O ξ O ξ O ξ2 O ξ12 O ξ13 1 2 1 ε = 0.1 12q 3 6.53 7q 7 6.17 6q 7 1 2 1 ε = 0.01 26q 3 14.08 10q 7 10.98 8q 7 1 2 1 ε = 0.001 56q 3 30.33 13q 7 19.52 12q 7 37 3.5 Superconvergence In this section, we study the superconvergence property of the DG solution at Radau points. 3.5.1 Superconvergence at Radau points In this subsection, we study the superconvergence property of DG solutions at Radau points. The local basis functions of P k (Ij ) are chosen to be the Lagrangian polyno- mials, based on the following k + 1 points ζk,l xj+rl = xj + ∆x, l = 0, · · · , k 2 where {ζk,l } are the roots of the Radau polynomial Pk+1 (x) − Pk (x) and Pk (x) is the Legendre polynomial of degree k. Figure 3.5 shows the local basis functions of P 1 based on Radau points. φj− 61 φj+ 12 @ @ @ @ c c c @c c c xj− 1 @ 6 xj− 23 xj− 21 xj+ 21 xj+ 23 Figure 3.5: Local basis function of P 1 based on Radau points. Table 3.2 displays {ζk,l} for k = 1, 2, 3. Note that xj+rk = xj+ 1 is the downwind 2 point of the element Ij . Remark 3.5.1. The basis functions in this section are different from those discussed 38 Table 3.2: {ζk,l }kl=0 : the roots of Pk+1 (x) − Pk (x) ζk,l k = 1 ζ1,0 = − 31 , √ ζ1,1 = 1 √ k = 2 ζ2,0 = − 1+5 6 , ζ2,1 = − 1−5 6 , ζ2,2 = 1 ζ3,0 = −0.82282408097459210520890771246109 ζ = −0.18106627111853057827014749586234 k = 3 3,1 ζ3,2 = 0.57531892352169411205048377975200 ζ3,3 = 1 in Section 3.4 and the amplification matrices are different, but they have the same eigenvalues. These eigenvalues are also the same as the ones derived in [55] using an orthonormal basis of Legendre polynomials. This can be proved as follows: Assume  k+1  k+1 the amplification matrices with the basis functions φlj l=1 and ϕlj l=1 are G1 and G2 , respectively. Then G1 = T −1 G2 T  k+1  k+1 where T is the change-of-basis matrix from φlj l=1 to ϕlj l=1 and obviously T is invertible. G1 and G2 are similar matrices and thus they have the same characteristic polynomials. As in Section 3.3, the DG method with the above basis can be written in a matrix form. The standard Fourier analysis is used to find the explicit solution of the DG scheme as discussed in Section 3.4. We do not repeat the details of this analysis procedure. Denote kerl k∞ := max U(xj+rl , t) − u(xj+rl , t) , 1≤j≤N ke− 1 k∞ := max U(x− 1 , t) − u(x− 1 , t) . + 2 1≤j≤N j+ 2 j+ 2 The followings are the error estimate results. 39 • k=1 r 1 16 2 ker0 k∞ = + t ξ 3 + O(ξ 4), 72 81 r 1 16 2 ke− k +1 ∞ = + t ξ 3 + O(ξ 4); 2 72 9 • k=2 t 5 ker0 k∞ = 0.00112979589711327123927891362988ξ 4 − ξ + O(ξ 6), 7200 t 5 ker1 k∞ = 0.00082979589711327123927891362988ξ 4 + ξ + O(ξ 6), q 7200 2 t + 144/25 ke− k + 21 ∞ = ξ 5 + O(ξ 6); 7200 • k=3 ker0 k∞ = 4.58526548614820367915101051257673 × 10−5ξ 5 + 1.48442413127661313280202433713330 × 10−6 ξ 7 −2.95523546192529651029864723686194 × 10−7 tξ 8 + O(ξ 9), ker1 k∞ =4.81347804829324175268316712031863 × 10−5 ξ 5 − 1.74662334589086702133200039227532 × 10−6ξ 7 + 1.12983256993051227300311527010632 × 10−8 tξ 8 + O(ξ 9), ker2 k∞ =2.60817329279298238245105944584508 × 10−5 ξ 5 + 1.86838652276109018545310450926279 × 10−6 ξ 7 + 3.23693699988152665996155459585616 × 10−7 tξ 8 + O(ξ 9), q 2 t + 576/49 7 ke− k +1 ∞ = ξ 2 1412200 3.20791891379522690550503036041469 × 10−7 t 8 + p 2 ξ + O(ξ 9). 576 + 49t 40 3.5.2 Superconvergence at the Downwind Point In this subsection, we take a particular look at the error at the downwind point of each cell. Let Mkf be the number of points per wavelength needed to guarantee the error ke− k ≤ ε when using P k polynomials. According to the results derived in the +1 ∞ 2 previous subsection, even when t ≫ 1, we have • k=1 r  3 1 16 2 t πq 4π ke− k + 12 ∞ = + t ξ 3 + O(ξ 4) ≃ ξ = ≤ ε. 72 9 72 36 M1f • k=2 q 2  5 t + 144/25 t 5 πq 6π ke− k + 21 ∞ = 5 6 ξ + O(ξ ) ≃ ξ = ≤ ε. 7200 7200 3600 M2f • k=3 q 2  7 t + 576/49 t πq 8π ke− k + 12 ∞ = 7 8 ξ + O(ξ ) ≃ ξ7 = ≤ ε. 1412200 1411200 705600 M3f Therefore,  π  13  q  31 M1f ≥ 4π , 36 ε  π  51  q  15 M2f ≥ 6π , 3600 ε  π  17  q  71 M3f ≥ 8π . 705600 ε 41 Table 3.3 shows the lower bound of Mkf as a function of q for three representative values of ε. Table 3.3: Lower bound of Mpf as a function of q M1f M2f M3f 1 1 1 ε = 0.1 12q 3 7q 5 6q 7 1 1 1 ε = 0.01 25q 3 11q 5 8q 7 1 1 1 ε = 0.001 55q 3 18q 5 11q 7 3.6 Fully Discrete Schemes The previous sections deal with the semi-discrete schemes. In practice, the time discretization also plays an important role. An analysis of fully discretized approxi- mation will be presented in this section with the basis functions discussed in Section 3.4. 3.6.1 The Case of P 1 In this subsection, we present the details of the fully discretized error analysis for the piecewise linear case. Let u˜ be the approximate solution of the differential-difference equations (3.11) obtained by using p-stage explicit Runge-Kutta methods of order p (RK(p, p)) [32]. Then     iωxj− 1 u˜j− 1 e 4 t  4  n   =R  , n= (3.45) iωxj+ 1 ∆t u˜j+ 1 e 4 4 42 where ∆t2 2 ∆tp p R = 1 + ∆tG + G +···+ G (3.46) 2! p! with the amplification matrix G defined in (3.20). Let Q = [V1 , V2 ] be the matrix with G’s eigenvectors as columns. Clearly, we have    λ1 0  Q−1 GQ = Λ =  . (3.47) 0 λ2 e = Q−1 RQ. From (3.47) and (3.46), we have Denote R p e = Q−1 RQ = 1 + ∆tΛ + · · · + ∆t Λp R p!   ∆tp p  1 + ∆tλ 1 + · · · λ 0   p! 1  = .  ∆tp p  0 1 + ∆tλ2 + · · · λ p! 2 Thus, e −1 )n = QR Rn = (QRQ en Q−1  n ∆tp p  1 + ∆tλ 1 + · · · λ 0  (3.48)  p! 1  −1 = Q  Q .  ∆tp p  0 1 + ∆tλ2 + · · · λ p! 2 Substituting (3.48) into (3.45), we can determine the explicit solution of the RKDG scheme for solving (3.8). Comparing this with the exact solution (3.9) will give us ∆t the quantitative error estimates, based on the assumption cfl = a ∆x . By a simple Taylor expansion, the error between the DG solution with an RK(2,2) 43 time discretization and the exact solution is ξ2  ξ2  at − 3i(−1+6 cfl)t u˜j− 1 − U(xj− 1 , t) = 1 − 4icfl2 t e−it − 1 − 6 cfl + 18 cfl2 cfl ∆x e 1−6 cfl+18 cfl2 4 4 24 24 ξ3  + −i − 8t − 24cfl2 + 72cfl3 t e−it 576 ξ3  at − 3i(−1+6 cfl)t + i 1 − 6 cfl + 18 cfl2 cfl ∆x e 1−6 cfl+18 cfl2 576 ξ3  at − 3i(−1+6 cfl)t + t 1 − 6 cfl + 18 cfl2 cfl ∆x e 1−6 cfl+18 cfl2 + · · · 24 where ξ = ω∆x and t = ωat.  cflat∆x In order to let the terms with coefficients 1 − 6 cfl + 18 cfl2 go to 0 as ∆x → 0, we assume 1 − 6 cfl + 18 cfl2 < 1, i.e. cfl < 31 . Then   −it ξ2 iξ 2 2 iξ 3 ξ3 ξ3 2 ξ3 3 4 u˜j− 1 − U(xj− 1 , t) = e − cfl t − − t− cfl t + cfl t + O(ξ ) . 4 4 24 6 576 72 24 8 Thus ke− 1 k∞ = max u˜j− 1 − U(xj− 1 , t) 4 1≤j≤N 4 4 q  3 (3.49) 1 2t + 7cfl2 t − 18cfl3 t ξ 2 = 1 + 16cfl4 t ξ 2 − p + O(ξ 4). 24 144 1 + 16cfl t 4 2 We keep the same notations here as in the semi-discrete case. Similarly, q  1 2t + 7cfl2 t − 18cfl3 t ξ 3 4 2 2 ke+ 1 k∞ = 1 + 16cfl t ξ + p + O(ξ 4) (3.50) 4 24 144 1 + 16cfl t 4 2 44 and e1 = max{ke− 1 k∞ , ke+ 1 k∞ } 4 4 q  (3.51) 1 4 2 22t + 7cfl2 t − 18cfl3 t ξ 3 = 1 + 16cfl t ξ + p + O(ξ 4). 24 144 1 + 16cfl t 4 2 Since it takes a minimum of two points per wavelength to uniquely specify a wave, the largest wave number that can be represented on the N cells (reminder: 2N two points are used for each cell) is ω = 2 = N. Therefore, we consider simple waves with 0 < ω ≤ N. If ω is much less than N in magnitude, then there are many grid points per wavelength, and the solution is well resolved. In this case, ξ is a small number, and we can neglect O(ξ 3) and obtain q 1 2 e1 ≃ 1 + 16cfl4 t ξ 2 . (3.52) 24 If ξ is large, corresponding to fewer grid points per wavelength in the solution, then the P 1 solution is not a good approximation of the exact solution and approximation polynomials of degree > 1 are required. Similarly, when using an RK(3,3) time discretization, the error estimates under the assumption1 cfl < 0.418, 2 ξ ξ 3  ke− 1 k∞ = − t 1 + 3cfl + O(ξ 4), 3 (3.53) 4 24 72 ξ2 ξ3  ke+ 1 k∞ = + t 1 + 3cfl3 + O(ξ 4). (3.54) 4 24 72 −2+ cflat∆x 1 in order to let the terms with coefficients 1 − 6 cfl + 18 cfl2 − 36cfl3 → 0 as ∆x → 0. 45   At t = O 1ξ , ξ2 ξ3  e1 ≃ + t 1 + 3cfl3 . (3.55) 24 72 The necessary number of points per wavelength required to ensure a specific error ε can be obtained by rewriting (3.52) and(3.55) in terms of M1 and q and setting e1 ≤ ε. 3.6.2 The Case of P 2 For k = 2 and k = 3, the analysis is similar. In this subsection, we present the fully discretized error estimates when piecewise P 2 polynomials are used. When using an RK(3,3) time discretization, under the assumption cfl < 0.327, we have r 1 7222932 6 2 3 ke− 1 k∞ = 1+ cfl t ξ + O(ξ 4). (3.56) 3 1296 3125 q 1 2 ke0 k∞ = 1 + 100cfl6 t ξ 3 + O(ξ 4), (3.57) 240 r 23 72900 6 2 3 ke+ 1 k∞ = 1+ cfl t ξ + O(ξ 4). (3.58) 3 6480 529 Thus, q 1 2 e2 = max{ke− 1 k∞ , ke0 k∞ , ke+ 1 k∞ } ≃ 1 + 100cfl6 t ξ 3 . (3.59) 3 3 240 46 When using an RK(4,4) time discretization, under the assumption cfl < 0.351, we have 1 3 t ke− 1 k∞ = ξ + cfl4 ξ 4 + O(ξ 5), (3.60) 3 1296 120 1 3 t ke0 k∞ = ξ + cfl4 ξ 4 + O(ξ 5), (3.61) 240 120 23 3 t ke+ 1 k∞ = ξ − cfl ξ + O(ξ 5). 4 4 (3.62) 3 6480 120   At t = O 1ξ , 1 3 t e2 ≃ ξ + cfl4 ξ 4 . (3.63) 240 120   1 While at t = O ξ2 , we need more terms of the Taylor series to approximate e2 . In this case, e2 ≃ 0.0041666666666667ξ 3 + 0.0083333333333333cfl4 ξ 4 t − 0.0003201517489712ξ 5 + ξ 6 −6.0626102292768959 × 10−6 t − 0.0000680298353909cfl4 t  +0.0008873456790123cfl5 t − 0.0029761904761905cfl6 t  (3.64) 2 + ξ 7 −0.00013365952729388 + 2.3148148148148148 × 10−6 t 2 2 −0.0000354938271605cfl4 t + 0.0002314814814815cfl5 t 2 2 +0.0001360596707819cfl8 t − 0.0017746913580247cfl9 t  2 +0.005787037037037cfl10 t . The necessary number of points per wavelength required to ensure a specific error ε can be obtained by rewriting (3.59), (3.63) and (3.64) in terms of M2 and q and setting e2 ≤ ε. 47 3.6.3 The Case of P 3 For the cases k = 1 and k = 2, the cfl condition is consistent with the CFLL2 table in [26]. However, for k = 3, we cannot obtain the cfl condition using the same procedure by Taylor expansion and Mathematica due to limited memory. In this subsection, we provide the error estimates based on a fixed cfl = 0.1. When using an RK(4,4) time discretization, we have q 2 ke− 3 k∞ = ξ 3.2541560538021135 × 10−8 + 6.944444444444444 × 10−13 t + O(ξ 5), 4 8 q 2 ke− 1 k∞ = ξ 1.3128715187393769 × 10−9 + 6.944444444444444 × 10−13 t + O(ξ 5), 4 8 q 2 ke+ 1 k∞ = ξ 1.3227188998669514 × 10−7 + 6.944444444444444 × 10−13 t + O(ξ 5), 4 8 q 2 ke+ 3 k∞ = ξ 9.3484606662364834 × 10−9 + 6.944444444444444 × 10−13 t + O(ξ 5). 4 8 And q 4 2 e3 ≃ ξ 1.3227188998669514 × 10−7 + 6.9444444444444444 × 10−13 t . (3.65) When using an RK(5,5) time discretization, we have ke− 3 k∞ = 0.0001803927951389ξ 4 + 1.3888888888888889 × 10−8 tξ 5 + O(ξ 6), 8 ke− 1 k∞ = 0.0000362335689484ξ 4 + 1.3888888888888889 × 10−8 tξ 5 + O(ξ 6), 8 ke+ 1 k∞ = 0.0003636920262897ξ 4 + 1.3888888888888889 × 10−8 tξ 5 + O(ξ 6), 8 ke+ 3 k∞ = 0.0000966874379960ξ 4 + 1.3888888888888889 × 10−8 tξ 5 + O(ξ 6). 8 48   At t = O ξ12 , e3 ≃ 0.0003636920262897ξ 4 + 1.3888888888888889 × 10−8 tξ 5 . (3.66) The necessary number of points per wavelength required to ensure a specific error ε can be obtained by rewriting (3.65) and (3.66) in terms of M3 and q and setting e3 ≤ ε. Table 3.4 shows the Mk as a function of q (we only present the leading term here), derived from the error estimates in each subsection. These estimates is computed at   t = O 1ξ . Here we fix cfl = 0.3 for k = 1, cfl = 0.2 for k = 2 and cfl = 0.1 for k = 3. Table 3.4: Leading term of the lower bound of Mk as a function of q M1 M2 M3 RK(p, p) RK(2,2) RK(3,3) RK(3,3) RK(4,4) RK(4,4) RK(5,5) 1 1 1 1 1 1 ε = 0.1 12q 2 12q 3 5q 3 3q 4 2q 4 2q 5 1 1 1 1 1 1 ε = 0.01 39q 2 26q 3 11q 3 8q 4 4q 4 2q 5 1 1 1 1 1 1 ε = 0.001 122q 2 57q 3 24q 3 10q 4 8q 4 4q 5 3.6.4 A Measure of Work per Wavelength So far, we have expressed ek as a function of Mk and q. Now the work per wavelength to obtain the error ek when integrating to time t using RK(p, p) is given by t Wk = 2(k + 1) × Mk × × p, ∆t which is a product of the number of operations per mesh point per time step 2(k +1), 2π(k+1) the number of points per wavelength Mk = ω∆x , the total number of time steps 49 t ∆t ∆t and the number of stages per time step p. Using ωat = 2πq and cfl = a ∆x , we obtain t Wk = 2(k + 1)Mk p ∆t 2πq = 2(k + 1)Mk p ωa∆t 2πq = 2(k + 1)Mk p ω∆x cfl q 2π = 2(k + 1)Mk p cfl ω∆x q Mk = 2(k + 1)Mk p cfl k + 1 2 = p qMk2 . cfl 3.7 Numerical Results In this section, we provide numerical experiments to demonstrate the predicted re- sults presented in the previous sections. To verify the predicted results for the semi-discrete scheme, we adopt SSPRK(9,9) [32] to make the temporal error negligible compared to the spatial error. This time discretization method for solving du/dt = Lu, where L is a spatial discretization operator (for this SSP method to be ninth order, L needs to be linear) is defined as 50 follows: u(i) = u(i−1) + ∆t Lu(i−1) , i = 1, · · · , 8 7 X  u(9) = α9,k u(k) + α9,8 u(8) + ∆t L u(8) , k=0 where 16687 2119 103 α9,0 = , α9,0 = , α9,2 = , 45360 5760 560 53 11 1 α9,3 = , α9,4 = , α9,5 = , 864 720 320 1 1 1 α9,6 = , α9,7 = , α9,8 = . 2160 10080 362880 ∆t Uniform meshes are used in the calculation and the CFL condition cfl = a ∆x is chosen from the table in [26]. Tables 3.5-3.7 list the error between DG solution and exact solution for fixed N t with different cell size N, where N = ω . The numerical results are computed using P k polynomials and RK(9,9) time discretization. The predicted results are computed using the formulas derived in Section 3.4. The predicted errors agree with the computed results very well. These tables also show that the orders of the errors are different at different final time t. This is because the higher order terms of the error dominate for large t. Figures 3.6-3.8 show the time evolution of the L∞ error and ek error when using P k polynomials with three different time discretization methods. One is the same order as the spatial error, another one is one order higher than the spatial error and the last one is the RK(9,9) which makes the temporal error negligible comparing with 51 Table 3.5: When using P k polynomials on a uniform mesh of N cells, t = 1 N Numerical results Predicted by analysis k=1 e1 order e1 order 20 4.4639E-03 4.5430E-03 40 1.0763E-03 2.05 1.0819E-03 2.07 80 2.6333E-04 2.03 2.6375E-04 2.04 160 6.5072E-05 2.02 6.5096E-05 2.02 k=2 e2 order e2 order 20 1.2796E-04 1.2825E-04 40 1.6139E-05 2.99 1.6119E-05 2.99 80 2.0179E-06 3.00 2.0177E-06 3.00 160 2.5231E-07 3.00 2.5230E-07 3.00 k=3 e3 order e3 order 20 3.6972E-06 3.5497E-06 40 2.2928E-07 4.01 2.2153E-07 4.00 80 1.3858E-08 4.05 1.3840E-08 4.00 160 8.6484E-10 4.00 8.6494E-10 4.00 Table 3.6: When using P k polynomials on a uniform mesh of N cells, t = 100 N Numerical results Predicted by analysis k=1 e1 order e1 order 20 4.5526E-02 4.7177E-02 40 6.3632E-03 2.84 6.4111E-03 2.88 80 9.2878E-04 2.78 9.2990E-04 2.79 160 1.4831E-04 2.65 1.4837E-04 2.65 k=2 e2 order e2 order 20 1.3370E-04 1.3463E-04 40 1.6120E-05 3.05 1.6164E-05 3.06 80 2.0172E-06 3.00 2.0180E-06 3.00 160 2.5226E-07 3.00 2.5230E-07 3.00 k=3 e3 order e3 order 20 3.5443E-06 3.5570E-06 40 2.2147E-07 4.00 2.2148E-07 4.01 80 1.3831E-08 4.00 1.3838E-08 4.00 160 8.6229E-10 4.00 8.6490E-10 4.00 52 Table 3.7: When using P k polynomials on a uniform mesh of N cells, t = 1000 N Numerical results Predicted by analysis k=1 e1 order e1 order 20 3.5039E-01 4.3476E-01 40 5.3261E-02 2.72 5.4858E-02 2.99 80 6.9551E-03 2.94 6.9858E-03 2.97 160 9.0479E-04 2.94 9.0535E-04 2.95 k=2 e2 order e2 order 20 4.4844E-04 8.2153E-04 40 2.0799E-05 4.43 2.1489E-05 5.26 80 2.0575E-06 3.34 2.0589E-06 3.38 160 2.5111E-07 3.03 2.5261E-07 3.03 k=3 e3 order e3 order 20 3.7047E-06 3.7487E-06 40 2.2252E-07 4.06 2.2298E-07 4.07 80 1.3902E-08 4.00 1.3850E-08 4.01 160 8.6496E-10 4.01 8.6499E-10 4.00 the spatial error. The numerical L∞ error is computed by taking a uniform partition of each element with 20 points. These figures show that the time discretization method needs to be at least one order higher than the DG method, in order to get the same result as the semi-discrete case. The order of the RK method higher than k + 1 makes little difference when the DG method is using P k polynomials. Also, the figures demonstrate that the predicted errors in general agree with the computed results very well. Tables 3.8-3.10 list the numerical errors and their orders at Radau points for different final time t. These tables verify the results derived in Section 3.5: the errors at the downwind point of each element are superconvergent of order 2k + 1 and at other Radau points are superconvergent of order k + 2 when using piecewise P k polynomials. These tables also show the relation between the error and t: for k = 1, errors are proportional to t while for k = 2, 3, the leading terms of the errors   do not depend on t, thus the errors are almost the same for t under O 1ξ . We have 53 The error e1 ∞ Numerical L error 0.16 0.16 RK(2,2) Numerical RK(2,2) Predicted 0.14 0.14 RK(3,3) Numerical RK(2,2) RK(3,3) Predicted 0.12 RK(3,3) 0.12 RK(9,9) Numerical RK(9,9) Semi-dicrete error 0.10 0.10 0.08 0.08 0.06 0.06 0.04 0.04 0.02 0.02 0.00 0.00 0 50 100 150 200 250 300 0 100 200 300 400 time time Figure 3.6: Time evolution of L∞ error (left) and e1 error (right) for P 1 DG solution with N = 40, cfl = 0.3. 7E-05 The error e2 7E-05 ∞ Numerical L error RK3,3) Numerical 6E-05 RK(3,3) Predicted 6E-05 RK(4,4) Numerical RK(4,4) Predicted 5E-05 RK(3,3) 5E-05 RK(9,9) Numerical RK(4,4) Semi-dicrete error RK(9,9) 4E-05 4E-05 3E-05 3E-05 2E-05 2E-05 1E-05 1E-05 0 0 0 100 200 300 400 0 100 200 300 400 time time Figure 3.7: Time evolution of L∞ error (left) and e2 error (right) for P 2 DG solution with N = 80, cfl = 0.2. 3.5E-09 Numerical L∞ error 3.5E-09 The error e3 3.0E-09 3.0E-09 RK(4,4) Numerical RK(4,4) Predicted RK(5,5) Numerical 2.5E-09 RK(4,4) 2.5E-09 RK(5,5) Predicted RK(5,5) RK(9,9) Numerical RK(9,9) Semi-dicrete error 2.0E-09 2.0E-09 1.5E-09 1.5E-09 1.0E-09 1.0E-09 0 200 400 600 800 1000 200 400 600 800 1 time time Figure 3.8: Time evolution of L∞ error (left) and e3 error (right) for P 3 DG solution with N = 160, cfl = 0.1. 54 also tested the P 2 case at final time t = 500, the error max {kerl k∞ } achieves fifth 0≤l≤1 order. This is because as t increases to 500, the fifth order term dominates. This verifies the fact that during different time intervals, the dominant terms are different. We have also used a non-uniform mesh which is 60% random perturbation of the uniform mesh. For example, the right end of the cell Ij is now xj+ 1 + 60%(rj+ 1 − 2 2 0.5)∆x, where xj+ 1 , ∆x are taking values under the uniform mesh and rj+ 1 is a 2 2 random number from the uniform distribution over the range (0, 1). Tables 3.11- 3.13 list the errors and their orders at Radau points in this case for different final time t, when using piecewise P k polynomials. We can see that the error at Radau points is of order around k + 2, including the downwind point. We have also tested the cases of P 2 and P 3 using 1% random perturbation of the uniform perturbation. The numerical results still do not show the strong superconvergence of order 2k + 1. Table 3.8: Error at Radau points using P 1 DG and RK(9,9) on a uniform mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 4.75E-04 4.25E-03 4.17E-02 40 5.91E-05 3.01 5.37E-04 2.98 5.36E-03 2.96 ker0 k∞ 80 7.38E-06 3.00 6.73E-05 3.00 6.72E-04 3.00 160 9.22E-07 3.00 8.42E-06 3.00 8.41E-05 3.00 20 6.97E-04 4.28E-03 4.18E-02 ke− k 40 8.87E-05 2.97 5.39E-04 2.99 5.35E-03 2.97 + 21 ∞ 80 1.12E-05 2.99 6.77E-05 2.99 6.72E-04 2.99 160 1.40E-06 3.00 8.48E-06 3.00 8.41E-05 3.00 Tables 3.14 and 3.15 list the necessary number of points per wavelength required to guarantee ek ≤ ε for various RKDG schemes and values of ε. Denote RKDG(p, k+ 1) to be the scheme using the RK(p, p) method for time discretization and the DG method with piecewise P k polynomials for space discretization. For these schemes with RK(9,9) method, the results of semi-discrete case are used as the predicted results. These two tables also verify our conclusions about the error estimates with 55 Table 3.9: Error at Radau points using P 2 DG and RK(9,9) on a uniform mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 1.05E-05 1.23E-05 4.97E-05 max {kerl k∞ } 40 6.73E-07 3.96 6.37E-07 4.27 1.83E-06 4.76 0≤l≤1 80 4.25E-08 3.98 3.88E-08 4.07 7.30E-08 4.65 160 2.67E-09 3.99 2.56E-09 3.92 3.27E-09 4.48 20 1.09E-06 4.34E-06 4.23E-05 ke− k 40 3.44E-08 4.98 1.36E-07 4.99 1.33E-06 5.00 + 21 ∞ 80 1.08E-09 5.00 4.26E-09 5.00 4.15E-08 5.00 160 3.37E-11 5.00 1.33E-10 5.00 1.30E-09 5.00 Table 3.10: Error at Radau points using P 3 DG and RK(9,9) on a uniform mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 1.64E-07 1.46E-07 1.48E-07 max kerl k∞ 40 4.70E-09 5.12 4.60E-09 4.99 4.60E-09 5.01 0≤l≤2 80 1.44E-10 5.02 1.44E-10 5.00 1.44E-10 5.00 160 4.50E-12 5.00 4.50E-12 5.00 4.50E-12 5.00 20 1.86E-08 2.24E-09 2.13E-08 ke− k 40 2.12E-10 6.45 1.76E-11 6.99 1.67E-10 7.00 + 21 ∞ 80 5.09E-13 8.70 1.38E-13 7.00 1.31E-12 7.00 160 3.75E-16 10.41 1.08E-15 7.00 1.02E-14 7.00 320 2.85E-18 7.04 8.43E-18 7.00 7.98E-17 7.00 640 2.23E-20 7.00 6.61E-20 6.99 6.23E-19 7.00 Table 3.11: Error at Radau points using P 1 DG and RK(9,9) on a random mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 1.82E-03 6.73E-03 6.33E-02 40 2.39E-04 2.93 7.68E-04 3.13 7.50E-03 3.08 ker0 k∞ 80 3.23E-05 2.89 9.05E-05 3.08 8.87E-04 3.08 160 3.11E-06 3.37 1.22E-05 2.89 1.20E-04 2.89 320 4.98E-07 2.64 1.49E-06 3.03 1.42E-05 3.08 20 1.19E-03 6.81E-03 6.34E-02 40 1.50E-04 2.98 7.80E-04 3.13 7.51E-03 3.08 − k ke+ 1 ∞ 80 2.27E-05 2.72 8.96E-05 3.12 8.87E-04 3.08 2 160 2.36E-06 3.27 1.22E-05 2.88 1.20E-04 2.89 320 3.51E-07 2.75 1.43E-06 3.09 1.42E-05 3.08 56 Table 3.12: Error at Radau points using P 2 DG and RK(9,9) on a random mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 3.20E-05 3.14E-05 1.31E-04 40 2.63E-06 3.61 2.58E-06 3.61 4.63E-06 4.82 ker0 k∞ 80 1.55E-07 4.08 1.99E-07 3.70 2.11E-07 4.45 160 1.32E-08 3.56 1.37E-08 3.86 1.11E-08 4.25 320 7.70E-10 4.10 7.85E-10 4.12 7.25E-10 3.94 20 1.81E-05 1.39E-05 9.31E-05 40 4.87E-07 5.21 5.20E-07 4.74 2.81E-06 5.05 ke− 1 k∞ +2 80 6.06E-08 3.01 3.66E-08 3.83 9.28E-08 4.92 160 2.34E-09 4.69 1.25E-09 4.87 2.74E-09 5.08 320 1.33E-10 4.14 6.55E-11 4.26 1.01E-10 4.76 Table 3.13: Error at Radau points using P 3 DG and RK(9,9) on a random mesh of N cells t=1 t = 10 t = 100 N error order error order error order 20 3.90E-07 4.16E-07 4.71E-07 40 4.73E-08 3.04 3.55E-08 3.55 3.74E-08 3.65 ker0 k∞ 80 1.10E-09 5.43 8.65E-10 5.36 1.04E-09 5.17 160 2.71E-11 5.34 3.18E-11 4.76 2.94E-11 5.14 320 1.25E-12 4.44 1.43E-12 4.47 1.40E-12 4.39 20 2.29E-07 2.71E-08 4.51E-08 40 1.80E-08 3.67 8.42E-09 1.69 3.91E-09 3.53 ke− 1 k∞ +2 80 2.97E-10 5.92 1.16E-10 6.19 4.84E-11 6.34 160 8.42E-12 5.14 3.56E-12 5.02 1.25E-12 5.27 320 2.49E-13 5.08 1.05E-13 5.08 4.71E-14 4.73 57 different time discretization methods, which states that the RK method needs to be at least one order higher than the DG method, in order to get the same result as the semi-discrete case; and there is little difference among the RK methods of order higher than k + 1 when the DG method is using piecewise P k polynomials. Table 3.14: Necessary number of points per wavelength for ek ≤ ε when time period q=1 ε = 0.1 ε = 0.01 ε = 0.001 RKDG(p, k + 1) Numerical Predicted Numerical Predicted Numerical Predicted RKDG(2,2) 15.02 14.63 42.87 42.47 130.06 129.81 RKDG(3,2) 13.43 14.09 34.32 34.62 92.60 92.73 RKDG(9,2) 13.14 13.82 33.83 34.14 91.86 91.99 RKDG(3,3) 8.00 6.78 15.74 14.62 31.69 31.49 RKDG(4,3) 7.61 6.66 14.37 14.20 30.17 30.46 RKDG(9,3) 7.63 6.53 14.46 14.08 30.00 30.33 RKDG(4,4) 6.16 6.17 11.10 10.98 19.25 19.52 RKDG(5,4) 6.16 6.17 11.10 10.98 19.24 19.52 RKDG(9,4) 6.16 6.17 11.10 10.98 19.24 19.52 Table 3.15: Necessary number of points per wavelength for ek ≤ ε when time period q = 10 ε = 0.1 ε = 0.01 ε = 0.001 RKDG(p, k + 1) Numerical Predicted Numerical Predicted Numerical Predicted RKDG(2,2) 39.46 38.60 122.60 122.06 386.36 385.97 RKDG(3,2) 26.74 27.38 60.85 61.03 140.86 140.93 RKDG(9,2) 26.08 26.72 59.48 59.67 138.16 138.23 RKDG(3,3) 13.47 11.27 26.07 24.27 53.43 52.30 RKDG(4,3) 11.38 13.58 18.83 19.86 32.85 33.11 RKDG(9,3) 11.51 13.53 18.94 19.64 32.47 32.32 RKDG(4,4) 8.00 6.19 12.65 11.00 20.45 19.57 RKDG(5,4) 8.00 6.19 12.64 10.99 20.35 19.53 RKDG(9,4) 8.00 8.70 12.64 12.65 20.36 20.35 Tables 3.16-3.17 list the work per wavelength of the scheme RKDG(p, k + 1), wk , required to obtain the error ε when time period q = 1 and q = 10. We should point out that we could have increased the order of the RK method used in the scheme RKDG(k + 2, k + 1) without significantly affecting the number of points per 58 wavelength, but the numbers of stages per timestep would increase considerably and so does the work per wavelength. It is clear that even for short time period, high order methods are the most appropriate choice when accuracy is the primary consideration. When using P k polynomials for the DG method, RK(k + 1, k + 1) needs less work than the RK(k + 2, k + 2) for the short time period. While for long time periods, the better choice is RK(k + 2, k + 2). Table 3.16: When time period q = 1 Type of Multiplies Stages Wk Order Scheme per Meshpoint per Timestep ε = 0.1 ε = 0.01 ε = 0.001 RKDG(2,2) 2 4 2 2,900 24,000 224,700 RKDG(3,2) 2 4 3 4,000 24,000 172,000 RKDG(3,3) 3 6 3 1,400 6,400 29,800 RKDG(4,3) 3 6 4 1,800 8,100 37,100 RKDG(4,4) 4 8 4 3,000 9,600 30,500 RKDG(5,4) 4 8 5 3,800 12,000 38,100 Table 3.17: When time period q = 10 Type of Multiplies Stages Wk Order Scheme per Meshpoint per Timestep ε = 0.1 ε = 0.01 ε = 0.001 RKDG(2,2) 2 4 2 201,400 1,995,300 19,891,700 RKDG(3,2) 2 4 3 149,900 745,000 3,972,500 RKDG(3,3) 3 6 3 38,100 176,800 820,400 RKDG(4,3) 3 6 4 51,800 157,800 438,600 RKDG(4,4) 4 8 4 30,600 96,900 306,300 RKDG(5,4) 4 8 5 38,300 120,800 381,500 In Figure 3.9, we compare the necessary number of points per wavelength required to obtain a fixed error 0.01 between two second order accurate methods. One is the DG method with P 1 polynomials and the other is the second order finite difference method. It can be seen clearly that, even for a short time period, the DG method needs much fewer points per wavelength. Also, the necessary number of points per 59 wavelength of the finite difference scheme grows faster with time than the RKDG method. This is because the leading term of the error for the RKDG method does not depend on t, while the finite difference method does. This conclusion also holds for higher order schemes between the DG method and the finite difference method with the same order. In [65], Swartz and Wendroff computed the necessary number of intervals per wavelength required to obtain a fixed error for the finite element method using smooth splines as basis functions and various time discretization methods, such as the trapezoidal method and the leap-frog method. With the same order of time discretization method and space discretization method, the RKDG method requires fewer intervals per wavelength than the finite element method with smooth splines as basis functions to obtain a specified error. Moreover, the leading term of the error of this finite element method using smooth splines depends on t. This means the necessary number of intervals per wavelength for the finite element method using smooth splines also grows faster with time than the RKDG method. ε=0.01 nd 400 2 FD 1 P DG 300 200 100 0 0 10 20 30 40 q Figure 3.9: Necessary points per wave for kek∞ ≤ ε Remark 3.7.1. All the numerical results and the predicted results take point value collocation as initial conditions in this chapter. The usual way of taking initial conditions in a finite element method is via an L2 projection. However this does not affect the results in our chapter. See Figure 3.10. The following are quatitative eror 60 for k = 1: 1 2 t 3 ke− 1 k∞ = ξ − ξ + O(ξ 4) 4 32 72 5 t ke+ 1 k∞ = ξ 2 − ξ 3 + O(ξ 4) 4 96 72 0.006 0.03 point value point value 2 2  L projection  L projection  0.004  0.02      ||e-1/4||∞          0.002  ||e-1/4||∞ 0.01                        0.00     0 0 100 200 300 0 20 40 60 80 100 time time 0.04 0.010 point value point value 2 2 0.03  L projection 0.008  L projection    0.006  0.02   ||ej+1/4||∞                 0.004  ||ej+1/4||∞ 0.01                     0.002    0.00 0.000 0 100 200 300 0 50 100 time time Figure 3.10: Comparison of time evolution of the errors with two different initial interpolations. Remark 3.7.2. We also did this Fourier analysis procedure using central flux. Both our theoritical results and numerical results is only k-th order accurate when using DG methods with piecewise P k polynomials. 61 3.8 Concluding Remarks In this chapter, by Fourier analysis, we have derived the quantitative error estimates of the DG methods using piecewise P k polynomials with 1 ≤ k ≤ 3, for solving time dependent linear wave problems. We have proved the superconvergence property of the DG solution at Radau points. The error of the DG solution shows order 2k + 1 superconvergence at the downwind point of each element and order k + 2 supercon- vergence at other Radau points. We also provide a fully discretized error analysis with various Runge-Kutta methods. We have computed the necessary number of points per wavelength required to obtain a fixed error for several RKDG schemes. A very important implication of the discussion in this chapter is that the dominant terms of the error for the RKDG schemes are different during different time intervals. This further justifies the advantage of choosing DG methods for long time simulation of linear wave equations. The technique of Fourier analysis discussed in this chapter can be extended to the hyperbolic systems with constant coefficients. Theoretically, it can also be extended to the multidimensional case, however the algebraic manipulations may become pro- hibitively complicated. Chapter 4 A Simple Weighted Essentially Non-Oscillatory Limiter for Runge-Kutta Discontinuous Galerkin Methods 62 63 4.1 Introduction In this chapter, we consider the following hyperbolic conservation law ut + f (u)x = 0, (4.1) u(x, 0) = u0 (x), and its two-dimensional version, where u and f (u) can be either scalars or vectors. We investigate a simple limiter using weighted essentially non-oscillatory (WENO) methodology for the Runge-Kutta discontinuous Galerkin (RKDG) methods, with the goal of obtaining a robust and high order limiting procedure to simultaneously maintain uniform high order accuracy in smooth regions and control spurious nu- merical oscillations near discontinuities. The idea of this limiter is to reconstruct the entire polynomial based on the polynomials of the DG solution in the target and neighboring cells, instead of reconstructing point values or moments based on cell averages or lower order moments. The main difficulty in solving (4.1) is that solutions may contain discontinuities even if the initial conditions are smooth. DG methods can compute solutions to (4.1), which are either smooth or have weak shocks and other discontinuities, with- out further modification. However, for problems containing strong discontinuities, the scheme will generate significant oscillations and even nonlinear instabilities. We often need to apply nonlinear limiters to control these oscillations. Many such lim- iters exist in the literature, such as the minmod type total variation bounded (TVB) limiter [24, 23, 21, 20, 25], the moment-based limiter [13] and the more recent im- proved moment limiter [14]. Although these limiters can control spurious numerical oscillations near discontinuities, they tend to degrade accuracy when mistakenly used 64 in smooth regions of the solution. It is usually difficult to design limiters to achieve both high order accuracy and a non-oscillatory property near discontinuities. Qiu and Shu [52] and Zhu et al. [76] have made such an attempt using WENO methodol- ogy [34, 47, 38, 29, 35, 46, 48, 56] as limiters for the DG methods. They use the usual WENO reconstruction based on cell averages of neighboring cells as in [38, 35, 56], to reconstruct the values of the solutions at certain Gaussian quadrature points in the target cells, and then rebuild the solution polynomials from the original cell average and the reconstructed values at the Gaussian quadrature points through numerical integration of the moments. This limiter needs to use the information from not only the immediate neighboring cells but also neighbors’ neighbors, making it complicated to implement in multi-dimensions, especially for unstructured meshes [76, 35, 73]. The effort in [49, 51] attempts to construct Hermite type WENO approximations, which use the information of not only the cell averages but also the lower order mo- ments such as slopes, to reduce the spread of reconstruction stencils. However for higher order methods the information of neighbors’ neighbors is still needed. In this chapter, we use the WENO methodology to design a new and simpler limiter for the RKDG methods. We do not reconstruct the point values or moments individually and separately, but attempt to reconstruct the entire polynomial in one shot, using the information only from the target cell and its immediate neighbors. As we will see later, this approach simultaneously removes the problem of negative weights and reduces considerably the complexity of implementation, especially for multi-dimensional meshes including unstructured meshes. Comparing with previous WENO type limiters for DG schemes, the limiter developed in this chapter uses most fully the information of the complete polynomials which are already available for DG methods in the target and neighboring cells. Because of this richness of available information, the choice of linear weights are much less restrictive. Essentially, any 65 choice of positive linear weights which add up to one is adequate for accuracy. We refer to [29, 41, 28] and the review paper [61] for a detailed discussion of such practice in choosing linear weights for WENO procedures in finite volume WENO schemes. For our new WENO limiter, following the practice in [28], we give the central cell with a larger linear weight compared to the neighboring cells since in smooth regions the central stencil should provide the most stable reconstruction together with the highest quality in accuracy. 4.2 DG Algorithm Formulation In this section, we give an overview of the algorithm formulation of the RKDG method for solving (4.1) in the one-dimensional case. As discussed in Section 2.1, with a slight abuse of notation, the semi-discrete DG method for solving (4.1) is defined as follows: find the unique function u = u(t) ∈ Vhk such that, for j = 1, · · · , N, Z Z ut v dx − f (u)vx dx + fˆj+ 1 v(x− j+ 1 ) − fˆj− 1 v(x+ j− 1 )=0 (4.2) 2 2 2 2 Ij Ij holds for all test functions v ∈ Vhk . fˆj+ 1 is the the so-called monotone numerical 2 fluxes (approximate or exact Riemann solvers in the system case). To discretize the time variable, we use the TVD third order Runge-Kutta method (2.3). Other TVD, or SSP time discretizations [30] can of course also be used. For simplicity, we consider the forward Euler time discretization of the semi- discrete scheme (4.2). Starting from a solution un ∈ Vhk at time level n (for the initial condition, u0 is taken as the L2 projection of the analytical initial condition 66 into Vhk ). We would like to “limit” it to obtain a new function un,new before advancing it to next time level. That is: find un+1 ∈ Vhk , such that, for j = 1, · · · , N, Z Z un+1 − un,new v dx − f (un,new )vx dx + fˆj+ n,new 1 v(x− j+ 1 ) − fˆj− n,new 1 v(x+ j− 1 )=0 Ij ∆t Ij 2 2 2 2 holds for all test functions v ∈ Vhk . The limiting procedure to go from un to un,new will be discussed in the following section. 4.3 A New WENO Limiter In this section, we present the details of our new WENO limiting procedure for the RKDG methods. As in [52], we also adopt the following framework: 1. Identify the troubled cells, namely, those cells which might need the limiting procedure. 2. Replace the solution polynomials in the troubled cells with reconstructed poly- nomials, which keep the original cell averages, maintain the original high order of accuracy, but are less oscillatory. 4.3.1 Identify the Troubled Cells In this subsection, we discuss the identification of troubled cells. This part is not the emphasis of our chapter and we do not attempt to compare the pros and cons of various troubled cell identification procedures. We will simply use the TVB minmod limiter [24, 23, 21, 20, 25] to identify troubled cells. We emphasize that the goal of our 67 limiter is to be insensitive to the troubled-cell indicators. That is, if more troubled cells are identified than they actually exist, we will increase the computational cost of the algorithm but should otherwise maintain the original high order accuracy. To detect the troubled cells using the limiter described above, we first denote the cell average of the solution u as Z 1 u¯j = u dx (4.3) ∆xj Ij and further denote u˜j = u− j+ 1 − u¯j , u˜˜j = u¯j − u+ j− 1 . (4.4) 2 2 u˜j and u˜˜j are modified either by the usual minmod limiter [33], (mod) (mod)  u˜j = m (˜ uj , ∆+ u¯j , ∆− u¯j ) , u˜˜j = m u˜˜j , ∆+ u¯j , ∆− u¯j , (4.5) where ∆+ u¯j = u¯j+1 − u¯j , ∆− u¯j = u¯j − u¯j−1 , and the minmod function m is defined in (2.9) or by the TVB modified minmod function (2.11). Then, we declare that whenever one of the minmod functions (4.5) gets enacted (returns other than the first argument), this cell is marked as a troubled cell and subject to WENO reconstructions. Of course, if too few cells are identified as trou- bled cells, oscillations and possible instability may not be avoided. If too many cells are identified as troubled cells, the computational cost associated with the second step will increase. Therefore, troubled-cell indicator is a very important issue for 68 WENO limiters. However our main concern in this chapter is how to design the new WENO limiter. We refer the readers to [50] about the comparison among different troubled-cell indicators. We have given the details of identifying troubled cells using the TVB minmod limiters for the one-dimensional scalar case. For two-dimensional scalar case, we use the TVB minmod limiter defined in [20]; for the one- and two-dimensional sys- tems, we use the characteristic-wise TVB minmod limiters defined in [21] and [25], respectively. 4.3.2 Reconstruction of the New Polynomials in the Trou- bled Cells Using a WENO Limiter: Scalar Case In this subsection, we present the details of the reconstruction procedure for the new polynomials in the troubled cells using our WENO limiter for scalar conservation laws. The idea of this WENO limiter is to reconstruct a new polynomial on the troubled cell Ij which is a convex combination of polynomials on this cell and its immediate neighboring cells, with necessary adjustments to keep the original cell av- erage on the target cell. The nonlinear weights in the convex combination coefficients follows the classical WENO procedure. We start with the one-dimensional scalar case. Assume that the cell Ij is a troubled cell. Denote the DG solution polynomial of u on the cells Ij−1 , Ij , Ij+1 as p0 (x), p1 (x), p2 (x), respectively (Figure 4.1). In order to make sure that the reconstructed polynomial maintains the original cell average of p1 in the target cell 69 p1 (x) p0 (x) p2 (x) Ij−1 Ij Ij+1 Figure 4.1: DG polynomials. Ij , we make the following modifications: p˜0 (x) = p0 (x) − p¯0 + p¯1 , p˜2 (x) = p2 (x) − p¯2 + p¯1 , (4.6) where Z Z Z 1 1 1 p¯0 = p0 (x) dx, p¯1 = p1 (x) dx, p¯2 = p2 (x) dx. ∆xj Ij ∆xj Ij ∆xj Ij The final nonlinear WENO reconstruction polynomial pnew 1 (x) is now defined by a convex combination of these modified polynomials: pnew 1 (x) = ω0 p ˜0 (x) + ω1 p1 (x) + ω2 p˜2 (x). (4.7) From (4.6) and (4.7), it is easy to prove that pnew 1 has the same cell average as p1 if the weights satisfy ω0 + ω1 + ω2 = 1. pnew 1 also maintains the order of accuracy of p1 can be directly proved by the following lemma. Lemma 4.3.1. Assume u is sufficiently smooth in [a, b]. Denote a partition of [a, b] in N cells as follows: a = x 1 < x 3 < · · · < xN + 1 = b. 2 2 2 1 Denote Ij = [xj− 1 , xj+ 1 ], xj = 2 (xj− 1 + xj+ 1 ), ∆xj = xj+ 1 − xj− 1 and h = 2 2 2 2 2 2 max1≤j≤N ∆xj . We assume that the mesh is regular, namely, there is a constant 70 c > 0 independent of h such that ∆xj ≥ ch, 1 ≤ j ≤ N. (4.8) If |u(x) − p(x)| ≤ C1 hk+1 , ∀x ∈ Ij where p(x) is a polynomial of degree k, then |u(x) − p(x)| ≤ C2 hk+1 , ∀x ∈ Ij+1 . Here C1 and C2 depend on u and its derivatives but are independent of h. Proof Denote {ξi }ki=0 as the set of the following k + 1 equally spaced points in the cell Ij :   2i − k ξi = xj + ∆xj , i = 0, · · · , k. (4.9) 2(k + 2) Denote uI (x) to be the Lagrange interpolation polynomial of u(x) based on {ξi }. Then k X uI (x) = u(ξi )li (x), i=0 where k Y x − ξi li (x) = . m=0, m6=i ξm − ξi By Theorem B.2.1 in Appendix B.1, we have k u(k+1) (ζ) Y u(x) − uI (x) = (x − ξi ) (4.10) (k + 1)! i=0 71 where min{ξ0 , ξ1 , · · · , ξn , x} ≤ ζ ≤ max{ξ0 , ξ1 , · · · , ξn , x}. On the other hand, p(x) can be written in the following Lagrange form k X p(x) = p(ξi )li (x) i=0 Notice that, by the assumption, we have |u(ξi) − p(ξi )| ≤ C1 hk+1 , i = 0, · · · , k (4.11) Therefore, according to (4.8), (4.9), (4.10) and (4.11), for x ∈ Ij+1 , we have |u(x) − p(x)| ≤ |u(x) − uI (x)| + |uI (x) − p(x)| u(k+1) (ζ) Y k X k = (x − ξi) + (u(ξi) − p(ξi))li (x) (k + 1)! i=0 i=0 (k+1) k u (ζ) X ≤ k+1 k+1 2 h + C1 hk+1 |li (x)| (k + 1)! i=0 (k+1) u k k (ζ) k+1 k+1 X Y x − ξ i ≤ 2 h + C1 hk+1 (k + 1)! ξ − ξi i=0 m=0, m6=i m X 2k hk k |M| k+1 k+1 ≤ 2 h + C1 hk+1  k (k + 1)! i=0 ∆xj k+2 X 2k hk k |M| k+1 k+1 ≤ 2 h + C1 hk+1  (k + 1)! ch k i=0 k+2   2k+1 M C1 (k + 1)(2k + 4)k = + hk+1 (k + 1)! ck where M = maxx∈[a,b] u(k+1) (x). 2k+1 M C1 (k+1)(2k+4)k This proves the lemma with C2 = (k+1)! + ck . 72 Lemma 4.3.2. Under the same assumption as Lemma 4.3.1. If 3 ku(·) − p(·)kL2 (Ij ) ≤ C1 hk+ 2 , where p(x) is a polynomial of degree k, then 3 ku(·) − p(·)kL2 (Ij+1 ) ≤ C2 hk+ 2 . Here C1 and C2 depend on u and its derivatives but are independent of h. Proof Denote {ξi }ki=0 as the set of the shifted k + 1 Gauss points in the cell Ij : zi ξi = xj + ∆xj , i = 0, · · · , k (4.12) 2 where {zi }ki=0 are the roots of the Legendre polynomial of degree k. Here we also denote uI (x) to be the Lagrange interpolation polynomial of u(x) based on {ξi }. Then k X uI (x) = u(ξi )li (x), i=0 and k uk+1 (ζ1 ) Y u(x) − uI (x) = (x − ξi) (4.13) (k + 1)! i=0 where min{ξ0 , ξ1 , · · · , ξn , x} ≤ ζ1 ≤ max{ξ0 , ξ1 , · · · , ξn , x}. We also rewrite p(x) in the following Lagrange form k X p(x) = p(ξi )li (x) i=0 73 Notice that, by the error estimate of Gauss quadrature formula, we have Z ku(·) − p(·)k2L2 (Ij ) = (u(x) − p(x))2 dx Ij k (4.14) X (u2 )(2k+2) (ζ2 )((k + 1)!)4 = (u(ξi ) − p(ξi ))2 ωi + 3 ∆x2k+3 j i=0 (2k + 3)[(2k + 2)!] where ωi is the correponding weights of ξi for Legendre Gauss quadrature rule. ku(·) − p(·)kL2 (Ij+1 ) ≤ ku(·) − uI (·)kL2 (Ij+1 ) + kuI (·) − p(·)kL2 (Ij+1 ) := I + II By (4.13), we have Z 2 I = (u(x) − uI (x))2 dx Ij+1 Z k !2 u (ζ1 ) Y (k+1) ≤ (x − ξi ) dx Ij+1 (k + 1)! i=0 22(k+1) M 2 2k+3 ≤ h ((k + 1)!)2 This gives us 2k+1M k+ 3 I≤ h 2 (4.15) (k + 1)! By (4.8),(4.12) and (4.14), Z 2 II = (uI (x) − p(x))2 dx Ij+1 k X Z 2 = (u(ξi) − p(ξi )) li2 (x)dx i=0 Ij+1 k Z k !2 X Y x − ξi = (u(ξi) − p(ξi ))2 dx i=0 Ij+1 ξm − ξi m=0, m6=i 74 k Z k !2 X Y 2h ≤ (u(ξi) − p(ξi))2 dx i=0 Ij+1 ξm − ξi m=0, m6=i k k !2 X Y 2h = (u(ξi) − p(ξi ))2 zm −zi ∆xj i=0 m=0, m6=i 2 ∆xj k k !2 X Y 4 ≤ (u(ξi) − p(ξi))2 h i=0 m=0, m6=i c(zm − zi ) k k !2 X Y 4 = (u(ξi) − p(ξi ))2 ωi √ h i=0 m=0, m6=i c ωi (zm − zi ) k X 2 ≤M h (u(ξi) − p(ξi))2 ωi i=0   (u2 )(2k+2) (ζ2 )((k + 1)!)4 = M22 h 2 ku(x) − p(x)kL2 (Ij ) − 2k+3 ∆xj (2k + 3)[(2k + 2)!]3  2 (2k+2)  (u ) (ζ 2 )((k + 1)!)4 ≤ M22 h ku(x) − p(x)kL2 (Ij ) + 2 3 ∆xj 2k+3 (2k + 3)[(2k + 2)!]   M3 ((k + 1)!)4 ≤ M22 h C12 h2k+3 + h2k+3 (2k + 3)[(2k + 2)!]3   M22 M3 ((k + 1)!)4 = M22 C12 + h2k+4 (2k + 3)[(2k + 2)!]3 | {z } C32 Qk 4 where M2 = max0≤i≤k √ m=0, m6=i c ωi (zm −zi ) and M3 = maxx∈Ij (u2 )(2k+2) (x). Hence II ≤ C3 hk+2 . (4.16) The conclusion of this lemma can be obtained by (4.15) and (4.16). 75 Following [38, 35, 12, 41], the normalized nonlinear weights are defined as ω ¯j ωj = P , (4.17) lw¯l where the non-normalized nonlinear weights w¯j are functions of the linear weights γj and the so-called smoothness indicators βj as follows: γj ω ¯j = . (4.18) (ε + βj )r We use ε = 10−6 and r = 2 in all computations in this chapter. As in [38, 12], we use the following smoothness indicator for the one-dimensional case: k Z X  2 dl βj = ∆xj2l−1 pj (x) dx. (4.19) l=1 Ij dxl For the two-dimensional case and more details about this smoothness indicator, we refer to [38, 12, 59]. Notice that, because we have used the complete information of the three polynomials p0 (x), p1 (x), p2 (x) in the three cells Ij−1 , Ij , Ij+1 , we do not have extra requirements on the linear weights in order to maintain the original high order accuracy. The linear weights can be chosen to be any set of positive numbers adding up to one. Since for smooth solutions the central cell is usually the best one, we put a larger linear weight on the central cell than on the neighboring cells, i.e. γ1 ≫ γ0 , γ1 ≫ γ2 . γ1 γ1 Lower values of the ratio , γ0 γ2 yield better results on discontinuities while larger values are usually better for smooth solutions. In our numerical tests we take γ0 = 0.001, γ1 = 0.998, γ2 = 0.001, (4.20) 76 which maintains the original high order in smooth regions and can keep essentially non-oscillatory shock transition. We summarize the WENO limiting procedure for one-dimensional scalar con- servation laws as follows. Assuming that DG solution at time level n is un , for j = 1, · · · , N, 1. Use the minmod limiter described in Subsection 4.3.1 to detect whether Ij is a troubled cell or not. 2. If Ij is not a troubled cell, then un,new |Ij = un |Ij . If Ij is a troubled cell, then (a) Denote un on the cells Ij−1 , Ij , Ij+1 as p0 (x), p1 (x), p2 (x), respectively and modify p0 (x), p2 (x) to p˜0 (x), p˜2 (x) using (4.6). (b) Determine the linear weights by (4.20). (c) Compute the smoothness indicators βl using (4.19) for l = 0, 1, 2. (d) Compute the normalized nonlinear weights wl using (4.17) and (4.18) for l = 0, 1, 2. (e) The reconstruction polynomial is given by (4.7), i.e. un,new |Ij = ω0 p˜0 (x) + ω1 p1 (x) + ω2 p˜2 (x). For the two-dimensional scalar case with rectangular meshes considered in this chapter, the limiting procedure is similar as described above. The WENO reconstruc- tion polynomial on the troubled cell Ii,j is a convex combination of the polynomials on this cell and its four neighboring cells {Ii,j−1, Ii,j+1, Ii−1,j , Ii+1,j } suitably mod- ified to maintain the original cell average on the target cell Ii,j . For our numerical tests, we put a larger linear weight 0.996 on the troubled cell Ii,j and the neighboring 77 cells {Ii,j−1, Ii,j+1, Ii−1,j , Ii+1,j } get the smaller linear weight 0.001 (Figure 4.2). The nonlinear weights of the convex combination follow the classical WENO procedure, with smoothness indicators still computed as the sums of L2 -norm squares of all the derivatives of the respective polynomials, as defined in (4.21) X Z βl = (∆xi ∆yj )|α|−1 (D α pl (x, y))2 dx dy (4.21) 1≤|α|≤k Ii,j ∂α where α = (α1 , α2 ), |α| = α1 + α2 and D α = . ∂xα1 ∂y α2 Ii,j+1 0.001 Ii−1,j Ii,j Ii+1,j 0.001 0.996 0.001 Ii,j−1 0.001 Figure 4.2: Linear weights for two dimensional rectangular mesh. 4.3.3 WENO Limiting Procedure for Systems In this subsection, we present the details of the WENO limiting procedure for sys- tems. One-dimensional systems. Consider equation (4.1) where u and f are vectors 78 with m components. In order to achieve better non-oscillatory qualities, the WENO reconstruction limiter is applied with a local characteristic field decompositions, see, e.g. [59] for more details. ∂f Denote the Jacobian matrix by Aj = ∂u u ¯j . Denote the left and right eigenvectors (p) (p) (p) (q) of Aj by lj , rj , p = 1, · · · , m, normalized so that lj · rj = δpq . Let R(¯ uj ) be the m × m matrix with the right eigenvectors as columns, i.e.   (1) (2) (m) R(¯ u j ) = rj , rj , · · · , rj . (4.22) Clearly, R−1 (¯ uj ) is a m × m matrix with the left eigenvectors as rows, that is  T (1) (2) (m) R−1 (¯ uj ) = lj , lj , · · · , lj . (4.23) The WENO limiting procedure for one-dimensional system case is given as fol- lows. Assume that the DG solution at time level n is un (for simplicity, we will use the notation u for the following description), for j = 1, · · · , N, uj ) and R−1 as defined in (4.22) and (4.23). 1. Compute R = R(¯ 2. Compute ∆+ v¯j = R−1 (¯uj+1 − u¯j ), ∆− v¯j = R−1 (¯uj − u¯j−1),     v˜j = R−1 u− j+ 1 − u¯j and v˜˜j = R−1 u¯j − u+j− 1 respectively. 2 2 (mod) (mod) 3. Compute v˜j ˜ vj , ∆+ v¯j , ∆− v¯j ) and v˜˜j = m(˜ ˜ v˜˜j , ∆+ v¯j , ∆− v¯j ) with = m( the modified minmod function m ˜ define in (2.11) for each component of the vectors. 79 4. If (mod) (mod) v˜j = v˜j & v˜˜j = v˜˜j , (4.24) then Ij is not a troubled cell and unew |Ij = u|Ij . Otherwise, (i) Denote the DG polynomial u on the cells Ij+l as uj+l . Project uj+l into the characteristic fields vj+l = R−1 uj+l , for l = −1, 0, 1. Note that vj+l is a m-component vector and each component is a polynomial. (ii) Perform Steps (a), (b), (c), (d), (e) as in the one-dimensional scalar case for each component of vj which is a troubled component, namely the corresponding component makes the condition (4.24) not satisfied. The updated vector vj is denoted as vjnew . (iii) unew |Ij = Rvjnew . Two-dimensional systems. Considering the following two-dimensional system: ut + f (u)x + g(u)y = 0, (4.25) where u, f (u) and g(u) are vectors. The DG scheme with Euler forward time dis- cretization for solving (4.25) on a rectangular mesh is Z un+1 − u v dx dy = Ii,j ∆t Z Z Z f (u)vx dx dy − fˆi+ 1 (y)v(xi+ 1 , y) dy + − fˆi− 1 (y)v(x+ i− 1 , y) dy (4.26) 2 2 2 2 Ii,j Ij Ij Z Z Z − + + g(u)vy dx dy − gˆj+ 1 (x)v(x, yj+ 1 ) dx + gˆi− 1 (x)v(x, yj− 1 ) dx, 2 2 2 2 Ii,j Ii Ii where u = un and fˆ and gˆ are monotone numerical fluxes. 80 For two-dimensional systems, we need to be more careful when using the WENO limiting procedure with a local characteristic field decompositions, since there are two Jacobian matrices corresponding to fluxes in the x and y directions respectively, and therefore two sets of eigenspaces. We perform the characteristic-wise WENO reconstruction in the x-direction and y-direction separately. Assume that Ii,j is a troubled cell, 1. In the x-direction, we choose the polynomials on the cells Ii−1,j , Ii,j , Ii+1,j to reconstruct a new polynomial ux,new i,j using the characteristic-wise WENO lim- ∂f iting procedure with the Jacobian matrix ∂u as in the one-dimensional system case. 2. Similarly, in the y-direction, we choose the polynomials on the cells Ii,j−1, Ii,j , Ii,j+1 to reconstruct a new polynomial uy,new i,j using the characteristic-wise ∂g WENO limiting procedure with the Jacobian matrix ∂u .  3. unew |Ii,j = 1 2 ux,new i,j +uy,new i,j . After limiting the DG solution, we advance it to the next time level for the Euler forward time discretization by Z un+1 − unew v dx dy = Ii,j ∆t Z Z x,new f (u )vx dx dy + g(uy,new )vy dx dy Ii,j Ii,j Z Z − fˆ(ux,new )i+ 1 (y)v(x− i+ 1 , y) dy + fˆ(ux,new )i− 1 (y)v(x+ i− 1 , y) dy 2 2 2 2 Ij Ij Z Z y,new − − gˆ(u )j+ 1 (x)v(x, yj+ 1 ) dx + gˆ(uy,new )i− 1 (x)v(x, yj− + 1 ) dx. 2 2 2 2 Ii Ii 81 The SSP Runge-Kutta time discretization is just a convex combination of such Euler forward steps. 4.4 Numerical Experiments In this section, we provide numerical experiments to demonstrate the performance of the WENO limiter for RKDG methods described in Section 4.3. Even though the advantage of simplicity is most evident for higher-dimensional unstructured meshes, we show in this chapter only results for one and two dimensional structured meshes. For all the computational results, we use the local Lax-Friedrichs flux. For the one-dimensional examples, the CFL number is set to be 0.3 for the P 1 case, 0.15 for the P 2 case and 0.1 for the P 3 case (for the P 3 case ∆t is further reduced in the accuracy test). We have used both uniform and nonuniform meshes in the numerical experiments, obtaining similar results. The nonuniform meshes are obtained from a 20% random perturbation of each node of the uniform mesh. Take the one-dimensional case as an example. The cell boundary point is now xj+ 1 + 20%(rj+ 1 − 0.5)∆x, where xj+ 1 2 2 2 and ∆x are taking values from the uniform mesh and rj+ 1 is a random number from 2 the uniform distribution over the range (0, 1). We will only show results with some accuracy tests on nonuniform meshes as representative tests. For all the accuracy tests (Tables 4.1-4.12), in order to see the effect of the WENO limiter on the accuracy of the RKDG method, we use the TVB minmod limiter with a small TVB constant M = 0.01 to identify troubled cells, resulting in many good cells being identified as troubled cells. For Figures 4.3-4.14, the solid lines are for 82 the exact solutions or grid converged solutions, and the symbols “+” are for the numerical solutions (just one point per cell is plotted). 4.4.1 Scalar Conservation Laws Example 4.4.1. We consider the 1D transport equation: ut + ux = 0, 0 ≤ x ≤ 2π (4.27) with the initial condition u(x, 0) = sin x and periodic boundary conditions. The exact solution of (4.27) is u(x, t) = sin(x − t). Table 4.1 lists the errors and orders of the accuracy for the RKDG method with a WENO limiter compared to the original RKDG method. We can see that the WENO limiter maintains both the designed order of accuracy and the magnitude of the errors of the original RKDG method. Table 4.1: 1D Transport equation. t = 2π. Uniform mesh with N cells DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 3.20E-03 1.73E-02 2.66E-02 7.33E-02 40 7.82E-04 2.03 4.76E-03 1.87 5.05E-03 2.40 2.90E-02 1.34 1 P 80 1.98E-04 1.98 1.24E-03 1.94 4.76E-04 3.41 4.28E-03 2.76 160 5.02E-05 1.98 3.15E-04 1.97 5.33E-05 3.16 4.23E-04 3.34 320 1.27E-05 1.99 7.96E-05 1.99 1.27E-05 2.07 8.27E-05 2.35 20 8.73E-05 5.09E-04 1.21E-04 5.10E-04 40 1.10E-05 2.99 6.44E-05 2.98 1.30E-05 3.21 6.44E-05 2.99 P2 80 1.38E-06 2.99 8.07E-06 3.00 1.49E-06 3.13 8.07E-06 3.00 160 1.73E-07 3.00 1.01E-06 3.00 1.79E-07 3.06 1.01E-06 3.00 320 2.17E-08 3.00 1.26E-07 3.00 2.20E-08 3.03 1.26E-07 3.00 83 Example 4.4.2. We consider the 2D transport equation: ut + ux + uy = 0, 0 ≤ x, y ≤ 2π (4.28) with the initial condition u(x, y, 0) = sin(x + y) and periodic boundary conditions. The exact solution of (4.28) is u(x, y, t) = sin(x + y − 2t). Table 4.2 lists the errors and orders of the accuracy for the RKDG method with a WENO limiter compared to the original RKDG method. Similarly, we can see that the WENO limiter maintains both the designed order of accuracy and the magnitude of accuracy of the original RKDG method. Table 4.2: 2D Transport equation. t = π. Uniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 6.69E-03 4.68E-02 1.76E-02 6.40E-02 40 1.30E-03 2.37 1.34E-02 1.80 3.77E-03 2.22 2.41E-02 1.41 1 P 80 3.09E-04 2.07 3.56E-03 1.91 7.16E-04 2.40 6.94E-03 1.80 160 7.86E-05 1.98 9.15E-04 1.96 8.28E-05 3.11 1.03E-03 2.75 320 2.01E-05 1.97 2.32E-04 1.98 2.01E-05 2.04 2.34E-04 2.15 10 2.79E-03 3.16E-02 2.69E-03 2.99E-02 20 3.40E-04 3.04 4.03E-03 2.97 3.36E-04 3.00 4.05E-03 2.89 40 4.29E-05 2.99 5.07E-04 2.99 4.28E-05 2.97 5.07E-04 3.00 P2 80 5.40E-06 2.99 6.35E-05 3.00 5.40E-06 2.99 6.35E-05 3.00 160 6.79E-07 2.99 7.95E-06 3.00 6.79E-07 2.99 7.95E-06 3.00 320 8.51E-08 3.00 9.95E-07 3.00 8.51E-08 3.00 9.95E-07 3.00 Example 4.4.3. We consider the Burgers equation:   u2 ut + = 0, 0 ≤ x ≤ 2π (4.29) 2 x with the initial condition u(x, 0) = 0.5+sin x and periodic boundary conditions. The exact solution is smooth up to t = 1, then it develops a moving shock which also 84 acts with the rarefaction waves. We can get the exact solution by Newton iteration. For details, see [34]. The errors at t = 0.5 when the solution is smooth are listed in Tables 4.3 and 4.4. We see that the WENO limiter maintains both the designed order of accuracy and the magnitude of the errors of the original RKDG method. In Figure 4.3, we show the RKDG solutions with a WENO limiter at t = 1.5 using 80 cells. We can see that the schemes of all orders perform well in capturing this discontinuity without oscillations. Table 4.3: 1D Burgers equation. t = 0.5. Uniform mesh with N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 3.81E-03 3.82E-02 7.65E-03 6.86E-02 40 9.46E-04 2.01 1.03E-02 1.90 1.81E-03 2.08 2.01E-02 1.77 1 P 80 2.35E-04 2.01 2.65E-03 1.95 3.06E-04 2.57 3.10E-03 2.70 160 5.89E-05 2.00 6.73E-04 1.98 5.93E-05 2.37 6.73E-04 2.20 320 1.47E-05 2.00 1.70E-04 1.98 1.47E-05 2.01 1.70E-04 1.98 20 2.73E-04 5.07E-03 2.68E-04 5.08E-03 40 4.22E-05 2.69 8.96E-04 2.50 4.17E-05 2.68 8.96E-04 2.50 P2 80 6.17E-06 2.77 1.60E-04 2.48 6.17E-06 2.76 1.60E-04 2.48 160 8.86E-07 2.80 2.55E-05 2.65 8.92E-07 2.79 2.55E-05 2.65 320 1.25E-07 2.82 3.79E-06 2.75 1.27E-07 2.81 3.79E-06 2.75 20 1.85E-05 3.64E-04 2.17E-05 3.66E-04 40 1.00E-06 4.20 3.69E-05 3.30 1.03E-06 4.40 3.69E-05 3.31 P3 80 6.08E-08 4.05 2.22E-06 4.06 6.17E-08 4.06 2.22E-06 4.06 160 3.76E-09 4.02 1.44E-07 3.95 3.83E-09 4.01 1.44E-07 3.95 320 2.33E-10 4.01 9.19E-09 3.97 2.45E-10 3.97 9.19E-09 3.97 1.5 +++++ 1.5 +++++ 1.5 ++++++ +++ + +++ + +++ +++ +++ +++ ++ ++ ++ +++ +++ +++ ++ ++ ++ ++ ++ ++ 1 ++ 1 ++ 1 ++ ++ ++ ++ ++ ++ ++ +++ +++ +++ ++ ++ ++ ++ ++ ++ + ++ + ++ ++ ++ ++ ++ 0.5 ++ 0.5 ++ 0.5 ++ u u u +++ +++ + +++ +++ ++ ++ + ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ 0 ++ 0 ++ 0 ++ ++ ++ ++ ++ ++ ++ +++ +++ +++ + + ++ ++ ++ ++ ++ ++ +++ + +++ + +++ +++ +++ +++ +++++++ -0.5 +++++ -0.5 +++++ -0.5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.3: Burgers equation at t = 1.5 with N = 80 cells. 85 Table 4.4: 1D Burgers equation. t = 0.5. Nonuniform mesh with N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 3.84E-03 4.01E-02 6.35E-03 7.39E-02 40 9.85E-04 1.96 1.28E-02 1.65 1.93E-03 1.72 2.07E-02 1.84 1 P 80 2.44E-04 2.02 3.09E-03 2.05 3.16E-04 2.61 3.16E-03 2.71 160 6.14E-05 1.99 8.28E-04 1.90 6.22E-05 2.35 8.28E-04 1.93 320 1.53E-05 2.00 2.36E-04 1.81 1.53E-05 2.02 2.36E-04 1.81 20 2.86E-04 4.87E-03 2.79E-04 4.88E-03 40 4.43E-05 2.69 7.68E-04 2.66 4.38E-05 2.67 7.68E-04 2.67 2 P 80 6.37E-06 2.80 1.61E-04 2.25 6.37E-06 2.78 1.61E-04 2.25 160 9.08E-07 2.81 2.83E-05 2.51 9.16E-07 2.80 2.83E-05 2.51 320 1.30E-07 2.81 4.12E-06 2.78 1.32E-07 2.80 4.12E-06 2.78 20 1.85E-05 3.91E-04 1.85E-05 3.95E-04 40 1.01E-06 4.20 3.20E-05 3.61 1.03E-06 4.17 3.20E-05 3.62 3 P 80 6.42E-08 3.97 2.79E-06 3.52 6.55E-08 3.97 2.79E-06 3.52 160 4.40E-09 3.87 2.80E-07 3.32 4.42E-09 3.89 2.80E-07 3.32 320 2.72E-10 4.02 1.66E-08 4.07 2.77E-10 4.00 1.66E-08 4.07 Example 4.4.4. A two-dimensional version of Example 4.4.3     u2 u2 ut + + = 0, 0 ≤ x, y ≤ 2π (4.30) 2 x 2 y is tested with initial condition u(x, y, 0) = 0.5 + sin(x + y) and periodic boundary conditions. The exact solution is one-dimensional depending only on ξ = x + y; however, our meshes are rectangular in the (x, y) coordinates, and thus this example is a truly two-dimensional test problem. As in Example 4.4.3, we collect the L1 and L∞ errors at t = 0.25 (smooth solution) in Tables 4.5 and 4.6. At t = 0.5, a shock begins to form. We compute the solutions of the RKDG methods using P k polynomials with a WENO limiter with 80 × 80 meshes until t = 0.75 and plot the solution on the diagonal cells in Figures 4.4 and 4.5. Again, we can see that the WENO limiter obtains uniform high order accuracy and sharp, non-oscillatory shock transitions for the RKDG methods. Example 4.4.5. Our last scalar example is the Buckley-Leverett problem that is 86 Table 4.5: 2D Burgers equation at t = 0.25. Uniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 7.56E-03 1.22E-01 1.11E-02 2.07E-01 40 1.92E-03 1.98 3.55E-02 1.78 3.09E-03 1.85 6.57E-02 1.65 P1 80 4.79E-04 2.00 9.36E-03 1.92 6.60E-04 2.23 9.36E-03 2.81 160 1.20E-04 2.00 2.39E-03 1.97 1.23E-04 2.43 2.39E-03 1.97 320 2.99E-05 2.00 6.04E-04 1.98 2.99E-05 2.04 6.04E-04 1.98 20 8.62E-04 4.37E-02 8.62E-04 4.37E-02 40 1.16E-04 2.89 6.08E-03 2.84 1.16E-04 2.89 6.08E-03 2.84 2 P 80 1.50E-05 2.96 9.60E-04 2.66 1.50E-05 2.96 9.60E-04 2.66 160 1.90E-06 2.98 1.35E-04 2.83 1.90E-06 2.98 1.35E-04 2.83 320 2.43E-07 2.97 1.83E-05 2.88 2.44E-07 2.96 1.83E-05 2.88 20 1.30E-04 9.32E-03 1.56E-04 9.03E-03 40 9.28E-06 3.80 7.54E-04 3.63 1.00E-05 3.96 7.54E-04 3.58 3 P 80 5.91E-07 3.97 6.05E-05 3.64 6.21E-07 4.01 6.05E-05 3.64 160 3.73E-08 3.99 3.98E-06 3.92 4.50E-08 3.79 3.98E-06 3.92 320 2.35E-09 3.99 2.56E-07 3.96 2.66E-09 3.83 2.56E-07 3.96 Table 4.6: 2D Burgers equation at t = 0.25. Nonuniform mesh with N × N cells DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 7.70E-03 1.38E-01 1.24E-02 2.82E-01 40 1.99E-03 1.96 4.62E-02 1.58 3.48E-03 1.84 8.75E-02 1.69 1 P 80 4.93E-04 2.01 1.21E-02 1.93 6.86E-04 2.35 1.21E-02 2.85 160 1.24E-04 2.00 3.38E-03 1.84 1.28E-04 2.42 3.38E-03 1.84 320 3.10E-05 1.99 8.50E-04 1.99 3.10E-05 2.05 8.50E-04 1.99 20 8.95E-04 5.48E-02 8.92E-04 5.48E-02 40 1.21E-04 2.89 8.18E-03 2.74 1.21E-04 2.88 8.18E-03 2.74 2 P 80 1.55E-05 2.96 1.37E-03 2.57 1.57E-05 2.95 1.37E-03 2.57 160 1.95E-06 3.00 1.98E-04 2.79 1.99E-06 2.97 1.98E-04 2.79 320 2.50E-07 2.96 2.77E-05 2.84 2.60E-07 2.94 2.77E-05 2.84 20 1.35E-04 1.27E-02 1.60E-04 1.27E-02 40 9.63E-06 3.81 1.08E-03 3.56 1.04E-05 3.95 1.08E-03 3.56 3 P 80 6.18E-07 3.96 1.04E-04 3.37 6.89E-07 3.91 1.04E-04 3.37 160 3.90E-08 3.99 6.99E-06 3.90 6.18E-08 3.48 1.25E-05 3.05 320 2.35E-09 3.99 2.56E-07 3.96 2.66E-09 3.83 2.56E-07 3.96 87 (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.4: 2D Burgers solution at t = 0.75 with 80 × 80 cells. 1.5 ++ ++ 1.5 ++ ++ 1.5 ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + + 1 + + 1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.5 ++ + 0.5 ++ + 0.5 ++ + u u u + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + 0 + + + + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ ++ ++ ++ ++ ++ ++++ ++++ ++++ ++++ ++++ ++++ -0.5 -0.5 -0.5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.5: 2D Burgers solution that cuts along the diagonal at t = 0.75 with 80 × 80 cells. 88 1 ++ ++ 1 +++ 1 +++ + + + + + + + + + + + + + + 0.8 + ++ 0.8 + ++ 0.8 + ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ +++ ++ ++ 0.6 +++ +++ 0.6 +++ +++ 0.6 +++ ++++ ++++ +++ +++ +++ +++ u u u ++++ +++ +++ ++ ++ + 0.4 0.4 + 0.4 + + + + + 0.2 0.2 + 0.2 + + + ++++ +++++ +++++ ++++++ ++++++ ++++++ + ++++++ + +++++ +++++ 0 ++++++++++++++++++++++++++++ ++++++++++++++++ 0 ++++++++++++++++++++++++++++ ++++++++++++++++ 0 ++++++++++++++++++++++++++++ ++++++++++++++++ -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.6: Buckley-Leverett problem at t = 0.4 with 80 cells. governed by the equation   4u2 ut + = 0, (4.31) 4u2 + (1 − u)2 x with the initial condition u = 1 for − 21 ≤ x ≤ 0 and u = 0 elsewhere. The exact solution is a shock-rarefaction-contact discontinuity mixture. The solution is computed up to t = 0.4. Figure 4.6 shows the numerical solutions of RKDG methods with a WENO limiter using N = 80 cells. Again, all schemes perform similarly well for this example. 4.4.2 Euler System in One Dimension The 1D Euler system is given by ut + f (u)x = 0, (4.32) 89 T where u = (ρ, ρv, E)T and f (u) = (ρv, ρv 2 + p, v(E + p)) . Here ρ is the density, v is the velocity, E is the total energy, and p is the pressure, with   1 2 p = (γ − 1) E − ρv ; (4.33) 2 γ = 1.4 is used in the computation. For details of the Jacobian, its eigenvalues, eigenvectors, etc., see [34, 54]. We consider the following typical examples. Example 4.4.6. This example is to test the order of accuracy for the RKDG meth- ods with our WENO limiter. The initial condition is set to be ρ(x, 0) = 1 + 0.2 sin x, v(x, 0) = 1, p(x, 0) = 1 and the boundaries are periodic. The exactly solution is ρ(x, t) = 1 + 0.2 sin(x − t), v(x, t) = 1 and p(x, t) = 1. Tables 4.7 and 4.8 list the errors and orders of accuracy for the RKDG method with a WENO limiter compared to the RKDG method without a limiter. Again, we can see that the WENO limiter maintains both the designed order of accuracy and the magnitude of accuracy of the original RKDG method. Table 4.7: 1D Euler system with initial condition ρ(x, 0) = 1 + 0.2 sin x, v(x, 0) = 1, p(x, 0) = 1 at t = 2π. Uniform mesh with N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 4.68E-04 2.98E-03 6.81E-03 1.60E-02 40 1.23E-04 1.93 7.80E-04 1.93 5.10E-04 3.74 3.10E-03 2.37 P1 80 3.20E-05 1.95 1.99E-04 1.97 3.52E-05 3.86 2.63E-04 3.56 160 8.17E-06 1.97 5.04E-05 1.98 8.20E-06 2.10 5.27E-05 2.32 320 2.07E-06 1.98 1.27E-05 1.99 2.07E-06 1.99 1.27E-05 2.05 20 2.81E-05 1.56E-04 4.45E-05 1.95E-04 40 3.63E-06 2.95 2.04E-05 2.94 4.69E-06 3.25 2.04E-05 3.26 P2 80 4.57E-07 2.99 2.57E-06 2.99 5.25E-07 3.16 2.57E-06 2.99 160 5.72E-08 3.00 3.22E-07 3.00 6.15E-08 3.09 3.22E-07 3.00 320 7.15E-09 3.00 4.02E-08 3.00 7.41E-09 3.05 4.02E-08 3.00 Example 4.4.7. We consider here two well-known problems of the Euler equation 90 Table 4.8: 1D Euler system with initial condition ρ(x, 0) = 1 + 0.2 sin x, v(x, 0) = 1, p(x, 0) = 1 at t = 2π. Nonuniform mesh with N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 4.99E-04 3.88E-03 6.75E-03 1.59E-02 40 1.36E-04 1.88 1.08E-03 1.84 5.57E-04 3.60 3.24E-03 2.30 P1 80 3.49E-05 1.96 2.83E-04 1.94 4.09E-05 3.77 3.23E-04 3.32 160 8.96E-06 1.96 6.81E-05 2.05 9.00E-06 2.18 7.19E-05 2.17 320 2.28E-06 1.98 1.73E-05 1.98 2.28E-06 1.98 1.73E-05 2.06 20 3.05E-05 1.85E-04 5.40E-05 2.35E-04 40 3.95E-06 2.95 2.57E-05 2.85 5.45E-06 3.31 2.57E-05 3.20 2 P 80 4.99E-07 2.98 3.28E-06 2.97 5.06E-07 3.43 3.28E-06 2.97 160 6.20E-08 3.01 4.13E-07 2.99 6.86E-08 2.88 4.13E-07 2.99 320 7.80E-09 2.99 5.42E-08 2.93 8.34E-09 3.04 5.42E-08 2.93 (4.32) which have the following Riemann type initial conditions:    uL , x < 0, u(x, 0) = (4.34)   uR , x > 0. The first one is the Sod problem [64]. The initial data are (ρL , vL , pL ) = (1, 0, 1); (ρR , vR , pR ) = (0.125, 0, 0.1). (4.35) The second one is the Lax problem [43], with the initial data (ρL , vL , pL ) = (0.445, 0.698, 3.528); (ρR , vR , pR ) = (0.5, 0, 0.571). (4.36) The numerical results with the WENO limiter are in Figures 4.7-4.10. To save space, we only show the plots for density. The figures for velocity and pressure are not shown. For Figures 4.7 and 4.9, we use the TVB constant M = 0.01. We can see that there are no oscillations near the discontinuities, however we also observe rather severe smearing, especially for the contact discontinuity, due to this strong limiter. For the Sod problem, relaxing the limiting by taking M = 30 improves the 91 smearing at the price of slight over- and under-shoots, comparing Figure 4.8 with 4.7. The Lax problem is more sensitive to the parameter M, we relax the limiting by taking M = 7, 30, 50 for k = 1, k = 2, k = 3 respectively, see, Figure 4.10. 1 ++++++++++++++++++++++++++ 1 +++++++++++++++++++++++++++ 1 +++++++++++++++++++++++++++ ++ + + + + + + + + + + + + + + 0.8 + + 0.8 + + 0.8 + + + + + + + + + + + + + + + + + ++ ++ ++ density density density 0.6 ++ 0.6 ++ 0.6 ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ +++++++++++++++ +++++++++++++++++ +++++++++++++++++ +++ + + 0.4 + 0.4 + 0.4 + + + + ++ +++++++++++++ ++++++++++++++ +++++++++++++ ++ ++ + + 0.2 + 0.2 + 0.2 + ++++++++++++++ ++++++++++++++ ++++++++++++++ 0 0 0 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.7: Sod problem. t = 2. M = 0.01. N = 100. Density. 1 +++++++++++++++++++++++++++ 1 +++++++++++++++++++++++++++ 1 +++++++++++++++++++++++++++ + + + + + + + + + + + + + + + 0.8 + + 0.8 + + 0.8 + + + + + + + + + + + + + + + + + + + + density density density 0.6 ++ 0.6 ++ 0.6 ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ +++++++++++++++++++ ++++++++++++++++++++ ++++++++++++++++++++ 0.4 + 0.4 0.4 + + + +++++++++++++++ + ++++++++++++++++ ++++++++++++++++ 0.2 0.2 0.2 +++++++++++++++ +++++++++++++++ +++++++++++++++ 0 0 0 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.8: Sod problem. t = 2. M = 30. N = 100. Density. 1.4 1.4 1.4 +++++ ++ + ++ + + +++++ +++++ + +++ + + + + 1.2 + 1.2 + 1.2 + + 1 1 + 1 + + + density density density + 0.8 0.8 0.8 + + + + 0.6 0.6 + 0.6 ++ + ++++++++++++++ ++++++++++++++++ ++++++++++++++++ + ++++++++++++++++++ ++++++++++++++++++ + +++++++++++++++++++ + 0.4 ++++ 0.4 +++++ 0.4 ++++ +++++ +++++ +++++ +++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++ 0.2 0.2 0.2 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.9: Lax problem. t = 1.3. M = 0.01. N = 100. Density. Example 4.4.8. To demonstrate the advantage of higher order methods, we use 92 1.4 1.4 1.4 + ++++++++ + + +++++++++ +++++++++++ + 1.2 + 1.2 1.2 + + 1 1 1 + + density density density + 0.8 0.8 0.8 + + 0.6 + 0.6 + 0.6 + + ++++++++++++++++ ++++++++++++++++ +++++++++++++++++ ++++++++++++++++++ + +++++++++++++++++++ +++++++++++++++++++ 0.4 ++++ 0.4 ++++ 0.4 ++++ +++++ ++++ ++++ +++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++ 0.2 0.2 0.2 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1, M = 7 (b) k = 2, M = 30 (c) k = 3, M = 50 Figure 4.10: Lax problem. t = 2. N = 100. Density. the Euler equation (4.32) with initial condition (ρL , vL , pL ) = (3.857143, 2.629369, 10.333333), when x < −4, (4.37) (ρR , vR , pR ) = (1 + 0.2 sin(5x), 0, 1), when x ≥ −4. This example was used in [63]. It describes the interaction of a Mach 3 shock with a density wave. A Mach 3 shock is initially located at x = −4 and moves to the right. A sine wave is superimposed to the density to the right region to the shock. It contains both shocks and fine structures in smooth regions. Our results are shown in Figures 4.11-4.12. The solid lines are the referenced “exact” solution, which is a converged solution computed by the fifth order finite difference WENO scheme [38] with 2000 grid points. Again, we explore the effect of the TVB constant M in the minmod limiter to identify troubled cells. If M is adjusted adequately, schemes of all orders can perform extremely well, see Figure 4.12. Example 4.4.9. We consider the interaction of blast waves of the Euler equation (4.32) with the initial condition      uL , 0 ≤ x < 0.1,  u(x, 0) = uM , 0.1 ≤ x < 0.9, (4.38)      uR , 0.9 ≤ x < 1, 93 5 5 5 + ++ ++ + + ++ ++ ++ 4 ++++++++++++++++++++++++++++++++++++++++++ ++++++ +++++++ +++++++ ++++++ +++++++ 4 ++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++ +++++++++++ +++++++++++ +++++ +++++ ++ +++++++ ++++++++ +++++++++ ++++ + 4 ++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ ++++++++ +++++++ +++ ++++++ ++ ++++++ +++ ++++++ +++ +++++++ +++++ +++ + + + + +++++ ++++ + ++++++ + +++++ + +++ ++ ++ ++++ +++++ +++++ +++++ + +++++ ++++++++ ++++++++ +++ ++ + + +++++ ++++++ ++++++ +++++ ++ + ++ + + + ++ ++ + + + 3 3 3 + density density density + + + 2 + 2 2 + + ++++++++ ++++ +++ ++++ ++++++ +++++ ++++++++++++++ ++++++++++++ + + +++++++++++++ ++++++++++++ 1 ++++++++ +++++++++ 1 +++ ++++ +++ +++ ++++++ ++++++ 1 +++ +++ +++ +++ ++++++ ++++++ ++ ++ 0 0 0 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.11: The shock density wave interaction problem. t = 1.8. M = 0.01. N = 200. 5 5 5 + + ++ ++ ++ + + ++ ++ + ++ ++++ + ++ + + + +++++ +++ 4 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ +++++++++++ +++++++ ++ +++++++ 4 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ +++++++++++ +++++ + ++++++ 4 +++++++++++++++++++++++++++++++++++++++++++++ +++++++++++ ++++ +++++++ ++++ + ++ + + + + ++ ++++ ++++ +++ + +++ ++ +++ ++++++ +++++ ++++++ +++++ +++ ++++ ++++ ++ + + ++ + + + + + + ++ + +++ ++ ++ + +++ +++ ++++ + + ++++ ++++++ +++++++ + + ++ + + +++ ++ + + + + + ++ + + + + ++ ++ + + ++ + + + + ++ 3 3 3 density density density + + + 2 2 2 + + + +++++++++ ++++++++ +++++++ +++++++ ++++ +++++ ++++ +++++ + +++++++ ++++++++ ++++ +++++ +++++ ++++ ++ 1 +++ ++++ +++++ ++++ ++ 1 1 ++++++++ ++++++++ +++++++++ +++++++++ ++++++++ ++++++++ 0 0 0 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 -5 -3 -1 1 3 5 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.12: The shock density wave interaction problem. t = 1.8. M = 300. N = 200. 94 where ρL = ρM = ρR = 1, vL = vM = vR = 0, pL = 103 , pM = 10−2 , pR = 102 . A reflecting boundary condition is applied to both ends. See [67, 34]. The computed density ρ is plotted at t = 0.038 against the reference “exact” solution, which is a converged solution computed by the fifth order finite difference WENO scheme [38] with 16000 grid points. The results are in Figures 4.13-4.14. For the P 3 case, we add the positivity-preserving limiter [66] to avoid negative density or negative pressure during the time evolution. We can see that the pictures are satisfactory, except for the smearing of contact discontinuities, which seems more serious for this problem. For this problem, we have not noticed any significant difference for taking M from 0.01 to 300. 6 6 6 +++ + ++ + ++++ ++ + +++++ + ++ + + + + ++ + + ++ + + + + ++++ + + + + + +++ + + + + +++ + + ++ + + + + ++ ++ ++++++ + ++++++++ + + ++++ + + ++ + + density ++ + density density 4 + 4 + + + 4 ++++ + + +++ + +++++ + + +++++ +++ ++ + + + + + + + + + + + + ++ + ++++ + 2 +++ 2 ++ + 2 +++ + + + + + + + +++ + ++ + + + + + + + + + + + + + + + + +++++++++++ + + +++++++++++ + + + + + + + +++++++++++ + + + + + + + + + ++ ++ + ++++++++++++++++++++++++++ ++ +++++++++++++++++++++++++++ + + + ++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.13: The blast wave problem. t = 0.038. M = 0.01. N = 400. 6 6 6 ++ +++ + +++ + ++ ++ + + ++++ + + ++ + ++ + + ++ + + + + + ++ + + + +++ + + ++ ++ + ++++++ + + ++ + + + ++ +++ + ++++ + ++++ + +++++ + + + ++ + + +++ + ++ + density density density 4 ++ + + 4 + + + + 4 ++++++ ++++ + ++++ ++++ ++++ + ++ + + + + + + + + + + + + + 2 + +++ 2 ++++ 2 ++++ ++ + ++ +++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +++++++++++ + + +++++++++++ + + + + + ++++++++++ + + + + + + + + + + + + + ++++++++++++++++++++++++++ + + +++++++++++++++++++++++++++ ++ + ++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) k = 1 (b) k = 2 (c) k = 3 Figure 4.14: The blast wave problem. t = 0.038. M = 200. N = 400. 95 4.4.3 Euler System in Two Dimension The 2D Euler system is given by ut + f (u)x + g(u)y = 0,       ρ ρu ρv             (4.39)  ρu   ρu2 + p   ρuv        u= , f (u) =   , g(u) =  .  ρv   ρuv   ρv 2 + p              E u(E + p) v(E + p) Here ρ is the density, (u, v) is the velocity, E is the total energy, and p is the pressure, with   1 2 2 p = (γ − 1) E − ρ(u + v ) ; (4.40) 2 γ = 1.4 is used in the computation. Example 4.4.10. This example is to test the order of accuracy for RKDG methods with our WENO limiter. The initial condition is set to be ρ(x, y, 0) = 1 + 0.2 sin(x + y), u(x, y, 0) = 0.7, v(x, y, 0) = 0.3, p(x, y, 0) = 1 and the boundary conditions are periodic. The exactly solution is ρ(x, y, t) = 1 + 0.2 sin(x + y − t), u(x, y, t) = 0.7, v(x, y, t) = 0.3 and p(x, y, t) = 1. We collect the L1 and L∞ errors at t = 2π in Tables 4.9 and 4.10. Again, we can see that the WENO limiter maintains both the designed order of accuracy and the magnitude of the error of the original RKDG method. Example 4.4.11. We consider the two-dimensional vortex evolution problem [35, 59], which is an idealized problem for the 2D Euler equations (4.39). The set up of this problem is as follows: The mean flow is ρ = 1, p = 1, and (u, v) = (1, 1) 96 Table 4.9: 2D Euler equation with initial condition ρ(x, y, 0) = 1 + 0.2 sin(x + y), u(x, y, 0) = 0.7, v(x, y0) = 0.3, p(x, y, 0) = 1 at t = 2π. Uniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 2.67E-03 6.33E-03 2.68 8.31E-03 2.36E-02 40 3.35E-04 3.00 1.84E-03 1.78 7.63E-04 3.45 3.92E-03 2.59 P1 80 5.92E-05 2.50 5.38E-04 1.77 6.20E-05 3.62 5.97E-04 2.72 160 1.40E-05 2.08 1.44E-04 1.90 1.64E-05 1.92 1.54E-04 1.95 20 8.92E-05 7.41E-04 9.73E-05 7.42E-04 40 1.08E-05 3.04 1.06E-04 2.81 1.14E-05 3.10 1.06E-04 2.81 P2 80 1.29E-06 3.08 1.39E-05 2.93 1.33E-06 3.09 1.39E-05 2.93 160 1.56E-07 3.05 1.76E-06 2.98 1.60E-07 3.06 1.76E-06 2.98 Table 4.10: 2D Euler equation with initial condition ρ(x, y, 0) = 1 + 0.2 sin(x + y), u(x, y, 0) = 0.7, v(x, y0) = 0.3, p(x, y, 0) = 1 at t = 2π. Nonuniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 2.76E-03 7.20E-03 8.14E-03 2.62E-02 40 3.52E-04 2.97 2.70E-03 1.41 8.28E-04 3.30 4.51E-03 2.54 P1 80 6.21E-05 2.50 6.87E-04 1.98 7.01E-05 3.56 7.60E-04 2.57 160 1.49E-05 2.06 2.09E-04 1.72 1.49E-05 2.23 2.12E-04 1.85 20 9.25E-05 8.32E-04 1.08E-04 8.53E-04 40 1.12E-05 3.05 1.27E-04 2.71 1.22E-05 3.14 1.27E-04 2.74 P2 80 1.31E-06 3.09 1.59E-05 3.01 1.40E-06 3.13 1.59E-05 3.01 160 1.58E-007 3.05 2.05E-06 2.95 1.65E-07 3.09 2.05E-06 2.95 97 (diagonal flow). We add, to this mean flow, an isentropic vortex (perturbations in (u, v) and the temperature T = ρp , no perturbation in the entropy S = p ργ : ǫ 0.5(1−r2 ) (γ − 1)ǫ2 1−r2 (δu, δv) = e (−¯ y , x¯), δT = − e , δS = 0, (4.41) 2π 8γπ 2 x, y¯) = (x − 5, y − 5), r 2 = x¯2 + y¯2 , and the vortex strength ǫ = 5. Since the where (¯ mean flow is in the diagonal direction, the vortex movement is not aligned with the mesh direction. The computational domain is taken as [−5, 15] × [−5, 15], extended periodically in both directions. It is clear that the exact solution of the Euler equation with the above initial and boundary conditions is just the passive convection of the vortex with the mean velocity. We compute the solution to t = 2 for the accuracy test. The accuracy results are shown in Tables 4.11 and 4.12. Again, we can see that the WENO limiter maintains both the designed order of accuracy and the magnitude of the error of the original RKDG method. Table 4.11: 2D Euler System of Smooth Vortex Evolution at t = 2. Uniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 2.11E-03 1.63E-01 4.05E-03 3.33E-01 40 4.57E-04 2.20 4.04E-02 2.01 8.74E-04 2.21 7.52E-02 2.15 P1 80 8.88E-05 2.36 1.28E-02 1.66 1.37E-04 2.67 1.54E-02 2.29 160 1.67E-05 2.41 2.90E-03 2.14 2.05E-05 2.74 3.24E-03 2.25 20 5.44E-04 1.01E-01 1.95E-03 1.18E-01 40 6.19E-05 3.14 9.99E-03 3.34 8.72E-05 4.48 9.87E-03 3.58 P2 80 7.70E-06 3.01 1.28E-03 2.97 1.24E-05 2.81 1.59E-03 2.63 160 1.08E-06 2.84 1.78E-04 2.85 1.41E-06 3.14 2.21E-04 2.85 Example 4.4.12. We consider the double Mach reflection problem [67]. It contains strong shock waves and contact discontinuity which is a good example to test the numerical scheme to show the ability to capture strong shock wave and the resolution for small scale structure. The computational domain for this problem is chosen to be [0, 4] × [0, 1]. The reflecting wall lies at the bottom, starting from x = 61 . Initially 98 Table 4.12: 2D Euler System of Smooth Vortex Evolution at t = 2. Nonuniform mesh with N × N cells. DG without limiter DG with WENO limiter(M = 0.01) N L1 error order L∞ error order L1 error order L∞ error order 20 2.15E-03 1.67E-01 4.16E-03 3.31E-01 40 4.81E-04 2.16 4.08E-02 2.04 9.29E-04 2.16 7.69E-02 2.10 P1 80 9.89E-05 2.28 1.21E-02 1.75 1.51E-04 2.62 1.80E-02 2.10 160 1.91E-05 2.38 2.79E-03 2.12 2.34E-05 2.69 3.29E-03 2.45 20 5.68E-04 9.37E-02 2.37E-03 1.46E-01 40 6.67E-05 3.09 1.10E-02 3.09 1.03E-04 4.53 1.14E-02 3.68 P2 80 8.09E-06 3.04 1.13E-03 3.28 1.52E-05 2.75 1.63E-03 2.80 160 1.11E-06 2.86 1.67E-04 2.77 1.62E-06 3.23 2.42E-04 2.76 a right-moving Mach 10 shock is positioned at x = 61 , y = 0 and makes a 60◦ angle with the x-axis. For the bottom boundary, the exact post-shock condition is imposed for the part from x = 0 to x = 61 , and a reflective boundary condition is used for the rest. At the top boundary, the flow values are set to describe the exact motion of a Mach 10 shock. We compute the solution up to t = 0.2. As in [67], only results in [0, 3] × [0, 1] are shown. Two different uniform meshes, with 480 × 120 and 960 cells, and three different values of the TVB constant, M = 0.01, M = 100 and M = 200, are used in the numerical experiments. The density is plotted in Figure 4.15 for M = 0.01 and in Figure 4.16 for M = 100. In all the plots, we use 32 contours equally distributed from ρ = 1.35 to 23. It is not easy to observe any significant difference among these results in the picture. However, if we show a “blown-up” portion around the double Mach region, as in Figures 4.17 and 4.18, it is clear that one observes an increased resolution with an increasing k on the same mesh. Also, the resolution is slightly better as M increases from M = 0.01 to M = 200; however, this difference is not significant. Finally, we list in Table 4.13 the percentage of troubled-cells among all the cells. 99 Table 4.13: Percentage of troubled cells in the double Mach problem k=1 k=2 N M = 0.01 M = 100 M = 200 M = 0.01 M = 100 M = 200 480 × 120 20.31 5.65 4.53 25.84 18.78 16.11 960 × 240 14.35 4.69 3.79 23.31 19.62 18.05 1 0.8 0.6 y 0.4 0.2 0 1 2 3 X 1 0.8 0.6 y 0.4 0.2 0 1 2 3 X Figure 4.15: Double Mach reflection problem. M = 0.01. 960 × 240 cells. Thirty two equally spaced density contours from 1.35 to 23. Top: k = 1. Bottom: k = 2. 100 1 0.8 0.6 y 0.4 0.2 0 1 2 3 X 1 0.8 0.6 y 0.4 0.2 0 1 2 3 X Figure 4.16: Double Mach reflection problem. M = 100. 960 × 240 cells. Thirty two equally spaced density contours from 1.35 to 23. Top: k = 1. Bottom: k = 2. 101 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X Figure 4.17: Double Mach reflection problem. 480 × 120 cells. M = 0.01 (top), M = 100 (middle) and M = 200 (bottom). Thirty two equally spaced density contours from 1.35 to 23. Left: k = 1. Right: k = 2. 102 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X 0.6 0.6 0.4 0.4 y y 0.2 0.2 0 0 2 2.2 2.4 2.6 2.8 2 2.2 2.4 2.6 2.8 X X Figure 4.18: Double Mach reflection problem. 960 × 240 cells. M = 0.01 (top), M = 100 (middle) and M = 200 (bottom). Thirty two equally spaced density contours from 1.35 to 23. Left: k = 1. Right: k = 2. Chapter 5 Conclusion 103 104 This dissertation is based on [26, 74, 75] and focuses on the numerical solutions of conservation laws using discontinuous Galerkin methods. It contains three parts. In the first part, we have given an overview of DG methods for solving con- servation laws. Since the basis functions can be completely discontinuous, these methods have flexibility which is not shared by typical finite element methods, such as the allowance of arbitrary triangulation with hanging nodes, complete freedom in changing the polynomial degrees in each element independent of that in the neigh- bor (p adaptivity), and extremely local data structure (elements only communicate with immediate neighbors regardless of the order of accuracy of the scheme) and the resulting embarrassingly high parallel efficiency. In the second part, using Fourier analysis, we have developed quantitative error estimates for semi-discrete and fully discrete DG schemes with various Runge-Kutta methods, and investigated superconvergence properties at Radau points. Based on these quantitative errors, we have computed the necessary number of points per wavelength required to obtain a fixed error for several RKDG schemes. According to this analysis, assume that we know for a given problem the maximum wave num- ber needed to adequately describe the solution, then we can choose the number of points such that the wave is well resolved. This can reduce the cost of computation while preserving the accuracy. If more detailed information is known about the spec- tral distribution of the Fourier coefficients, then this can be used to obtain sharper comparisons of efficiency by weighting the error function appropriately. In the third part, we have designed a new and simpler limiter using WENO methodology for DG methods, with the goal of obtaining a robust limiting procedure to simultaneously obtain uniform high order accuracy and sharp, non-oscillatory shock transitions. The idea of this simple WENO limiter is that the reconstruction 105 polynomial on the troubled cell is a convex combination of the polynomials on this cell and its immediate neighboring cells, with the nonlinear weights of the linear combination following the classical WENO procedure. This limiting procedure can obtain both uniform high order accuracy and sharp, non-oscillatory shock transitions for the RKDG methods. The numerical performance is similar to previous WENO limiters which are much more complicated to implement. Improving the procedure for identifying troubled cells and implementing the limiter for more general meshes and three-dimensional problems constitute ongoing work. Appendix A A.1 Cardano’s Method This method1 is to find the roots of a cubic function in the following form x3 + ax2 + bx + c = 0. (A.1) The three roots are a a a x1 = u1 + u2 − , x2 = ωu1 + ω ¯ u2 − , x2 = ω ¯ u1 + ωu2 − 3 3 3 where s r √ 1 3 3 q q 2 p3 ω=− + i, u1,2 = − ± + 2 2 2 4 27 with a2 2a3 − 9ab p = b− , q =c+ . 3 27a3 1 http://en.wikipedia.org/wiki/Cubic equation 106 107 A.2 Ferrari’s Method This method2 is to find the roots of a quartic function in the following form x4 + ax3 + bx2 + cx + d = 0. (A.2) Its solution can be found by means of the following calculations: 3a2 a3 ac 3a4 a2 b ac α=− + b, β= − + c, γ=− + − + d. 8 8 2 256 16 4 If β = 0, then s p a −α ±t α2 − 4γ x = − ±s (for β = 0 only). (A.3) 4 2 Otherwise, continue with α2 α3 αγ β 2 P = − − γ, Q = − + − , 12 108 3 8 r Q Q2 P 3 R=− ± + , (either sign of the square root will do) 2 4 27 √ 3 U = R, (there are 3 complex roots, any one of them will do)    −5α + U − P 6 3U if U 6= 0 y= √ ,   − 5 α + U − 3Q if U = 0 6 p W = α + 2y. q  2β a ±s W ∓t − 3α + 2y ±s W x=− + . 4 2 2 http://en.wikipedia.org/wiki/Quartic function Appendix B B.1 Lagrange Interpolation polynomial This section1 documents the interpolation polynomial in the Lagrange form. B.1.1 Interpolation Given a set of (n + 1) data points {(x0 , y0 ), (x1 , y1), · · · , (xn , yn )}, the points defined by (xi )0≤i≤n are called points of interpolation. The points (yi )0≤i≤n are the values of interpolation. To interpolate a function f, we define the values of interpolation as follows: yi = f (xi ), ∀i = 0, · · · , n. 1 http://www.math-linux.com/spip.php?article71 108 109 B.2 Lagrange interpolation polynomial The purpose here is to determine the unique polynomial of degree n, Pn which verifies Pn (xi ) = f (xi ), ∀i = 0, . . . , n. The polynomial which meets this equality is Lagrange interpolation polynomial n X Pn (x) = lk (x)f (xk ) k=0 where lk are polynomials of degree n that form a basis of Pn n Y x − xi x − x0 x − xk−1 x − xk+1 x − xn lk (x) = = ··· ··· i=0, i6=k xk − xi xk − x0 xk − xk−1 xk − xk+1 xk − xn B.2.1 Properties of Lagrange interpolation polynomial and Lagrange basis • They are the lk polynomials which verify the following property:    1 i=k lk (xi ) = δki = , ∀i = 0, . . . , n.   0 i= 6 k • They form a basis of the vector space Pn of polynomials of degree at most equal to n. n X αk lk (x) = 0. k=0 110 By setting: x = xi , we obtain: n X n X αk lk (xi ) = αk δki = 0 =⇒ αi = 0. k=0 k=0 The set (lk )0≤k≤n is linearly independent and consists of n + 1 vectors. It is thus a basis of Pn . • Finally, we can easily see that: n X n X Pn (xi ) = lk (xi )f (xk ) = δki f (xk ) = f (xi ). k=0 k=0 B.2.2 Existence and uniqueness of Lagrange interpolation polynomials Existence. The proof is shown above. Actually, it corresponds to the construction of Lagrange interpolation polynomial with regard to the basis (lk )0≤k≤n of Pn . Uniqueness. Consider two elements Pn and Qn of Pn which verify Pn (xi ) = Qn (xi ) = f (xi ), ∀i = 0, . . . , n. Let Rn = (Pn − Qn ) ∈ Pn . The polynomial Rn has (n + 1) roots which are exactly the (xi )0≤i≤n since Rn (xi ) = Pn (xi ) − Qn (xi ) = f (xi ) − f (xi ) = 0, ∀i = 0, . . . , n. Rn has then (n + 1) roots and Rn ∈ Pn . Therefore, Rn = 0 =⇒ Pn = Qn . 111 We have thus shown the existence and uniqueness of Lagrange interpolation polyno- mial. B.2.3 Error in Lagrange interpolation Assume that f ∈ C n+1 ([a, b]) and x ∈ [a, b]. Let I be the closed set defined by I = [min(x, x0 ), max(x, xn )] (the smallest closed set containing x and xi , ∀ i). Theorem B.2.1. n f n+1 (ξ) Y ∀x ∈ [a, b], f (x) − Pn (x) = (x − xi ), (n + 1)! i=0 where min{x0 , x1 , · · · , xn , x} ≤ ξ = ξ(x) ≤ max{x0 , x1 , · · · , xn , x}. Proof: There are two possible ways to prove this theorem: 1) the one encountered in the polynomial interpolation in the form of Newton and, 2) the following proof. If x = xi , the problem is over: f (x) − Pn (x) = f (xi ) − Pn (xi ) = 0. This means f (x)−Pn (x) has n+1 roots. Therefore, we can rewrite it in the following form: n Y f (x) − Pn (x) = Φ(x) (x − xi ). i=0 We also define the function g: n Y g(x, t) = f (t) − Pn (t) − (t − xi )Φ(x). i=0 112 g(x, ·) is (n + 1) times differentiable and evaluates to zero at the (n + 2) points x0 , x1 , . . . , xn , x of the interval I. By successively applying Rolles theorem, g (n+1) (x, ·) evaluates to zero at a point ξ ∈ I: g (n+1) (x, ξ) = 0. The (n + 1)th derivative of g(x, ·) can be easily calculated: g (n+1) (x, t) = f (n+1) (t) − (n + 1)!Φ(x). By setting t = ξ, we have: g (n+1) (x, ξ) = f (n+1) (ξ) − (n + 1)!Φ(x) = 0. Therefore f (n+1) (ξ) Φ(x) = . (n + 1)! We finally conclude that n f n+1 (ξ) Y f (x) − Pn (x) = (x − xi ). (n + 1)! i=0 Since ξ is dependent on x, x0 , · · · , xn , ξ is a function of x, i.e. ξ = ξ(x). Furthermore, by Rolle’s theorem, we have min{x0 , x1 , · · · , xn , x} ≤ ξ = ξ(x) ≤ max{x0 , x1 , · · · , xn , x}. 113 Corollary B.2.1. Assume that f ∈ C n+1 ([a, b]) and x ∈ [a, b]. Yn (x − xi ) i=0 ∀x ∈ [a, b], |f (x) − Pn (x)| ≤ sup |f n+1 (x)|. (n + 1)! x∈[a,b] Alternatively stated: (b − a)n+1 ∀x ∈ [a, b], |f (x) − Pn (x)| ≤ sup |f n+1(x)|. (n + 1)! x∈[a,b] Bibliography [1] S. Adjerid, M. Aiffa and J. Flaherty. High-order finite element methods for singularly perturbed elliptic and parabolic problems. SIAM Journal on Applied Mathematics, 55:520–543, 1995. [2] S. Adjerid, K. Devine, J. Flaherty and L. Krivodonova. A posteriori error esti- mation for discontinuous Galerkin solutions of hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 191:1097–1112, 2002. [3] S. Adjerid and A. Klauser. Superconvergence of discontinuous finite element solutions for transient convection-diffusion problems. Journal of Scientific Com- puting, 22–23:5–24, 2005. [4] S. Adjerid and T.C. Massey. A posteriori discontinuous finite element error esti- mation for two-dimensional hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 191:5877–5897, 2002. [5] S. Adjerid and T.C. Massey. Superconvergence of discontinuous Galerkin solu- tions for a nonlinear scalar hyperbolic problem. Computer Methods in Applied Mechanics and Engineering, 195:3331–3346, 2006. [6] S. Adjerid and T. Weinhart. Discontinuous Galerkin error estimation for linear symmetrizable hyperbolic systems. Mathematics of Computation, 80:1335–1367, 2011. 114 115 [7] S. Adjerid and T. Weinhart. Discontinuous Galerkin error estimation for linear symmetric hyperbolic systems. Computer Methods in Applied Mechanics and Engineering, 198:3113–3129, 2009. [8] M. Ainsworth. Dispersive and dissipative behavior of high order discontinuous Galerkin finite element methods. Journal of Computational Physics, 198:106– 130, 2004. [9] M. Ainsworth, P. Monk and W. Muniz. Dispersive and dissipative proper- ties of discontinuous Galerkin finite element methods for the second-order wave equation. Journal of Scientific Computing, 27:5–40, 2006. [10] M. Bakker. A note on C 0 Galerkin methods for two-point boundary problems. Numerische Mathematik, 38:447–452, 1982. [11] M. Bakker. One-dimensional Galerkin methods and superconvergence at inte- rior nodal points. SIAM Journal on Numerical Analysis, 21:101–110, 1984. [12] D.S. Balsara and C.-W. Shu. Monotonicity preserving weighted essentially non- oscillatory schemeswith increasingly high order of accuracy. Journal of Compu- tational Physics, 160:405–452, 2000. [13] R. Biswas, K. Devine and J. Flaherty. Parallel, adaptive finite element methods for conservation laws. Applied Numerical Mathematics, 14:255–283, 1994. [14] A. Burbeau, P. Sagaut and Ch.H. Bruneau. A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods. Journal of Computa- tional Physics, 169:111–150, 2001. [15] P. Castillo. A superconvergence result for discontinuous Galerkin methods ap- plied to elliptic problems. Computer Methods in Applied Mechanics and Engi- neering, 192:4675–4685, 2003. 116 [16] F. Celiker and B. Cockburn. Superconvergence of the numerical traces for dis- continuous Galerkin and hybridized methods for convection-diffusion problems in one space dimension. Mathematics of Computation, 76:67–96, 2007. [17] Y. Cheng and C.-W. Shu. Superconvergence and time evolution of discontinuous Galerkin finite element solutions. Journal of Computational Physics, 227:9612– 9627, 2008. [18] Y. Cheng and C.-W. Shu. Superconvergence of local discontinuous Galerkin methods for one-dimensional convection-diffusion equations. Computers and Structures, 87:630–641, 2009. [19] Y. Cheng and C.-W. Shu. Superconvergence of discontinuous Galerkin and local discontinuous Galerkin schemes for linear hyperbolic and convection-diffusion equations in one space dimension. SIAM Journal on Numerical Analysis, 47:4044–4072, 2010. [20] B. Cockburn, S. Hou and C.-W. Shu. The Runge-Kutta local projection dis- continuous Galerkin finite element method for conservation laws IV: The mul- tidimensional case. Mathematics of Computation, 54:545–581, 1990. [21] B. Cockburn, S.-Y. Lin and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One- dimensional systems. Journal of Computational Physics, 84:90–113, 1989. [22] B. Cockburn, M. Luskin, C.-W. Shu and E. S¨ uli. Enhanced accuracy by post- processing for finite element methods for hyperbolic equations. Mathematics of Computation, 72:577–606, 2003. [23] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework. Mathematics of Computation, 52:411–435, 1989. 117 [24] B. Cockburn and C.-W. Shu. The Runge-Kutta local projection P 1 discontinu- ous Galerkin method for scalar conservation laws. Mathematical Modelling and Numerical Analysis, 25:337–361, 1991. [25] B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. Journal of Computational Physics, 141:199–224, 1998. [26] B. Cockburn and C.-W. Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16:173–261, 2001. [27] J. Douglas and T. Dupont. Some superconvergence results for Galerkin methods for the approximate solution of two-point boundary value problems. In J.J. Miller, editor, Topics in Numerical Analysis II, Academic Press, 1973, 89–93. [28] M. Dumbser and M. K¨aser. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems. Journal of Com- putational Physics, 221:693–723, 2007. [29] O. Friedrich. Weighted essentially non-oscillatory schemes for the interpola- tion of mean values on unstructured grids. Journal of Computational Physics, 144:194–212, 1998. [30] S. Gottlieb, D.I. Ketcheson and C.-W. Shu. High order strong stability preserv- ing time discretizations. Journal of Scientific Computing, 38:251–289, 2009. [31] S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes. Mathematics of Computation, 67:73–85, 1998. [32] S. Gottlieb, C.-W. Shu and E. Tadmor. Strong stability preserving high order time discretization methods. SIAM Review, 43:89–112, 2001. 118 [33] A. Harten. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics, 49:257–393, 1983. [34] A. Harten, B. Engquist, S. Osher and S.R. Chakravarthy. Uniformly high or- der accurate essentially non-oscillatory schemes, III. Journal of Computational Physics, 71:231–303, 1987. [35] C. Hu and C.-W. Shu. Weighted essentially non-oscillatory schemes on trian- gular meshes. Journal of Computational Physics, 150:97–127, 1999. [36] F.Q. Hu and H.L. Atkins. Eigensolution analysis of the discontinuous Galerkin method with nonuniform grids: I. one space dimension. Journal of Computa- tional Physics, 182:516–545, 2002. [37] F.Q. Hu and H.L. Atkins. Two-dimensional wave analysis of the discontinuous Galerkin method with non-uniform grids and boundaryconditions. In Proceed- ings of the Eighth AIAA/CEAS Aeronautics Conference and Exhibit, Brecken- ridge, Colorado. 2012. AIAA paper no. 2002-2514. [38] G. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes. Journal of Computational Physics, 126:202–228, 1995. [39] C. Johnson and J. Pitk¨aranta. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Mathematics of Computation, 46:1–26, 1986. [40] O. Karakashian and C. Makridakis. A space-time finite element method for the nonlinear Schr¨odinger equation: the discontinuous Galerkin method. Mathe- matics of Computation, 67:479–499, 1998. [41] M. K¨aser and A. Iske. ADER schemes on adaptive triangular meshes for scalar conservation laws. Journal of Computational Physics, 205:486–508, 2005. 119 [42] H.-O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24:199–215, 1972. [43] P. D. Lax. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Communications on Pure and Applied Mathematics, 7:159–193, 1954. [44] P. Lesaint and P.A. Raviart. On a finite element method for solving the neu- tron transport equation. In C. de Boor, editor, Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, 1974, 89–145. [45] R.J. LeVeque. Numerical Methods for Conservation Laws. Birkhauser Verlag, Basel, 1990. [46] D. Levy, G. Puppo and G. Russo. Central WENO Schemes for hyperbolic systems of conservation laws. Mathematical Modelling and Numerical Analysis (M 2 AN), 33:547–571, 2001. [47] X. Liu, S. Osher and T. Chen. Weighted essentially non-oscillatory schemes. Journal of Computational Physics, 115:200–212, 1994. [48] J. Qiu and C.-W. Shu. On the construction, comparison, and local characteristic decomposition for high-order central WENO schemes. Journal of Computational Physics, 183:187–209, 2002. [49] J. Qiu and C.-W. Shu. Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: one-dimensional case. Journal of Computational Physics, 193:115–135, 2003. [50] J. Qiu and C.-W. Shu. A comparison of troubled-cell indicators for Runge- Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters. SIAM Journal on Scientific Computing, 27:995–1013, 2005. 120 [51] J. Qiu and C.-W. Shu. Hermite WENO schemes and their application as lim- iters for Runge-Kutta discontinuous Galerkin method II: two dimensional case. Computers and Fluids, 34:642–663, 2005. [52] J. Qiu and C.-W. Shu. Runge-Kutta discontinuous Galerkin method using WENO limiters. SIAM Journal on Scientific Computing, 26:907–929, 2005. [53] W.H. Reed and T.R. Hill. Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, Los Alamos, 1973. [54] P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 135:250–258, 1997. [55] S.J. Sherwin. Dispersion analysis of the continuous and discontinuous Galerkin formulations. In B. Cockburn, G. Karniadakis, and C.-W. Shu, editors, Discon- tinuous Galerkin methods: Theory, Computational and Applications, Springer, 2000, 425–431. [56] J. Shi, C. Hu and C.-W. Shu. A technique of treating negative weights in WENO schemes. Journal of Computational Physics, 175:108–127, 2002. [57] C.-W. Shu. TVB uniformly high-order schemes for conservation laws. Mathe- matics of Computation, 49:105–121, 1987. [58] C.-W. Shu. Total-Variation-Diminishing time discretizations. SIAM Journal on Scientific and Statistical Computing, 9:1073–1084, 1988. [59] C.-W. Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced Numerical Approxima- tion of Nonlinear Hyperbolic Equations, B. Cockburn, C. Johnson, C.-W. Shu 121 and E. Tadmor (Editor: A. Quarteroni), Lecture Notes in Mathematics, volume 1697, Springer, 1998, 325-432. [60] C.-W. Shu. A survey of strong stability preserving high order time discretiza- tion methods. In D. Estep and S. Tavener, editors, Collected Lectures on the Preservation of Stability under Discretization, SIAM: Philadelphia, PA, 2002, 51–65. [61] C.-W. Shu. High order weighted essentially non-oscillatory schemes for convec- tion dominated problems. SIAM Review, 51:82–126, 2009. [62] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77:439–471, 1988. [63] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. Journal of Computational Physics, 83:32–78, 1989. [64] G. Sod. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics, 27:1–31, 1978. [65] B. Swartz and B. Wendroff. The relative efficiency of finite difference and finite element methods. SIAM Journal on Numerical Analysis, 11:979–993, 1974. [66] C. Wang, X. Zhang, C.-W. Shu and J. Ning. Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations. Journal of Compu- tational Physics, 231:653–665, 2012. [67] P. Woodward and P. Colella. The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54:115 – 173, 1984. [68] M. Zhang and C.-W. Shu. An analysis of three different formulations of the discontinuous Galerkin method for diffusion equations. Mathematical Models and Methods in Applied Sciences (M 3 AS), 13:395–413, 2003. 122 [69] M. Zhang and C.-W. Shu. An analysis of and a comparison between the dis- continuous Galerkin and the spectral finite volume methods. Computers and Fluids, 34:581–592, 2005. [70] M. Zhang and C.-W. Shu. Fourier analysis for discontinuous Galerkin and related methods. Chinese Science Bulletin, 54:1809–1816, 2009. [71] Q. Zhang and C.-W. Shu. Error estimates to smooth solutions of Runge–Kutta discontinuous Galerkin methods for scalar conservation laws. SIAM Journal on Numerical Analysis, 42:641–666, 2004. [72] Q. Zhang and C.-W. Shu. Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin method for symmetrizable systems of conservation laws. SIAM Journal on Numerical Analysis, 44:1703–1720, 2006. [73] Y.-T. Zhang and C.-W. Shu. Third order WENO scheme on three dimen- sional tetrahedral meshes. Communications in Computational Physics, 5:836– 848, 2009. [74] X.H. Zhong and C.-W. Shu. Numerical resolution of discontinuous Galerkin methods for time dependent wave equations. Computer Methods in Applied Mechanics and Engineering, 200:2814–2827, 2011. [75] X.H. Zhong and C.-W. Shu. A simple weighted essentially nonoscillatory lim- iter for Runge-Kutta discontinuous Galerkin methods. Submitted to Journal of Computational Physics. [76] J. Zhu, J. Qiu, C.-W. Shu and M. Dumbser. Runge-Kutta discontinuous Galerkin method using WENO limiters II: Unstructured meshes. Journal of Computational Physics, 227:4330–4353, 2008.