High order numerical methods for hyperbolic equations: superconvergence, and applications to δ-singularities and cosmology by Yang Yang B.S., University of Science and Technology of China, Hefei, China 2009 M.S., Brown university, Providence, RI 2011 A Dissertation Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in the Division of Applied Mathematics at Brown University Providence, Rhode Island May 2013 c Copyright 2013 by Yang Yang ° This dissertation by Yang Yang is accepted in its present form by the 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 Jan Hesthaven, Ph.D., Reader Date Johnny Guzm´an, Ph.D., Reader Approved by the Graduate Council Date Peter Weber, Dean of the Graduate School iii The Vita of Yang Yang Yang Yang was born in Tianjin, China, on 03 March 1987, the son of Min Liu and Zhiyong Yang. After completing his work at Tianjin Nankai High School, he went on to the University of Science and Technology of China where he studied mathematics and received his Bachelor of Science in July 2009. After that, he entered the Division of Applied Mathematics at Brown University. Publications: • Y. Yang, I. Roy, C.-W. Shu and L.-Z. Fang, Effect of dust on Lyα photon transfer in optically thick halo , The Astrophysical Journal, 739 (2011), 91(11). • Y. Yang and C.-W. Shu, Analysis of optimal superconvergence of discontinuous Galerkin method for linear hyperbolic equations, SIAM Journal on Numerical Analysis, 50 (2012), 3110-3133. • Y. Yang and C.-W. Shu, Discontinuous Galerkin method for hyperbolic equa- tions involving δ-singularities: negative-order norm error estimates and appli- cations, Numerische Mathematik, to appear. • Y. Yang, I. Roy, C.-W. Shu and L.-Z. Fang, Angular distribution of Lyα reso- nant photons emergent from optically thick medium, submitted. • Y. Yang, D. Wei and C.-W. Shu, Discontinuous Galerkin method for Krause’s consensus models and pressureless Euler equations, submitted. Address: (Updated 22 March 2013) 182 George Street Providence, Rhode Island 02912 (401) 523-8384 This thesis was typed by the author. iv Acknowledgments First of all, I would like to thank my thesis advisor, Professor Chi-Wang Shu, who assigned such an interesting topic to me. Professor Shu is really patient, supportive, available and productive. I can discuss with him almost anytime I wish. Moreover, he also introduced other areas in applied mathematics and computation techniques to me, such as inverse problem and parallel programming, which benefited me during my on-site interview. Second, I want to thank Professor Li-Zhi Fang, who was a professor in physics department at the University of Arizona. Professor Fang gave me plenty of guidance on cosmology. However, very sad news came to me in April 6th 2012, that Professor Fang passed away. I was so astonished, since we finished our last project only one week earlier and started a new one. His enthusiasm in science and research is really impressive to me. Moreover, I also finished this project with Doctor Ishani Roy. With the help of Ishani, I got familiar with the project quickly. Next, I would like to thank my family, who supported me a lot. I understand that they missed me very much during my stay in US, even though they never said so. Moreover, I appreciate my girl friend, who never showed up these years. Therefore, I had nothing to do but concentrated myself on the projects everyday, including holidays, such as Christmas and Chinese new years. I also want to thank my group members, since the whole group formed a big family. Thank Xiong Meng and Doctor Qiang Zhang, who gave me suggestions on the project of superconvergence. Thank Doctor Xiangxiong Zhang, who shared his experience about the positivity-preserving limiter with me. Thank all the other members, who introduced other research areas in the group seminars. Finally, I would like to thank all of the people who have helped me to complete my thesis. v Contents Acknowledgments v I Discontinuous Galerkin methods for hyperbolic equa- tions 1 1 Introduction 2 2 Preliminaries 4 2.1 The DG scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Properties of the finite element space . . . . . . . . . . . . . . . . . . 7 2.4 Properties of the DG spatial discretization . . . . . . . . . . . . . . . 9 2.5 The error equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 Analysis of optimal superconvergence for linear hyperbolic equa- tions 13 3.1 Statement of the main result . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Proof of the main result . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2.1 The initial discretization . . . . . . . . . . . . . . . . . . . . . 17 3.2.2 Step 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.3 Step 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.4 Final estimate . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 vi 3.2.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Numerical tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4 Negative-order norm error estimates for linear hyperbolic equations involving δ-functions 44 4.1 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.1 Convolution kernel . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.2 An approximation result . . . . . . . . . . . . . . . . . . . . . 48 4.2 Singular initial condition . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2.1 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2.2 A proof of Theorem 4.2.1 . . . . . . . . . . . . . . . . . . . . . 50 4.2.3 The negative-order error estimate on Ω1 ⊂⊂ Ω\RT . . . . . . 54 4.3 Singular source term . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3.1 Duhamel’s principle . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3.2 Error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.4.1 Singular initial condition . . . . . . . . . . . . . . . . . . . . . 60 4.4.2 Singular source term . . . . . . . . . . . . . . . . . . . . . . . 64 4.5 Proof of Lemma 4.2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.5.1 The weight function . . . . . . . . . . . . . . . . . . . . . . . 68 4.5.2 The smooth solution . . . . . . . . . . . . . . . . . . . . . . . 70 4.5.3 Error representation and error equations . . . . . . . . . . . . 70 4.5.4 The final estimate . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5 Applications to Krause consensus models and pressureless Euler equations 74 5.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 vii 5.1.1 Limiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.1.2 High order time discretizations . . . . . . . . . . . . . . . . . 79 5.2 Krause’s consensus models . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2.1 Positivity-preserving high order schemes . . . . . . . . . . . . 79 5.2.2 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . 82 5.3 Pressureless Euler equations . . . . . . . . . . . . . . . . . . . . . . . 85 5.3.1 Numerical schemes in one dimension . . . . . . . . . . . . . . 85 5.3.2 Numerical schemes in two dimensions . . . . . . . . . . . . . . 89 5.3.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . 95 5.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 II Numerical cosmology 104 6 WENO solver of transfer equations of resonant photons 107 6.1 Basic theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.1.1 Radiative transfer equation of a dusty halo . . . . . . . . . . . 107 6.1.2 Integrated redistribution function . . . . . . . . . . . . . . . . 110 6.1.3 Eddington approximation . . . . . . . . . . . . . . . . . . . . 111 6.1.4 Re-distribution function . . . . . . . . . . . . . . . . . . . . . 113 6.1.5 Integral of the phase function . . . . . . . . . . . . . . . . . . 115 6.2 Numerical algorithm for equation (6.1) . . . . . . . . . . . . . . . . . 116 6.2.1 Conservation law . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.2.2 The WENO algorithm: approximations to the spatial derivatives117 6.2.3 High order numerical integration . . . . . . . . . . . . . . . . 120 6.2.4 Numerical boundary condition . . . . . . . . . . . . . . . . . . 120 6.2.5 Time Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.3 Numerical algorithm for equations (6.11) and (6.12) . . . . . . . . . . 121 6.3.1 Characteristic decomposition . . . . . . . . . . . . . . . . . . . 121 viii 6.3.2 Numerical Boundary Condition . . . . . . . . . . . . . . . . . 122 7 Effect of dust on Lyα photon transfer in optically thick halo 124 7.1 Basic theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.1.1 Dust models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.1.2 Numerical example: Wouthuysen-Field thermalization . . . . . 125 7.2 Dust effects on photon escape . . . . . . . . . . . . . . . . . . . . . . 127 7.2.1 Model I: scattering of dust . . . . . . . . . . . . . . . . . . . . 127 7.2.2 Model III: absorption of dust . . . . . . . . . . . . . . . . . . 128 7.2.3 Effective absorption optical depth . . . . . . . . . . . . . . . . 129 7.2.4 Escape coefficient . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.3 Dust effects on double-peaked profile . . . . . . . . . . . . . . . . . . 132 7.3.1 Dust and the frequency of double peaks . . . . . . . . . . . . 132 7.3.2 No narrowing and no widening . . . . . . . . . . . . . . . . . 134 7.3.3 Profile of absorption spectrum . . . . . . . . . . . . . . . . . . 135 7.4 Discussions and conclusions . . . . . . . . . . . . . . . . . . . . . . . 137 8 Angular distribution of Lyα resonant photons emergent from opti- cally thick medium 140 8.1 Transfer equations of resonant photons without dust . . . . . . . . . . 141 8.1.1 Test with Field’s analytical solution . . . . . . . . . . . . . . . 141 8.1.2 Time scale of the statistical equilibrium of the angular distri- bution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 8.2 Precision of the Eddington approximation . . . . . . . . . . . . . . . 144 8.2.1 Equations of the Eddington approximation . . . . . . . . . . . 144 8.2.2 Profiles of j and f . . . . . . . . . . . . . . . . . . . . . . . . 145 8.3 Angular distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 8.3.1 Frequency dependence . . . . . . . . . . . . . . . . . . . . . . 146 8.3.2 Dependence of the initial anisotropy . . . . . . . . . . . . . . 148 ix 8.3.3 Collimation of photons of the double peaks . . . . . . . . . . . 150 8.3.4 Large halo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 8.3.5 Effect of anisotropic scattering . . . . . . . . . . . . . . . . . . 153 8.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 x List of Tables 3.1 The error e at the Radau points for (3.34) when using P 1 and P 2 polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 The error ξ for equation (3.34) when using P 1 and P 2 polynomials. . 35 3.3 The cell average of the error e for equation (3.34) when using P 1 and P 2 polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 The error e at the Radau points for (3.35) when using P 1 and P 2 polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 The error ξ for (3.35) when using P 1 and P 2 polynomials. . . . . . . 39 3.6 The cell average of the error e for (3.35) when using P 1 and P 2 poly- nomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.7 Superconvergence results for equation (3.36) when using Q1 and Q2 polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.1 L2 -norm of the error between the numerical solution and the exact solution for (4.28) after post-processing in the region away from the singularity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 L2 -norm of the error between the numerical solution and the exact solution for (4.29) after post-processing in the region away from the singularity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.3 L2 -norm of the error between the numerical solution and the exact solution for (4.31) after post-processing in the region away from the singularity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 xi 4.4 L2 -norm of the error between the numerical solution with and the exact solution for (4.32) before and after post-processing in the region away from the singularity. . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1 L2 -norm of the error between the numerical density and the exact density for initial condition (5.23). . . . . . . . . . . . . . . . . . . . . 96 5.2 L2 -norm of the error between the numerical density and the exact density for initial condition (5.27). . . . . . . . . . . . . . . . . . . . . 100 7.1 Effective absorption optical depth τeffect . . . . . . . . . . . . . . . . . 131 8.1 critical effective optical depth τcrit . . . . . . . . . . . . . . . . . . . . 153 xii List of Figures 3.1 The support of φ: black polygons along the dashed line. . . . . . . . . 25 3.2 The support of P+ φ: the black boxes along the dashed line. . . . . . . 27 3.3 Comparison between |u − v| (solid line) and |¯ v − v˜| (dashed line). . . 37 4.1 Numerical solution for (4.28) at t = 0.5 with (right) and without (left) post-processing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Numerical solution (left) and the cut plot along x = y (right) for (4.29) at t = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.3 Solutions of u (left) and v (right) for (4.30) at t = 0.4. . . . . . . . . 64 4.4 Numerical solutions for (4.31) at t = 0.5 with (right) and without (left) post-processing. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.5 Numerical solutions for (4.32) at t = 0.5 plotted for the cell averages (left), six Gaussian points (middle) and the detailed zoom (right). In the left panel, the solid line is the exact solution and the symbols are the cell averages of the numerical solution. . . . . . . . . . . . . . . . 68 5.1 Numerical density for (5.13) at t = 1000 with N = 400 when using P 0 (left) and P 1 (right) polynomials. . . . . . . . . . . . . . . . . . . 83 5.2 Equal-angle-zoned mesh. . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.3 Numerical density ρ for (5.14) at t = 2000 with N = 200 when using P 0 polynomials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 xiii 5.4 Numerical density (left) and velocity (right) at t = 0.5 with P 1 poly- nomials for initial condition (5.24). . . . . . . . . . . . . . . . . . . . 97 5.5 Numerical density (left) and velocity (right) at t = 0.5 for initial condition (5.25). The solid line shows the exact solution while the symbols show the numerical solution. . . . . . . . . . . . . . . . . . . 98 5.6 Numerical density (left) and velocity (right) at t = 0.5 for initial condition (5.26). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.7 Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.28). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.8 Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.29). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.9 Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.30). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.10 Numerical density (left) and velocity field (right) at t = 0.4 with N = 50 for initial condition (5.31). . . . . . . . . . . . . . . . . . . . 103 7.1 The mean intensity j(η, r, x) at η = 500 and r = 100 for dust models I (left panel), II (middle panel) and III (right panel). The source is √ 2 S0 = 1 and φs (x) = (1/ π)e−x . The parameter a = 10−3 . In each panel, κ is taken to be 0, 10−4 , 10−3 and 10−2 . . . . . . . . . . . . . . 126 7.2 Flux f (η, r, x) of Lyα photons emergent from halos at the boundary R = 102 , and for the dust model I A = 1, g = 0.73. The parameter κ is taken to be 10−4 (left), 10−3 (middle) and 10−2 (right). The source √ 2 is S0 = 1 and φs (x) = (1/ π)e−x . The parameter a = 10−3 . . . . . . 128 7.3 Flux f (η, r, x) of Lyα photons emergent from halos at the boundary r = R = 102 . The parameters of the dust are A = 0 and κ = 10−4 (left), 10−3 (middle) and 10−2 (right). Other parameters are the same as in Figure 7.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 xiv 7.4 The time evolution of the total flux F (η) at the boundary of halos with R = τ0 = 102 (upper panels), and R = τ0 = 104 (lower Panels). √ 2 The source of S0 = 1 and φs (x) = (1/ π)e−x starts to emit photons at time η = 0. The parameters of dust are (A = 1, g = 0.73) (left); (A = 0.32, g = 0.73) (middle) and A = 0 (right). In each panel of R = 102 , κ is taken to be 10−4 , 10−3 and 10−2 . In the cases of R = 104 , κ is taken to be 10−4 , 10−3 . . . . . . . . . . . . . . . . . . . 130 7.5 Escaping coefficient fesc (η) as a function of the optical depth τ0 of halo at time η = 5 × 103 , 104 , and 3.2×104 from bottom to up. Dust is modeled by II, A = 0.32, g = 0.73, and κ = 10−3 . . . . . . . . . . 132 7.6 The two-peak frequencies x+ = |x− | as a function of aτ0 . The pa- rameter a is taken to be 10−2 (left) and 5 × 10−3 (right). Dust model III (pure absorption) is used, and κ is taken to be 10−3 . The dashed straight line gives log x± -log aτ with slope 1/3, which is to show the (aτ )1/3 -law of x± . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.7 ln[f (η, r, x, κ = 0)/f (η, r, x, κ)] as function of x for model II (up), and III (bottom), and κ = 10−3 (left) and 10−2 (right). Other parameters are the same as in Figure 7.2. . . . . . . . . . . . . . . . . . . . . . . 135 7.8 The spectrum of the flux emergent from halo of R = 102 (upper panels) and 104 (lower panels) with central source of equation (7.1) for the dust model I (left), II (middle) and III (right). Other parameters are the same as in Figure 7.2. . . . . . . . . . . . . . . . . . . . . . . 136 8.1 The WENO numerical solutions (solid lines) of equation (8.1) as- 2 suming the sources S is a) S = π −1/2 e−x for all µ (left); b) S = 2 2 10µ4 π −1/2 e−x (middle); and c) S = 14µ6 π −1/2 e−x (right) for µ > 0 and S=0 for µ < 0. The Field’s analytical solution is shown with dot lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xv 8.2 The WENO numerical solutions of angular distributions from equa- 2 tion (8.1) at x = 0, assuming the sources S are a) S = π −1/2 e−x for all 2 2 µ (left) and b) S = 10µ4 π −1/2 e−x (middle); and c) S = 14µ6 π −1/2 e−x (right) for µ > 0 and S=0 for µ < 0. . . . . . . . . . . . . . . . . . . . 144 8.3 A comparison of the solutions j and f with the Eddington approxi- mation (dashed curves) and the solutions of f without the Eddington approximation (solid curves). Relevant parameters are r = R = 102 , and a = 10−3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 8.4 A comparison of the solution s j and f with S=δ(µ−1/2) and δ(µ−1). Relative parameters are r = R = 100, a = 10−3 . . . . . . . . . . . . . 146 8.5 The µ distribution of photons emergent from a halo with radius R = 500. The frequencies are x = 0.0, 0.8, 1.6 and 2.4. The relevant pa- rameters of the calculation are η = 1.2 × 104 , γ = 0, and a = 10−3 . . . 147 8.6 The µ-distribution of halo with radius R =100 and source equation (8.9) with n=1, 2, 4, 6, and 8. For each x, the curves for different n overlap with each other. The frequencies are taken to be x = 0, 0.8, 1.6, 2.4. The parameters of the halos are η = 3500, γ = 0, and a = 10−3 .149 8.7 The µ-distributions at radial positions r = 100 (left), 300 (middle), 500 (right) of a halo with radius R = 1000. The source is given by equation (8.9) and n = 6. The frequencies are taken to be x = 0, 0.8, 1.6, 2.4. The dotted curves in the middle and right panels are µ-distributions at r = 100 x = 2.4. Other relevant parameters are η = 1.2 × 104 , γ = 0 and a = 10−3 . . . . . . . . . . . . . . . . . . . . 151 8.8 The µ-distributions with respect to the effective optical depth. Rele- vant parameters are η = 1.2 × 104 , γ = 0 and a = 10−3 . . . . . . . . . 152 xvi Part I Discontinuous Galerkin methods for hyperbolic equations 1 Chapter 1 Introduction We consider discontinuous Galerkin (DG) methods for solving hyperbolic equations ut + f (u)x = g(x, t), (x, t) ∈ R × (0, T ], (1.1) u(x, 0) = u0 (x), x ∈ R. The DG method is a class of finite element methods using completely discontin- uous piecewise polynomial space to represent as the numerical solutions and as the test functions. The method, first introduced in 1973 by Reed and Hill [86], was gen- eralized by Johnson and Pitk¨aranta to solve scalar linear hyperbolic equations with Lp -norm error estimates [64]. Subsequently, Cockburn et al. studied Runge-Kutta discontinuous Galerkin (RKDG) methods for hyperbolic conservation laws in a series of papers [35, 32, 33, 36]. As mentioned in [26], for linear hyperbolic equations, by using piecewise polynomial of degree k, the DG approximation is (k + 3/2)-th order superconvergent towards a particular projection of the exact solution. However, nu- merical experiments demonstrate that a rate of convergence of k + 2. We will use a dual argument to prove this property. One application of the DG methods is to solve hyperbolic equations involving δ-functions. It is well known that generic solutions of hyperbolic equations are not smooth. Discontinuities or even δ-singularities may appear in the solutions. The 2 3 DG methods have been shown to be L2 -stable for nonlinear hyperbolic equations with L2 -solutions which may contain discontinuities [60, 55]. In our work, we as- sume that the initial condition u0 , or the source term g(x, t), or the solution u(x, t) to (1.1) contains δ-singularities. Such problems appear often in applications, such as pressureless Euler equatons, and are difficult to approximate numerically. Many numerical techniques rely on modifications with smooth kernels which may smear such singularities, leading to large errors in the approximation. On the other hand, the DG methods are based on weak formulations and can be designed directly solve such problems without modifications, leading to very accurate results. We will pro- vide numerical examples to demonstrate this advantage. Moreover, we will give rigorous error estimates for the DG methods on some linear problems involving δ- singularities. As demonstrate that the DG approximations are high order accurate under suitable negative-order norms. This makes possible extraction of supercon- vergence solution by convolving the numerical solutions with suitable kernels. For nonlinear models equations, we will apply boundary-preserving (BP) limiters and prove the L1 -stability of the schemes. Chapter 2 Preliminaries In this chapter, we consider the hyperbolic equation (1.1) on the interval [0, 2π]. 2.1 The DG scheme First, we divide the computational domain Ω = [0, 2π] into N cells 0 = x 1 < x 3 < · · · < xN + 1 = 2π, 2 2 2 and denote ³ ´ Ij = xj− 1 , xj+ 1 2 2 as the cells. hj = xj+ 1 − xj− 1 denotes the length of each cell. We also define 2 2 h = hmax = maxj hj and hmin = minj hj as the lengths of the largest and smallest cells, respectively. We consider regular meshes, that is hmax ≤ Λhmin where Λ ≥ 1 is a constant during mesh refinement. Clearly, if Λ = 1, the mesh is uniform. Next, we define Vh = {v : v|Ij ∈ P k (Ij ), j = 1, · · · , N } 4 5 as the finite element space, where P k (Ij ) denotes the space of polynomials in Ij of degree at most k. We also define Hh1 = {φ : φ|Ij ∈ H 1 (Ij ), ∀j}. The DG scheme we consider is the following: find uh ∈ Vh , such that for any vh ∈ Vh ((uh )t , vh )j = (f (uh ), (vh )x )j − fˆj+ 1 vh− |j+ 1 + fˆj− 1 vh+ |j− 1 + (g(x, t), vh )j , (2.1) 2 2 2 2 R where (w, v)j = Ij wvdx, and vh− |j+ 1 = vh (x− j+ 1 ) denotes the left limit of the function 2 2 vh at xj+ 1 . Likewise for vh+ . Moreover, the numerical flux fˆ is a single valued function 2 defined at the cell interfaces and in general depends on the values of the numerical solution uh from both sides of the interfaces fˆj+ 1 = fˆ(uh (x− j+ 1 ), uh (x+ j+ 1 )). 2 2 2 In general, we use monotone fluxes. For the linear case f (u) = u, we consider the upwind fluxes fˆ = u− h . Then the numerical scheme (2.1) can be written as ((uh )t , vh )j = (uh , (vh )x )j − u− − − + h vh |j+ 1 + uh vh |j− 1 + (g(x, t), vh )j (2.2) 2 2 = −((uh )x , vh )j − [uh ]vh+ |j− 1 + (g(x, t), vh )j , (2.3) 2 where [uh ]j− 1 = uh (x+ j− 1 ) − uh (x− j− 1 ) is the jump of uh across xj− 1 . We use (2.2) 2 2 2 2 and (2.3) in Chapters 3 and 4 for the error estimates. Define Hj (uh , vh ) = (uh , (vh )x )j − u− − − + h vh |j+ 1 + uh vh |j− 1 , (2.4) 2 2 6 such that the DG scheme can be written as ((uh )t , vh )j = Hj (uh , vh ) + (g(x, t), vh )j . If we do not consider the source term (i.e. g(x, t) = 0), the scheme becomes ((uh )t , vh )j = Hj (uh , vh ). (2.5) 2.2 Norms We now define some norms that we use throughout the thesis. Denote kuk0,Ij as the standard L2 -norm of u on cell Ij . For any non-negative natural number `, we also define the norm and seminorm of the Sobolev space H ` (Ij ) as ( X ° ° )1/2 ° ` ° ° dα u °2 °d u° kuk`,Ij = ° ° , |u|`,Ij =° ° ° dxα ° ° dx` ° . 0≤α≤` 0,Ij 0,Ij For convenience, if ` = 0, the corresponding index will be omitted. We also define the L∞ -norm and seminorm by ° α ° ° ` ° °d u° °d u° kuk`,∞,Ij = max ° ° , |u|`,∞,Ij =° ° , 0≤α≤` ° dxα ° ° dx` ° ∞,Ij ∞,Ij where kuk∞,Ij is the standard L∞ -norm of u on cell Ij . Clearly, the L∞ -norm is stronger than the L2 -norm, and in one cell Ij , we have 1/2 kukIj ≤ hj kuk∞,Ij . (2.6) Moreover, we define the norms on D = ∪j∈Γ Ij for some index set Γ as follows: à !1/2 X kuk`,D = kuk2`,Ij , kuk`,∞,D = max kuk`,∞,Ij . j∈Γ j∈Γ 7 For convenience, if D = Ω = [0, 2π], the corresponding index will be omitted. Finally, the negative order Sobolev norm can be defined as R D u(x)φ(x)dx kuk−`,D = sup . φ∈C0∞ (D) kφk`,D 2.3 Properties of the finite element space In this section, we study the basic properties of the finite element space. Let us start with the classical inverse properties. Lemma 2.3.1. Assume u ∈ Vh , then there exists a constant C > 0 independent of h and u such that ° α ° °d u° ° ° −α ° dxα ° ≤ Ch kukD , α ≥ 1, (2.7) D ¯ ¯ ¯´ X ³¯¯ − ¯ ¯ + ¯ ¯ j+ 1 ¯ ¯ j− 1 ¯ ≤ Ch−1/2 kukD , u + u (2.8) 2 2 Ij ∈D where D can be the single cell Ij or the whole computational domain Ω. Proof: The proof is trivial. We just use the fact that the norms in finite dimensional spaces are equivalent. So we skip it here. We define P` (p) as the `-th order L2 -projection of p into Vh , such that (P` (p), u)j = (p, u)j , ∀u ∈ P ` (Ij ). (2.9) In addition, if ` ≥ 1, we can also define two Gauss-Radau projections P+ and P− as: (P+ (p), u)j = (p, u)j , ∀u ∈ P `−1 (Ij ) , and P+ (p)(x+ + j−1/2 ) = p(xj−1/2 ), (2.10) (P− (p), u)j = (p, u)j , ∀u ∈ P `−1 (Ij ) , and P− (p)(x− − j+1/2 ) = p(xj+1/2 ). (2.11) The projections P+ and P− are different from the exact collocation at different end 8 points of each cell. For the projection Ph , which is either Pk , P+ or P− , we denote the error operator by P⊥ h = I − Ph , where I is the identity operator. By a scaling argument, we obtain the following property [28]. Lemma 2.3.2. Suppose the function u(x) ∈ C k+1 (Ij ), then there exists a positive constant C independent of h and u, such that k+1 k+1 kP⊥ h ukIj ≤ Chj |u|k+1,Ij and kP⊥ h uk∞,Ij ≤ Chj |u|∞,k+1,Ij . (2.12) Moreover, one can also prove the following superconvergence property [5]. Lemma 2.3.3. Suppose u(x) ∈ C k+2 (Ij ), and xj is one of the downwind-biased Radau points in the cell Ij , then |(u − P− u)(xj )| ≤ Chk+2 j |u|k+2,∞,Ij . (2.13) However, if u is highly oscillatory or discontinuous, we cannot obtain any useful estimate of kP⊥ h uk by the two lemmas above. Therefore, we consider the following estimate. Lemma 2.3.4. Suppose u(x) is a bounded function, then kPh uk∞,Ij ≤ Ckuk∞,Ij , and kP⊥ h uk∞,Ij ≤ Ckuk∞,Ij . (2.14) Proof: For the simplicity of the presentation, we will only prove it for the P− projection. We consider the projection on the reference cell T = [−1, 1] and define a special norm in P k (T ) as ½ ¯Z ¯ ¾ ¯ 1 ¯ |||v||| = max |v(1)|, ¯¯ p ¯ v(s)s ds¯ : 0 ≤ p ≤ k − 1 . −1 It is not difficult to show this is indeed a norm. Since all norms in P k are equivalent, 9 we have kvk∞,T ≤ C|||v||| for any v ∈ P k (T ). Therefore, for any bounded function u, kP− uk∞,T ≤ C|||P− u||| = C|||u||| ≤ Ckuk∞,T . This proves the assertion on the reference cell. The general case follows from a standard scaling argument. By using (2.6) and Lemma 2.3.4, we obtain 1/2 1/2 kP⊥ ⊥ h ukIj ≤ hj kPh uk∞,Ij ≤ Chj kuk∞,Ij . Now, we consider the projection of functions depending not only on the spatial variable x but also on the time variable t. Suppose u(x, t) is a function differentiable and integrable with respect to t and assume t1 and t2 are two real values such that t1 < t2 . Then we have µZ t2 ¶ Z t2 Ph (ut (x, t)) = (Ph u(x, t))t , and Ph u(x, t)dt = (Ph u(x, t))dt. (2.15) t1 t1 As a result, we do not need to distinguish Ph (ut (x, t)) and (Ph u(x, t))t , and can simply denote them as Ph ut . 2.4 Properties of the DG spatial discretization In this subsection, we present some basic properties about the bilinear form Hj and the L2 -stability condition [30]. We consider the linear case, namely (1.1) with f (u) = u. The following lemma is given by Cockburn [30]. Lemma 2.4.1. Suppose uh is the DG numerical solution which satisfies (2.5) in each cell. By using the upwind flux, we have Z T X 2 kuh (T )k + [uh (t)]2j+1/2 dt ≤ kuh (0)k2 . (2.16) 0 1≤j≤N 10 Lemma 2.4.2. Suppose vh ∈ Vh and q(x) ∈ Hh1 , then the two Gauss-Radau projec- tions satisfy Hj (P⊥ − q(x), vh ) = 0, and Hj (vh , P⊥ + q(x)) = 0. (2.17) Proof: The proof is straight forward. So we skip it here. P P If we define (uh , vh ) = j (uh , vh )j and H(p, q) = j Hj (p, q), then Lemma 2.4.3. Suppose p(x) ∈ Hh1 and vh ∈ Vh , there holds H(P⊥ − p(x), vh ) = 0, and H(vh , P⊥ + p(x)) = 0. (2.18) Proof: The proof follows from Lemma 2.4.2 and the definition of H. So we skip it. 2.5 The error equation In this subsection, we also consider linear equation (i.e. f (u) = u in (1.1)) and proceed to construct the error equations. From (2.2) and definition (2.4), we have for any vh ∈ Vh ((uh )t , vh )j = Hj (uh , vh ) + (g(x, t), vh )j . Clearly, the exact solution u satisfies a similar equation (ut , vh )j = Hj (u, vh ) + (g(x, t), vh )j . Denote the error between the exact solution and the DG numerical solution to be e(t) = u(t) − uh (t). Then we have (et , vh )j = Hj (e, vh ), for any vh ∈ Vh . 11 Following the usual treatment in finite element analysis, we divide the error into the form e(t) = η(t) − ξ(t), where η(t) = u(t) − P− u(t), and ξ(t) = uh (t) − P− u(t). From Lemma 2.4.2, we obtain the error equations of the DG scheme. Suppose vh ∈ Vh then (et , vh )j = −Hj (ξ, vh ) = −(ξ, vhx )j + ξ − vh− |j+ 1 − ξ − vh+ |j− 1 (2.19) 2 2 = (ξx , vh )j + [ξ]vh+ |j− 1 (2.20) 2 because of upwinding. Equations (2.19) and (2.20) are fundamental in our analysis later. Let us finish this section by proving the following lemma. R Lemma 2.5.1. Suppose ξ¯ is the cell average of ξ, that is ξ¯ = ξ¯j = 1 hj Ij ξdx in cell Ij , for any j = 1, · · · , N . Then we have ¯ I ≤ Chj kξx kI ≤ Chj kPk et kI ≤ Chj ket kI . kξ − ξk (2.21) j j j j Proof: The right inequality is trivial and the left one follows from the Poincar´e inequality. So we only need to prove the middle one. Suppose Q is the Legendre polynomial of degree k in [-1,1] and define P = (−1)k Q. Then P satisfies the following three properties: (1) P is uniformly bounded: kP k∞,[−1,1] ≤ 1; (2) P evaluated at the left boundary is 1: P (−1) = 1; R1 (3) P is orthogonal to any polynomials with degree no greater than k−1: −1 P Rdx = 0 for any R(x) ∈ P k−1 ([−1, 1]). 2(x−xj ) Define Pj (x) = P ( hj ), then Pj also satisfies the corresponding three properties 12 in the cell Ij . In (2.20), we take vh = ξx − aPj , where a = ξx+ |j−1/2 is a real number, to obtain kξx k2Ij = (Pk et , ξx − aPj )j ¡ ¢ ≤ kPk et kIj kξx kIj + |a|kPj kIj ³ ´ −1/2 1/2 ≤ kPk et kIj kξx kIj + Chj kξx kIj hj ≤ CkPk et kIj kξx kIj , (2.22) where the constant C does not depend on j, h or u. Here, for the second step we use the Cauchy-Schwarz inequality, for the third one we use (2.6) and (2.8), and the last step is trivial. We finish the proof by dividing both sides of the above equation by kξx kIj . Chapter 3 Analysis of optimal superconvergence for linear hyperbolic equations In this chapter, we study one-dimensional linear hyperbolic conservation laws ut + ux = 0, (x, t) ∈ R × (0, T ], (3.1) u(x, 0) = u0 (x), x ∈ R, where the initial datum u0 is sufficiently smooth. We will consider both the peri- odic boundary condition u(0, t) = u(2π, t) and the initial-boundary value problem with u(0, t) = g(t). We use piecewise k-th degree polynomials to approximate the solution in each cell and prove that, under suitable initial discretization, the rate of convergence for the error between the DG solution and the exact solution is of order (k + 2)-th at the downwind-biased Radau points. Moreover, we also prove order (k + 2)-th superconvergence of the cell averages as well as the error between the DG solution and a particular type of projection of the exact solution. In [116], Zhang and Shu gave explicitly formulae for the DG solution in the case of piecewise linear functions for the linear convection equation on uniform meshes. 13 14 The leading error term is shown to be of a constant magnitude independent of time t. This motivates the division of the numerical error into two parts, one being the leading term and the other one being a superconvergent term. In [5, 6], Adjerid et al. proved the (k + 2)-th order superconvergence of the DG solutions at the downwind-biased Radau points for ordinary differential equa- tions. Later, Adjerid and Weihart [7, 8] investigated the local DG error for multi- dimensional first-order linear symmetric and symmetrizable hyperbolic systems. The authors showed that the projection of the local DG error is also (k + 2)-th order su- perconvergent at the downwind-biased Radau points by performing a local error analysis on Cartesian meshes. The global superconvergence is given by numerical experiments. In [7, 8], only initial-boundary value problems are considered, and the local DG error estimate is valid for t sufficiently large. Subsequently, Adjerid and Baccouch [4] investigated the global convergence of the implicit residual-based a pos- teriori error estimates, and proved that these estimates at a fixed time t converge to the true spatial error in the L2 norm under mesh refinement. In [25], Cheng and Shu proved (k + 32 )-th order superconvergence of the DG solution towards a particular projection of the exact solution. The authors considered the case of piecewise linear polynomials (k = 1) on uniform meshes with periodic boundary conditions for the linear conservation laws. Later Cheng and Shu also proved the same (k + 23 )-th order superconvergence when using piecewise k-th degree polynomials with arbitrary k on arbitrary regular meshes in [26]. In [26] the authors considered both periodic bound- ary conditions and initial-boundary value problems. However, the convergence rate obtained in [26] is not optimal. Numerical tests showed that the error of the DG solution towards the particular projection of the exact solution is (k + 2)-th order accurate, even on highly non-uniform meshes, when a special initial discretization is used. Recently, in [123] Zhong and Shu revisited the same problem and showed that the error between the DG numerical solution and the exact solution is (k + 2)-th order superconvergent at the downwind-biased Radau points and (2k + 1)-th order 15 superconvergent at the downwind point in each cell on uniform meshes with periodic boundary conditions for k = 1, 2 and 3. The proofs in [25, 123] use Fourier analysis and work only for uniform meshes and periodic boundary conditions. Moreover, Fourier analysis is difficult to perform for higher polynomial degree k since it relies explicitly on the structure of the algorithm matrices. In [26], a different framework to prove the superconvergence results that does not rely on Fourier analysis is offered and the results are valid for both periodic boundary conditions and initial-boundary value problems. In this chapter, we improve upon the result in [26]. A new tech- nique is adopted to obtain the optimal rate of superconvergence. The proof works for arbitrary regular meshes and schemes of any order. Even though the proof is given for the simple scalar equation (3.1), the same superconvergent results can be obtained for one-dimensional linear systems using similar points. 3.1 Statement of the main result Before proceeding to the main theorem, we first introduce a special initial discretiza- tion. We wish to require (uh )t = P− (ut ) and kP− u − uh kΩ = O(hk+2 ). (3.2) Notice that the special projection P− is used in the error estimates of the DG methods to derive optimal L2 -error bounds in the literature, e.g., in [118]. As in [26], we will prove that the numerical solution is closer to this special projection of the exact solution than to the exact solution itself. The exact way to discretize the initial data to achieve the property (3.2) will be given in Section 3.2.1. We can now state our main theorem. Theorem 3.1.1. Let u(x, t) ∈ C k+4 be the exact solution of the linear hyperbolic equation (3.1) and uh be the numerical solution of the DG scheme (2.5). The finite element space Vh is made up of piecewise polynomials of degree k ≥ 1 on regular 16 meshes, i.e. the ratio of the length of the largest cell to that of the smallest one is bounded during mesh refinement. At time T there holds the following estimate à N ! 12 1 X |(u − uh )(xj )|2 ≤ C(1 + T 2 )hk+2 kukk+4,∞,Ω , (3.3) N j=1 where Ω is the computational domain, and xj is any one of the downwind-biased Radau points in the cell Ij . The constant C does not depend on h, T or u. Remark 3.1.1. Theorem 3.1.1 is valid for both periodic boundary condition and initial-boundary value problems. Corollary 3.1.1. Suppose the conditions in the above theorem are satisfied, then we have ku − uh kL2 (Ω) ≤ C(1 + T 2 )hk+2 kukk+4,Ω , (3.4) kP− u − uh kL2 (Ω) ≤ C(1 + T 2 )hk+2 kukk+4,Ω , (3.5) where u − uh denotes the cell average of u − uh , and the constant C does not depend on h, T or u. 3.2 Proof of the main result In this section, we first discuss how to discretize the initial datum, then prove the main result, Theorem 3.1.1, and finally briefly discuss the application of the super- convergence results. The proof can be divided into several steps. Briefly, by using the triangle inequality, we separate |(u − uh )(xj )| into two parts, |(u − P− u)(xj )| and |ξ(xj )|. The superconvergence of the first term is given by Lemma 2.3.3 while the second one is more difficult to deal with and we separate this process into two steps. In the first step, we consider the estimates of et and ett . In the second step, we use a quadrature formula and consider the dual problem of (3.1). Besides the main theorem, we also prove Corollary 3.1.1 in this section. 17 3.2.1 The initial discretization In this subsection we consider the suitable discretization of the initial datum. As mentioned in Section 3.1, we would like to have the initial solution satisfy ξt = 0 and kξkΩ ≤ Chk+2 , see (3.2). We start from the requirement ξt = 0 and check whether a special numerical initial solution can be constructed which also satisfies the second requirement kξkΩ ≤ Chk+2 . Let us start from the following lemma. Taking vh = 1 in (2.19), we have R − Lemma 3.2.1. Ij et dx = 0, ∀ 1 ≤ j ≤ N if and only if ξj+ 1 is a constant which 2 does not depend on j. Denote the constant mentioned in the previous lemma as S. Clearly, such a constant gives us freedom to control kξkIj , as is shown in the following lemma. k+3/2 Lemma 3.2.2. Suppose ket kIj ≤ Chj , then S ≤ Chk+2 j if and only if kξkIj ≤ k+5/2 Chj . k+5/2 Proof: Suppose kξkIj ≤ Chj , then by Lemma 2.3.1 we have S ≤ Chk+2 j . On the other hand, suppose S ≤ Chk+2 j , then by Lemma 2.5.1 ξ¯j = ξj+ − ¯ − 1 1 − (ξ − ξj ) j+ 2 2 − ξ¯j kIj −1/2 ≤S+ Chj kξ 1/2 ≤ S + Chj ket kIj ≤ Chk+2 j . Therefore, kξkIj ≤ kξ¯j kIj + kξ − ξ¯j kIj ≤ Chj k+5/2 . k+3/2 Remark 3.2.1. The condition ket kIj ≤ Chj in Lemma 3.2.2 is true because we require ξt = 0. Actually, we can show ket kIj ≤ Chk+1 j |u|k+2,Ij . We will also use this estimate of et later. 18 There is a straightforward corollary of the above lemma. Corollary 3.2.1. Suppose the initial solution satisfies ξt = 0 and S ≤ Chk+2 , then kξkΩ ≤ Chk+2 . Now let us proceed to construct the initial solution uh from ξt = 0. R Lemma 3.2.3. Suppose Ij et = 0, then ξx is uniquely determined by Pk et in the cell Ij . Proof: Let vh+ |j− 1 = 0 in equation (2.20), then we have 2 (Pk et , vh )j = (ξx , vh )j . (3.6) Since the equation is linear, we only need to prove uniqueness. That is, suppose (Pk et , vh )j = 0, ∀ vh ∈ Vh and vh+ |j− 1 = 0, then ξx = 0. To prove this, let p(x) be 2 an arbitrary polynomial of degree no more than k and vh = p − p+ j− 1 . Then 2 (Pk et , p)j = (Pk et , p − p+ ) = 0. j− 1 j 2 This implies that Pk et = 0. By Lemma 2.5.1, we obtain ξx = 0. Remark 3.2.2. The expression of ut can be obtained by the partial differential equa- tion. Therefore it is not difficult to obtain Pk et from ξt = 0. − Now, the only thing left is to determine the value of the constant S = ξj− 1 . By 2 Corollary 3.2.1 we can simply take S = 0. However, such S does not satisfy the conservation of mass. If we consider periodic boundary condition we can select a R special S such that Ω ξ = 0. We first prove that such as S satisfies the property S ≤ Chk+2 . Actually, Z N X N ³ X ´ ξdx = ξ¯j hj = ¯ − 1 hj , S − (ξ − ξ)j+ Ω 2 j=1 j=1 19 which yields N X S|Ω| = ¯ − 1 hj . (ξ − ξ) (3.7) j+ 2 j=1 Then we obtain N à N !1/2 C X 3/2 C X C S≤ ket kIj hj ≤ h2k+5 j |u|k+2,Ω ≤ p hk+2 |u|k+2,Ω . (3.8) |Ω| j=1 |Ω| j=1 |Ω| In the first inequality in (3.8) we use Lemma 2.5.1 and (2.8). For the second inequal- ity we use the Cauchy-Schwartz inequality and the estimate ket kIj ≤ Chk+1 j |u|k+2,Ij P obtained in Remark 3.2.1. The last inequality follows from the fact that hj = |Ω| and hj ≤ h. Now we summarize how to implement the initial discretization. We divide the process into the following steps: (1) Let ξt = 0, then compute the value of et using the PDE. (2) Use Lemma 3.2.3 to find ξx . R (3) Compute ξ − ξ¯ in each cell from ξx and the fact that Ij ¯ (ξ − ξ)dx = 0. (4) Express S by using (3.7) or simply by taking S = 0. (5) Calculate ξ from the expressions of S and ξx . (6) Recover uh = ξ + P− u. From the process mentioned above, we observe that the initial solution is uniquely R − determined by the requirements ξt = 0 and Ω ξdx = 0 or ξj− 1 = 0. 2 3.2.2 Step 1 Now, we proceed to prove Theorem 3.1.1. The estimates of ket kΩ and kett kΩ follow from Lemma 2.3 in [26] with some minor changes. We skip the proofs and state the 20 results in the following two equations: kett (t)kΩ ≤ Chk+1 |u|k+3,Ω + Cthk+1 |u|k+4,Ω , (3.9) ket (t)kΩ ≤ Chk+1 |u|k+2,Ω + Cthk+1 |u|k+3,Ω . (3.10) Therefore, by Lemma 2.5.1 we have ¯ ≤ C(1 + t)hk+2 kukk+3,Ω . kξ − ξk (3.11) Before proceeding to the optimal error estimates of kξkΩ , we use a superconvergence result to prove the optimal error estimate of kek∞,Ω . Following [26], we can easily prove kξ(t)kΩ ≤ C(1 + t)hk+3/2 kukk+3,Ω . Since ξ is a polynomial of degree at most k in each cell, we have kξ(t)k∞,Ij ≤ Ch−1/2 kξ(t)kIj ≤ Ch−1/2 kξ(t)kΩ ≤ C(1 + t)hk+1 kukk+3,Ω . Notice that the right hand side of the above equation does not depend on j. We can therefore take the maximum on both sides to obtain kξ(t)k∞,Ω ≤ C(1 + t)hk+1 kukk+3,Ω . Finally, by Lemma 2.3.2, we obtain ke(t)k∞,Ω ≤ C(1 + t)hk+1 kuk∞,k+3,Ω . (3.12) 21 3.2.3 Step 2 Now, we proceed to estimate e(xj ). By Lemma 2.3.3, only ξ(xj ) is considered. Denote the downwind-biased Radau points of the cell Ij as xij , 0 ≤ i ≤ k. Also denote ψji to be a polynomial of degree k in cell Ij , such that   1 x = xi i ` j ψj (x` ) = . (3.13)  0 x = xij ` 6 2 By the Gauss-Radau quadrature, we have ξ(xij ) = w i hj (ξ, ψji ), where the constant wi is the weight of the quadrature at the ith downwind-biased Radau point on the reference interval [−1, 1]. Therefore, we only need to estimate (ξ, ψji ) for any 0 ≤ i ≤ k. Clearly, kψji k∞ ≤ C, where the positive constant C does not depend on i, j or h. Motivated by [34], we consider the dual problem of (3.1). For convenience, we denote by C a generic positive constant that does not depend on h, T or u, but may depend on Λ. Recall that Λ is the ratio of the length of the largest cell to that of the smallest one. We begin by considering the solution to the dual problem: (1) For the periodic boundary condition, find a function φij such that φij (·, t) satisfies φij t + φij x = 0, (x, t) ∈ R × (0, T ], φij (x, T ) = ψji (x), x ∈ R, (3.14) φij (0, t) = φij (2π, t), t ∈ [0, T ]. (2) For the initial boundary value problem, find a function φij such that φij (·, t) satisfies φij t + φij x = 0, (x, t) ∈ R × (0, T ], φij (x, T ) = ψji (x), x ∈ R, (3.15) φij (2π, t) = 0, t ∈ [0, T ]. For convenience, we drop the subscript j as well as the superscript i and denote ψ 22 as ψji and φ as φij . Following [34] Z T (e(T ), ψ) = (e, φ)(0) + (e, φ)t dt 0 Z T = (e, φ)(0) + [(et , φ) + (e, φt )]dt. (3.16) 0 We apply P+ to deal with the term (et , φ) + (e, φt ), with the definition of the projection given in (2.10). Recalling that e = η − ξ where the notations of ξ and η can be found in Section 2.5, we have (et , φ) + (e, φt ) = (et , P⊥ + φ) + (et , P+ φ) − (e, φx ) = (et , P⊥ + φ) + H(e, P+ φ) − (η, φx ) + (ξ, φx ) = (et , P⊥ + φ) − H(ξ, P+ φ) − (η, φx ) + H(ξ, φ) N X − ξ − [φ]j− 1 + ξ − φ− |N + 1 − ξ − φ+ | 1 2 2 2 j=2 N X = (et , P⊥ + φ) − (η, φx ) − ξ − [φ]j− 1 + ξ − φ− |N + 1 − ξ − φ+ | 1(3.17) , 2 2 2 j=2 where for the last equality we use Lemma 2.4.3. For the periodic boundary condition, the above turns out to be N X (et , φ) + (e, φt ) = (et , P⊥ + φ) − (η, φx ) − ξ − [φ]j− 1 . (3.18) 2 j=1 Integrating in t and noticing the fact that Z T N X ξ − [φ]j− 1 = 0, 2 0 j=1 23 since [φ(t)]i− 1 = 0 except for at most finitely many t, we have 2 Z T Z T Z T [(et , φ) + (e, φt )]dt = (et , P⊥ + φ)dt + (η, φt )dt. (3.19) 0 0 0 For the initial boundary value problem, keeping in mind the fact that ξ −1 = 0, (3.17) 2 becomes N X (et , φ) + (e, φt ) = (et , P⊥ + φ) − (η, φx ) − ξ − [φ]j− 1 + ξ − φ− |N + 1 . (3.20) 2 2 j=2 Integrating the above equation in t, and noticing the fact that Z T N X Z T − ξ [φ] j− 21 = 0, and ξ − φ− |N + 1 = 0, 2 0 j=1 0 since [φ(t)]i− 1 = 0 except for at most finitely many t, and φ− (t)|N + 1 = 0 when 2 2 t < T , we again obtain (3.19). We use integration by parts on the second term of the right hand side of (3.19), Z T Z T (η, φt )dt = (η, ψ)(T ) − (η, φ)(0) − (ηt , φ)dt. (3.21) 0 0 Plugging (3.21) into the second term on the right hand side of (3.19), then plugging (3.19) into the right hand side of (3.16), we obtain Z T Z T (e(T ), ψ) = (e, φ(0)) + (et , P⊥ + φ)dt + (η, ψ)(T ) − (η, φ)(0) − (ηt , φ)dt. 0 0 Noticing that P− u − uh = e − η, we have (P− u − uh , ψ)(T ) = Πj1 + Πj2 + Πj3 , 24 where Πj1 = (P− u − uh , φ)(0), Z T j Π2 = − (P⊥ − ut , φ)dt, 0 Z T Πj3 = (et , P⊥ + φ)dt. 0 For the first term, notice the fact that at t = 0, the support of φj is of length at least hmin . Therefore, each cell contains at most dΛe + 1 such φj , where dΛe denotes the smallest integer no smaller than Λ. In Section 3.2.1 we obtained the estimate kξ(0)kΩ ≤ Chk+2 |u|k+2,Ω , such that N X N X (Πj1 )2 = (P− u − uh , φj )2 (0) j=1 j=1 ≤ Ch(dΛe + 1)kξk2Ω ≤ Ch2k+5 |u|2k+2,Ω . (3.22) The estimate of Πj2 RT We proceed to estimate Πj2 = − 0 (P⊥ − ut , φ). For simplicity, only a periodic bound- ary condition is considered. However the estimate of the initial-boundary value problem can be obtained in exactly the same way. We extend our meshes onto the whole real line periodically, so the domain under consideration in this and the next subsections is R × [0, T ]. Clearly, the characteristic line which passes through (xj− 1 , T ), denoted by lj , is t = x + T − xj− 1 , 0 < t < T . We also assume that lj 2 2 and the cell boundary xi− 1 × [0, T ] intersect at t = tji . Denote the support of φ in 2 R × [0, T ] as Ωj , then the boundaries of the cells separate Ωj into several pieces, as © ª shown in Figure 3.1. Denote Ωji = Ωj ∩Ii ×[0, T ] and kj = min i : Ωji is not empty . Also define ∆j = {kj , kj + 1, · · · , j} to be the index set which contains the subscripts 25 of all the nonempty pieces. Then we can easily realize the following properties: (1) Ωj = ∪i∈∆j Ωji and |∆j | = j − kj + 1 ≤ d ThΛ e + 1. e j = {i ∈ ∆j |Ωj is not a parallelogram}, then |∆ (2) Denote ∆ e j | ≤ dΛe + 2 and i e j. j∈∆ (3) Among those which are not parallelograms, Ωjj is a triangle which lies in the e j = ∪ e Ωj , we have Ω region R × [T − h, T ], and by denoting Ω e j ∈ R × [0, 2h]. i∈∆j \j i (4) Suppose Ωji is a parallelogram then the vertices are (xi− 1 , tji ), (xi+ 1 , tji+1 ), (xi+ 1 , tj+1 i+1 ), 2 2 2 and (xi− 1 , tj+1 i ). 2 t=T Ωjj ti+1j j+1 ti+1 Ωi j j ti j+1 ti Ωk+1 j t=0 Ik Ik+1 Ii Ij Figure 3.1: The support of φ: black polygons along the dashed line. Now we can proceed to obtain the estimate Z T XZ (P⊥ − ut , φ) = P⊥ − ut φ dxdt. 0 i∈∆j Ωji 26 Consider the parallelogram Ωji . And noticing the fact Z à Z ! tji+1 ¡ ¢ tji+1 j+1 j+1 P⊥ − ut (ti ), φ dt = P⊥ − ut (ti ), φ dt tj+1 i tj+1 i à Z ! j+1 = P⊥ − ut (ti ), ψ dx Ij = 0. (3.23) Then we have Z Z tji+1 P⊥ − ut φ dxdy = (P⊥ − ut , φ)dt Ωji tj+1 i Z tji+1 j+1 = (P⊥ ⊥ − ut (t) − P− ut (ti ), φ)dt tj+1 Z i à ÃZ ! ! tji+1 t = P⊥ − utt (τ )dτ , φ dt tj+1 i tj+1 i Z tji+1 Z t ¡ ¢ = P⊥ − utt (τ ), φ dτ dt tj+1 i tj+1 i ≤ Chk+4 |u|k+3,∞,Ω . (3.24) e j . By using the third property of the partition of the Now, we consider Ωjj and Ω support of Ωj , we have Z P⊥ − ut φ dxdt ≤ Ch k+3 |u|k+2,∞,Ω , Ωjj and Z P⊥ − ut φ dxdt ≤ Ch k+3 |u|k+2,∞,Ω . ej Ω Combining the above, we obtain Z T (P⊥ − ut , φ) ≤ Ch k+3 |u|k+2,∞,Ω + CT hk+3 |u|k+3,∞,Ω . 0 27 The estimate of Πj3 We still consider periodic boundary condition and follow the procedure in the pre- vious section. However, there are two differences: (1) The support of P⊥ + φ, denoted as Tj , is different from Ωj . (2) We do not have the local estimates of kett kIj or ket kIj . To deal with the first one, we define Tij = Tj ∩ Ii × (0, T ). Clearly Tij is a rectangle covering Ωji and can be written as Tij = Ii × (t0 , t1 ), where t0 = inf{t : ∃x ∈ Ii s.t. (x, t) ∈ Ωji } and t1 = sup{t : ∃x ∈ Ii s.t. (x, t) ∈ Ωji }. If Ωji is a parallelogram, then t0 = tj+1 i and t1 = tji+1 (see Figure 3.2). We also denote Tej = ∪i∈∆e j \j Tij , then Tj = ∪i∈∆j Tij . Moreover, it is not difficult to obtain Tjj ⊂ Ij × (T − h, T ) and Tej ⊂ R × (0, 2h). Consider Tij such that Ωji is a parallelogram. And notice that t=T j Tj j ti+1 j+1 j ti+1 j Ti t i j+1 ti j Tk j t=0 Tk-1 Ik Ik+1 Ii Ij Figure 3.2: The support of P+ φ: the black boxes along the dashed line. 28 Z à Z ! tji+1 tji+1 (et (tj+1 i ), P⊥ + φ)dt = et (tj+1 i ), P⊥ + φ dt tj+1 i tj+1 i à Z ! = et (tj+1 i ), P⊥ + ψ dx Ij = 0. (3.25) Then we have Z Z tji+1 et P⊥ +φ dxdt = (et , P⊥ + φ)dt Tij tj+1 i Z tji+1 ¡ ¢ = et (t) − et (tj+1 i ), P⊥ + φ dt tj+1 Z i ÃZ ! tji+1 t = ett (τ )dτ, P⊥ + φ dt tj+1 i tj+1 i Z tji+1 3/2 ≤ Ch kett kIi dt. (3.26) tj+1 i Observe that sup{t : (x, t) ∈ Tej } ≤ 2h and inf{t : (x, t) ∈ Tjj } ≥ T − h. Therefore Z Z T 1/2 et P⊥ +φ dxdt ≤ Chj ket kIj dt, Tjj T −h and Z Z 2h X 1/2 et P⊥ +φ dxdt ≤ C ket kIi hi dt Tej 0 e j \j i∈∆  1/2 Z 2h X ≤ Ch1/2  ket k2Ii  dt. (3.27) 0 e j \j i∈∆ Combining the above, we obtain the estimate Πj3 ≤ CΓj1 + CΓj2 + CΓj3 , 29 where X Z tji+1 Γj1 =h 3/2 kett kIi dt, ej tj+1 i i∈∆j \∆ Z T 1/2 Γj2 = hj ket kIj dt, T −h  1/2 Z 2h X Γj3 = h1/2  ket k2Ii  dt. 0 e j \j i∈∆ As mentioned at the beginning, we do not have the local estimates of ket k or kett k, so we need to consider the summation with respect to j. First, we consider Γj1 . Keeping in mind the fact that, for any t ∈ (0, T ) and 1 ≤ i ≤ N , the information of kett (t)kIi is contained in at most dΛe + 1 instances of Γj1 , we have ÃZ !2 N X N X X tji+1 TΛ |Γj1 |2 ≤ h3 d e kett kIi dt j=1 j=1 h tj+1 ej i∈∆j \∆ i N X X Z tji+1 3 ≤ CT h kett k2Ii dt j=1 i∈∆j \∆ ej tj+1 i Z T 3 ≤ CT (dΛe + 1)h kett k2Ω dt 0 Z T ≤ CT h2k+5 (|u|k+3,Ω + t|u|k+4,Ω )2 dt ¡ 0 ¢ ≤ Ch2k+5 T 2 |u|2k+3,Ω + T 4 |u|2k+4,Ω . (3.28) 30 The second term is easy to deal with since N X N X Z T |Γj2 |2 ≤ h2 ket k2Ij dt j=1 j=1 T −h Z T = h2 ket k2Ω dt T −h Z T 2k+4 ≤ Ch (|u|k+2,Ω + t|u|k+3,Ω )2 dt T −h ¡ ¢ ≤ Ch2k+5 |u|2k+2,Ω + T 2 |u|2k+3,Ω dt. (3.29) R 2h The third term is also not difficult. Notice that for fixed i, 0 ket kIi is contained in at most dΛe + 1 instances of Γj3 , we have N X N X Z 2h X |Γj3 |2 ≤ h 2 ket k2Ij dt j=1 j=1 0 ˜ j \j i∈∆ Z 2h 2 ≤ (dΛe + 1)h ket k2Ω dt 0 Z 2h ≤ Ch2k+4 (|u|k+2,Ω + t|u|k+3,Ω )2 dt 2k+5 ¡0 2 ¢ ≤ Ch |u|k+2,Ω + h2 |u|2k+3,Ω dt. (3.30) Combining the above, we obtain N X N X ¡ ¢ |Πj3 |2 ≤C |Γj1 |2 + |Γj2 |2 + |Γj3 |2 ≤ C(1 + T 4 )h2k+5 kuk2k+4,Ω . j=1 j=1 Remark 3.2.3. By the same method mentioned in this subsection, we can also derive that N X |Πj2 |2 ≤ C(1 + T 2 )h2k+5 kuk2k+4,Ω . j=1 Here the upper bound is of T 2 instead of T 4 since kηt kΩ and kηtt kΩ do not grow in time. 31 3.2.4 Final estimate Now we proceed to the final estimate of |ξ(xj )|. We simply sum up all the previous estimates and obtain N X N X ¡ ¢ 2 |(ξ, ψj )| ≤ 3 |Πj1 |2 + |Πj2 |2 + |Πj3 |2 j=1 j=1 ¡ ¢ ≤ Ch2k+5 |u|2k+2,Ω + (1 + T 2 )kuk2k+4,Ω + (1 + T 4 )kuk2k+4,Ω ≤ Ch2k+5 (1 + T 4 )kuk2k+4,Ω . (3.31) Therefore, N N ¯ ¯2 1 X 2 1 X ¯¯ 2 ¯ ¯ |ξ(xj )| = (ξ, ψ ) N j=1 ¯ hj ¯ j N j=1 ≤ Ch2k+4 (1 + T 4 )kuk2k+4,Ω . By Lemma 2.3.3, N 1 X |(u − uh )(xj )|2 ≤ Ch2k+4 (1 + T 4 )kuk2k+4,∞,Ω . (3.32) N j=1 This completes the proof of the main theorem. Notice that, throughout the proof, we have not used any special property of ψ. Hence we can take ψ to be the indicator function of cell Ij , which yields the estimate of (3.4). Finally, (3.5) follows from (3.11) and (3.4). 3.2.5 Applications In Section 3.2.2, we proved optimal error estimates in L∞ by using the superconver- gence of ξ = uh − P− u. This can be considered as an application of the supercon- vergence result. Let us also briefly discuss some other applications. A notice is to use the superconvergence of the cell averages to construct a new a posterior error 32 indicator, in the same spirit as in [67] where the jump sizes at cell interfaces are used as an a posterior error indicator. To expose this, we denote v as the numerical solution instead of uh in this section and consider the cell Ij only. Following the superconvergence result, the cell average of the numerical solution vj is superclose to that of the exact solution uj . If we construct another numerical cell average vej which is not superconvergent, then the difference vej − vj is a good a posterior error indicator of the local error. We can construct vej in the following steps: (1) Extend the polynomial numerical solution from the two neighboring cells, de- noted by vj−1 and vj+1 , to the cell Ij . (2) Compute the cell averages of vj−1 and vj+1 in the cell Ij , and denote them by vg j−1 and vg j+1 respectively. (3) Define vej = θvg j−1 + (1 − θ)vg j+1 , where 0 ≤ θ ≤ 1. In general, the new cell average vej is only (k +1)-th order accurate. (4) The a posterior computable quantity vej − vj is asymptotically equal to the error vej − uj and is therefore a good indicator of the local error. Numerical evidence will be given in Section 3.3. Based on the superconvergence of the downwind-biased Radau points, we can construct another a posterior error indicator. We use vj and uj as the numerical and exact solutions in cell Ij , respectively. Denote xij , 0 ≤ i ≤ k as the downwind-biased Radau points and Sˆ = {xj− 1 , x0j , · · · , xkj }. Then we construct another numerical 2 approximation wj ∈ P k+1 ˆ For convenience, if xˆ ∈ (Ij ) which interpolates vj at x ∈ S. Sˆ is located at the cell interface, v(ˆ x) is denoted as the left limiter of the numerical approximation. Therefore, we have wj (x+ j− 1 ) = vj−1 (x− j− 1 ) and wj (x− j+ 1 ) = vj (x− j+ 1 ). 2 2 2 2 We define w(x) to be such a function such that for any j, w(x) agrees with wj in cell Ij . Obviously, w(x) is a continuous function. In (3.12), we mentioned that the error between the numerical solution and the exact solution is not superconvergent. If we can show that the new numerical approximation wj is superconvergent in the 33 L∞ -norm, the difference vj − wj is a good a posterior error indicator of the local ˆ Clearly, ej ∈ P k+1 (Ij ) be a polynomial thus interpolates uj at x ∈ S. error. Let w ej k∞,Ij ≤ Chk+2 . kuj − w (3.33) ˆ Theorem 3.1.1 yields |vj (ˆ On the other hand, for any xˆ ∈ S, x)| ≤ Chk+3/2 . x) − uj (ˆ Then we have |wj (ˆ x) − w x)| ≤ Chk+3/2 . Define a special norm in P k+1 (Ij ) as ej (ˆ ||v||Sˆ = max {|v(x)|} . x∈Sˆ It is not difficult to show this is indeed a norm. Since all norms in P k+1 are equivalent, we have kvk∞,Ij ≤ C||v||Sˆ for any v ∈ P k+1 (Ij ). Therefore, kwj − w ej ||Sˆ ≤ Chk+3/2 . ej k∞,Ij ≤ C||wj − w By using (3.33), we have kuj − wj k∞,Ij ≤ Chk+3/2 , which further yields kuj − wj k∞ ≤ Chk+3/2 . 3.3 Numerical tests The purpose of this section is to verify our main result, Theorem 3.1.1 as well as Corollary 3.1.1, and to present numerical evidence suggesting that the proved rate of superconvergence is optimal. In most cases, we consider random meshes (that is, each cell boundary point is randomly and independently perturbed from a uniform mesh up to a given percentage) and use Λ to denote the ratio of the length of the largest cell to that of the smallest one. 34 Example 3.1. We solve the following equation    u + ux = 0   t u(x, 0) = esin(x) . (3.34)     u(0, t) = u(2π, t) The exact solution to this problem is u(x, t) = esin(x−t) . We use ninth order SSP Runge-Kutta discretization in time [46] and take ∆t = 0.05hmin to reduce the time error. Non-uniform meshes are obtained by randomly and independently perturbing each node in a uniform mesh by up to 40%, and the example is tested with both P 1 and P 2 polynomials. The error in Theorem 3.1.1 at different downwind-biased Radau points at t = 1 on random meshes of N cells are computed. In Table 3.1, we observe (2k + 1)-th order superconvergence at the downwind point and (k + 2)-th order superconvergence at other Radau points. The initial solution is obtained by exactly the same way as mentioned in Section 3.2.1. The downwind-biased Radau points on the interval [-1,1] are − 13 and 1 for P 1 √ √ −1− 6 −1+ 6 polynomials, and are 5 , 5 and 1 for P 2 ones. Table 3.2 shows the rate of convergence of the error ξ. We observe that the order is k + 2, indicating that the estimate in (3.5) is sharp. Moreover, we also test the superconvergence for the cell average. Table 3.3 shows the result for Example 3.1 by using the method mentioned in Section 3.2.1 as well as with L2 - and P− -projection for the initial discretization. From the table, we find the convergent rates to be of order 2k + 1, k + 23 and at least k + 2, respectively, for the three different ways of numerical initial discretization. Now we follow the steps in Section 3.2.5 and construct the new numerical cell average. Following the same notations, θ is taken to be θ = 1, i.e. we consider the 35 Table 3.1: The error e at the Radau points for (3.34) when using P 1 and P 2 poly- nomials. 1st Radau point 2nd Radau point downwind point Polynomial N hmax Λ error order error order error order P1 50 0.202 6.056 1.86E-04 - 1.34E-04 - 100 0.111 6.169 3.76E-05 2.64 2.87E-05 2.55 200 5.408e-02 7.339 3.89E-06 3.17 2.92E-06 3.19 400 2.781e-02 6.956 4.90E-07 3.12 3.80E-07 3.07 P2 50 0.202 6.056 1.11E-06 - 1.09E-06 - 2.13E-07 - 100 0.111 6.169 1.70E-07 3.09 1.42E-07 3.37 1.78E-08 4.10 200 5.408e-02 7.339 1.03E-08 3.93 7.94E-09 4.04 4.28E-10 5.21 400 2.781e-02 6.956 7.74E-10 3.89 5.81E-10 3.93 1.37E-11 5.18 Table 3.2: The error ξ for equation (3.34) when using P 1 and P 2 polynomials. L2 norm of ξ P 1 Polynomial P 2 Polynomial N hmax Λ L2 error order L2 error order 50 0.202 6.056 4.09E-04 - 1.85E-06 - 100 0.111 6.169 8.67E-05 2.56 2.93E-07 3.05 200 5.408e-02 7.339 8.84E-06 3.19 1.70E-08 3.99 400 2.781e-02 6.956 1.11E-06 3.12 1.29E-19 3.88 36 Table 3.3: The cell average of the error e for equation (3.34) when using P 1 and P 2 polynomials. L2 -norm of the cell average of e P 1 polynomial P 2 polynomial initial discretization N hmax Λ L2 error order L2 error order R uht = P− ut 50 0.202 6.056 3.84E-04 - 5.18E-07 - (u − u)dx = 0 Ω h 100 0.111 6.169 8.23E-05 2.54 4.70E-08 3.96 200 5.408e-02 7.339 8.42E-06 3.19 1.08E-09 5.29 400 2.781e-02 6.956 1.06E-06 3.11 3.33E-11 5.23 L2 projection 50 0.202 6.056 3.87E-04 - 5.56E-06 - 100 0.111 6.169 1.19E-04 1.94 1.18E-06 2.56 200 5.408e-02 7.339 1.38E-05 3.02 1.11E-07 3.31 400 2.781e-02 6.956 2.02E-06 2.89 7.65E-09 4.02 P− projection 50 0.202 6.056 3.35E-04 - 1.07E-06 - 100 0.111 6.169 7.38E-05 2.50 9.30E-08 4.03 200 5.408e-02 7.339 7.72E-06 3.16 3.47E-09 4.60 400 2.781e-02 6.956 9.81E-07 3.10 1.72E-10 4.52 37 extension from the downwind cell to the right. We use P 2 polynomials on a uniform mesh with N = 100. Define piecewise constants s(x) such that in cell Ij vej − vj s(x) = sj = − 1, vej − uj where uj denotes the cell average of the exact solution u in cell Ij . We compute and observe that ksk∞,Ω = 7.90 × 10−4 , indicating that the computable quantity vej − vj is a good estimate of the local error vej − uj . From Figure 3.3, we observe that the computable quantity |¯ v − v˜| captures the profile of the local error |u − v| (computed at the middle point in each cell) well. 10-4 -5 10 Error -6 10 -7 10 -8 10 1 2 3 4 5 6 X Figure 3.3: Comparison between |u − v| (solid line) and |¯ v − v˜| (dashed line). 38 Example 3.2. We solve the following initial boundary value problem    u + ux = 0   t u(x, 0) = sin(x) . (3.35)     u(0, t) = sin(−t) The exact solution to this problem is u(x, t) = sin(x − t). We use third order SSP Runge-Kutta discretization in time, take ∆t = 0.1h2min to reduce the time error, and test the example with both P 1 and P 2 polynomials. The same quantities as in Example 3.1 on the same kind of random meshes of N cells are computed. The initial solution is obtained as given in Section 3.2.1. In Table 3.4 we observe that the error between the DG solution and the exact solution is (2k + 1)-th order superconvergent at the downwind point and (k + 2)-th order superconvergent at the other Radau points. Table 3.4: The error e at the Radau points for (3.35) when using P 1 and P 2 poly- nomials. 1st Radau point 2nd Radau point downwind point Polynomial N hmax Λ error order error order error order P1 50 0.202 6.056 6.06E-05 - 2.88E-05 - 100 0.111 6.169 8.92E-06 3.16 4.30E-06 3.14 200 5.408e-02 7.339 1.02E-06 3.04 4.73E-07 3.09 400 2.781e-02 6.956 1.25E-07 3.15 5.82E-08 3.15 P2 50 0.202 6.056 4.10E-07 - 3.11E-07 - 1.30E-08 - 100 0.111 6.169 3.20E-08 4.21 2.38E-08 4.24 5.51E-10 5.22 200 5.408e-02 7.339 1.84E-09 4.00 1.38E-09 3.99 1.45E-11 5.06 400 2.781e-02 6.956 9.36E-11 4.48 1.05E-10 3.87 4.31E-13 5.28 Table 3.5 shows the (k + 2)-th order superconvergence of the error ξ in the L2 - 39 norm, demonstrating that the estimate in (3.5) is sharp. Table 3.5: The error ξ for (3.35) when using P 1 and P 2 polynomials. L2 -norm of ξ P 1 polynomial P 2 polynomial N hmax Λ L2 error order L2 error order 50 0.202 6.056 1.24E-04 - 6.63E-07 - 100 0.111 6.169 1.87E-05 3.13 5.32E-08 4.17 200 5.408e-02 7.339 2.10E-06 3.06 3.03E-09 4.01 400 2.781e-02 6.956 2.58E-07 3.15 2.30E-10 3.88 As in Example 3.1, we also test the superconvergence for the cell average. Table 3.6 shows the result for Example 3.2 using the method in Section 3.2.1 as well as the L2 - and P− -projections for the initial discretization. We observe results similar to those in the periodic case. Example 3.3. We solve the following two-dimensional problem  u + u + u = 0 t x y , (3.36)  u(x, y, 0) = sin(x + y) with periodic boundary condition on the domain [0, 2π]2 . The exact solution is u(x, y, t) = sin(x + y − 2t). We use a random rectangular mesh defined as 0 = x 1 < · · · < xNx + 1 = 2π, 0 = y 1 < · · · < yNy + 1 = 2π 2 2 2 2 and Ii,j = [xi− 1 , xi+ 1 ] × [yj− 1 , yj+ 1 ]. 2 2 2 2 40 Table 3.6: The cell average of the error e for (3.35) when using P 1 and P 2 polyno- mials. L2 norm of the cell average of e P 1 polynomial P 2 polynomial initial discretization N hmax Λ L2 error order L2 error order uht = P− ut 50 0.202 6.056 1.11E-04 - 3.75E-08 - − − uh j+ 1 = P− uj+ 1 100 0.111 6.169 1.67E-05 3.12 1.63E-09 5.18 2 2 200 5.408e-02 7.339 1.90E-06 3.04 4.07E-11 5.16 400 2.781e-02 6.956 2.35E-07 3.14 1.14E-12 5.37 L2 projection 50 0.202 6.056 1.86E-04 - 2.59E-06 - 100 0.111 6.169 5.48E-05 2.02 3.59E-07 3.26 200 5.408e-02 7.339 6.86E-06 2.91 3.11E-08 3.43 400 2.781e-02 6.956 1.02E-06 2.87 2.86E-09 3.59 P− projection 50 0.202 6.056 1.01E-04 - 2.49E-07 - 100 0.111 6.169 1.53E-05 3.12 1.48E-08 4.66 200 5.408e-02 7.339 1.72E-06 3.06 5.73E-10 4.55 400 2.781e-02 6.956 2.06E-07 3.20 2.48E-11 4.72 41 We define the approximation space as © ª Vhk = uh : uh |Ii,j ∈ Qk (Ii,j ), 1 ≤ i ≤ Nx , 1 ≤ j ≤ Ny , where Qk (Ii,j ) denotes all the tensor product polynomials of degree at most k in x and in y on Ii,j . The Gauss-Radau projection P− is defined as follows: Z (P− u − u)vh dxdy = 0 Ii,j for any vh ∈ Vhk−1 , Z yj+ 1 2 (P− u(xi+ 1 , y) − u(xi+ 1 , y))wh (y)dy = 0 2 2 yj− 1 2 for any wh ∈ P k−1 , Z xi+ 1 2 (P− u(x, yj+ 1 ) − u(x, yj+ 1 ))zh (x)dx = 0 2 2 xi− 1 2 for any zh ∈ P k−1 , and P− u(xi+ 1 , yj+ 1 ) = u(xi+ 1 , yj+ 1 ). 2 2 2 2 We also use an upwind flux and ninth order SSP Runge-Kutta discretization in time with ∆t = 0.1hmin . We test the example with both Q1 and Q2 polynomials, and compute ξ = uh − P− u, the error of the cell average, as well as the error on the downwind-biased Radau points and the downwind point. For simplicity, we do not consider the Radau points one by one, but select the one which gives the largest error in each cell. The initial solution is given by the P− projection. It appears that similar superconvergence results are valid also in two-dimensions. However, the technique of the proof in this chapter, in particular the part related to 42 the special projections, does not seem to be easily extendable to two-dimensions. Table 3.7: Superconvergence results for equation (3.36) when using Q1 and Q2 poly- nomials. Q1 polynomial Q2 polynomial Error Nx × N y hmax Λ L2 error order L2 error order kP− u − uh k 10 × 10 0.997 3.722 2.65E-02 - 8.40E-04 - 20 × 20 0.552 4.809 4.41E-03 3.04 6.97E-05 4.22 40 × 40 0.270 7.339 5.12E-04 3.02 3.64E-06 4.13 80 × 80 0.133 6.326 6.33E-05 2.96 2.33E-07 3.89 ku − uh k 10 × 10 0.997 3.722 4.19E-02 - 6.47E-04 - 20 × 20 0.552 4.809 7.65E-03 2.88 5.26E-05 4.25 40 × 40 0.270 7.339 9.12E-04 2.98 2.20E-06 4.44 80 × 80 0.133 6.326 1.14E-04 2.94 1.16E-07 4.17 maxi,j |(u − uh )(x, y)| 10 × 10 0.997 3.722 2.28E-02 - 1.56E-03 - (x,y) is the Downwind- 20 × 20 0.552 4.809 5.65E-03 2.37 1.35E-04 4.14 biased Radau points in Ii,j 40 × 40 0.270 7.339 6.39E-04 3.05 8.88E-06 3.81 80 × 80 0.133 6.326 1.01E-04 2.62 7.15E-07 3.57 maxi,j |(u − uh )(x, y)| 10 × 10 0.997 3.722 1.99E-02 - 3.72E-04 - (x,y) is the downwind 20 × 20 0.552 4.809 2.79E-03 3.33 4.09E-05 3.74 point in Ii,j 40 × 40 0.270 7.339 3.93E-04 2.74 1.49E-07 4.64 80 × 80 0.133 6.326 4.65E-05 3.02 1.34E-07 3.42 3.4 Concluding remarks We have studied the behavior of the error between the upwinded DG solution and the exact solution for sufficiently smooth solutions of linear conservation laws. We prove that under suitable initial discretization, the error between the DG solution and the exact solution is (k + 2)-th order superconvergent at the downwind-biased Radau points. Moreover, numerical experiments show that the convergent rate is (2k + 1)-th order superconvergent at the downwind point as well as in the L2 -norm 43 of the cell average. We also prove that the DG solution is superconvergent with the rate k + 2 towards a particular projection of the exact solution. Chapter 4 Negative-order norm error estimates for linear hyperbolic equations involving δ-functions In this chapter, we develop and analyze DG methods for solving hyperbolic conser- vation laws involving δ-functions. δ-function has many equivalent definitions, and one of them is the following weak form: Z δ(x)v(x)dx = v(0), v(x) ∈ C(R). R As mentioned in Chapter 1, the DG methods are based on a weak form, and can be designed to approximate δ-functions directly. However, in DG methods, the test functions are completely discontinuous across the cell interfaces. If the δ-function is placed at the cell interface, then the numerical scheme is not well-defined. In [68], the author introduced the distribution theory for discontinuous test function, which extended the definition of δ-function as Z v(0+) + v(0−) δ(x)v(x)dx = , v(x) is piecewise continuous, R 2 44 45 which further completes the definition of the DG schemes. We consider negative-order norm error estimates. Such norms can be used to detect the oscillations of a function. In [34], Cockburn et al. proved high order superconvergence error estimates of DG methods including their divided differences for hyperbolic equations with smooth solutions in negative-order norms. They also demonstrated that the application of the post-processing techniques of Bramble and Schatz [17] can yield superconvergence in the strong L2 -norm. Other related works include [92, 103, 111, 91], where one-sided filter and local derivative post-processing were considered. We will extend the work in [34], and consider the general case where the exact solutions are not smooth. The first example of non-smooth solutions for hyperbolic equations is the follow- ing problem ut + ux = 0, (x, t) ∈ R × (0, T ], (4.1) u(x, 0) = u0 (x), x ∈ R, where the initial solution u0 (x) has compact support, with a discontinuity at x = 0, but is otherwise smooth. Clearly, the exact solution of (4.1) is discontinuous along the characteristic line x = t and the numerical DG solution has spurious oscillations around this discontinuity line, which we refer to as the pollution region. There are not too many works in the literature studying error estimates of DG methods for problems with discontinuous solutions. The first work in this direction seems to be that of Johnson et al. [63, 64, 65] for DG methods in both space and time. They have shown that, with linear space-time elements, the width of the pollution region is of the size at most O(h1/2 log 1/h). More recently, Cockburn and Guzm´an [31] and Zhang and Shu [117] revisited this problem with the RKDG methods and obtained similar results. Especially, in [31], the left boundary of the pollution region is shown to be at most O(h2/3 log(1/h)) from the singularity for piecewise linear DG method with second order Runge-Kutta time discretization on uniform meshes. The first problem we consider in this chapter is (4.1) with the initial condition u0 (x) having 46 a δ-singularity at x = 0. We consider a semi-discrete DG method and use the result in [117] to prove superconvergence results estimated in negative-order norms outside the pollution region. By convolving the DG solution with a suitable kernel, the post- processed approximation is (2k+1)-th order accurate in a region slightly smaller than the one above. The rate of convergence agrees with that in [34], in which the initial datum u0 (x) was assumed to be sufficiently smooth. Hyperbolic conservation laws with source terms have been analyzed by several authors [21, 47, 66, 99, 72]. In particular, in [47], the authors studied the following problem ut + f (u)x = g(x, t), (x, t) ∈ R × (0, T ], (4.2) u(x, 0) = u0 (x), x ∈ R, where f is a smooth convex function (f 00 (u) > 0 for all u) and g(x, t) = Gx (x, t) with G being a bounded, piecewise smooth function, and constructed L∞ -stable Godunov-type difference schemes. In [98], Santos and de Oliveira studied hyperbolic conservation laws whose source terms contain δ-singularities, and investigated the convergence of numerical discretization by using a finite volume scheme. Later, they considered a class of high resolution methods in [79]. In [78], Noussair studied the wave behavior of (4.2), where the source term also depends on u but not on the time variable t. We note that all these previous works did not provide any error estimates in the smooth region away from the singularities. In this chapter, we investigate a simpler case by assuming f (u) = u and g(x, t) = G0 (x) in the source term in (4.2), where G(x) is a step function which does not depend on the time variable t. We show that by convolving the DG solution with a suitable kernel, the post-processed approximation turns out to be (2k + 1)-th order superconvergent in the smooth region. 47 4.1 Post-processing 4.1.1 Convolution kernel Now, we proceed to describe the type of post-processing to be considered, following Bramble and Schatz [17]. Let χ be the indicator function of the interval (− 12 , 12 ), We define recursively the functions ψ (l) as ψ (1) = χ, ψ (n+1) = ψ (n) ∗ ψ (1) , for n ≥ 1. The numerical solution is post-processed by convolving it with a kernel K (ν,l) (x) which satisfies the following properties: (1) It has a compact support. (2) It reproduces polynomials p of degree ν − 1 by convolution: K (ν,l) ∗ p = p. (3) It is a linear combination of B-splines and is of the form X K (ν,l) (x) = kγν,l ψ (l) (x − γ). γ∈Z The weights kγν,l ∈ R are chosen so that (2) is satisfied. See [17, 34] for more details. (ν,l) (l) We also define KH (x) = K (ν,l) (x/H)/H and ψH (x) = ψ (l) (x/H)/H and it is not difficult to verify that (β−α) Dα (ψ (β) ) ∗ v = ψH α ∗ ∂H v, 1 where ∂H v(x) = H (v(x+ 12 H)−v(x− 12 H)). In general, we take H = nh, n = 1, 2, · · · . This property is important as it allows us to express derivatives of the convolution with the kernel in terms of simple difference quotients. 48 4.1.2 An approximation result Let us investigate the relationship between u − Kh2k+2,k+1 ∗ uh and the negative-order norm estimates of divided differences of the error u − uh . Theorem 4.1.1 (Bramble and Schatz [17]). Suppose the kernel Khν,l satisfies the properties listed in Section 4.1.1. Let v be a function in L2 (Ω1 ), where Ω1 is an open set in Ω, and u be a function in H ν (Ω1 ). Further assume Ω0 to be an open set in Ω1 such that Ω0 + 2supp(Khν,l ) ⊂⊂ Ω1 . Then we have hν X ku − Khν,l ∗ vkΩ0 ≤ C1 |u|ν,Ω1 + C1 C2 k∂hα (u − v)k−l,Ω1 , ν! 0≤α≤l P where C1 = γ∈Z |kγν,l | and C2 only depends on Ω0 , Ω1 , ν, and l. There is a straightforward corollary. Corollary 4.1.1. Suppose the conditions in Theorem 4.1.1 are satisfied. Further assume that k∂hα (u − v)k−l,Ω1 ≤ Chµ is valid for all α ≤ l and ν ≥ µ. Then we have ku − Khν,l ∗ vkΩ0 ≤ Chµ , where C only depends on Ω0 , Ω1 , ν, and l. 4.2 Singular initial condition Let us consider problem (4.1) and use upwind fluxes. We first state the main results in Theorem 4.2.1 and then give the proofs. We provide the negative-order norm error estimates in the whole space as well as in the region away from the singularities. 49 4.2.1 Main results The following lemma is the semi-discrete version of the result in Zhang and Shu [117]. For completeness we will give its proof in Section 4.5. Lemma 4.2.1. Let u be the exact solution of the initial value problem (4.1), where the initial condition u0 (x) ∈ C k+2 except for one singularity at x = 0. Let uh be the solution of the DG method (2.5) at time T , where the finite element space Vh is made up of the piecewise polynomials of degree k ≥ 1. Suppose h is the maximum cell length. Then there holds the following error estimate ku(T ) − uh (T )kΩ\RT ≤ Chk+1 , (4.3) where RT = (T − Ch1/2 log(1/h), T + Ch1/2 log(1/h)), and the bounding constant C > 0 does not depend on h. We will use Lemma 4.2.1 to prove the following theorem. Theorem 4.2.1. Suppose u ∈ C 2k+2 and the conditions of Lemma 4.2.1 are satisfied. Then by taking Ω0 + 2supp(Kh2k+2,k+1 ) ⊂⊂ Ω1 ⊂⊂ Ω\RT , we have ku(T ) − uh (T )k−(k+1) ≤ Chk , (4.4) ku(T ) − uh (T )k−(k+2) ≤ Chk+1/2 , (4.5) ku(T ) − uh (T )k−(k+1),Ω1 ≤ Ch2k+1 , (4.6) ku(T ) − Kh2k+2,k+1 ∗ uh (T )kΩ0 ≤ Ch2k+1 , (4.7) where the positive constant C does not depend on h. Here the mesh is assumed to be uniform for (4.7) but can be regular and non-uniform for the other three inequalities. Remark 4.2.1. To obtain (4.7), we have to assume the mesh is uniformly dis- tributed, that is hj = h, ∀ j. This is a result of the negative-order norm estimates 50 of the divided differences. Actually, we denote w = ∂h u and wh = ∂h uh . Clearly, w satisfies (4.1) with initial condition w(x, 0) = ∂h u(x, 0). If we shift the mesh by h 2 , then wh satisfies numerical scheme (2.5). By the same analysis for the proof of (4.6), we obtain k∂h (u − uh )k−(k+1),Ω1 = kw − wh k−(k+1),Ω1 ≤ Ch2k+1 . The estimates for higher order divided differences can be obtained by exactly the same line in this remark. Therefore, (4.7) follows directly from Corollary 4.1.1. Remark 4.2.2. The error estimates in the −(k + 1)-th order norm are used for problems with singular initial conditions while the estimates in the −(k + 2)-th order norm are used for problems with singular source terms. 4.2.2 A proof of Theorem 4.2.1 In this subsection, we give the discretization and prove the first three estimates in Theorem 4.2.1. Initial discretization From now on, we assume the δ-singularity of the initial datum is contained in cell Ii . For simplicity, we also assume the singularity is concentrated at 0, denoted as δ(x). We apply the L2 -projection Pk to discretize the initial condition to obtain kuh (0)k ≤ Ch−1/2 . At t = 0, for any function φ ∈ C0∞ (Ω), we have, for the cell Ii which contains the δ-singularity, (u − uh , φ)i = (u − uh , φ − Pk φ)i = (u, φ − Pk φ)i ≤ kφ − Pk φk∞,Ii 1 ≤ Chk+ 2 |φ|k+1,Ii . 51 In other cells, following the same analysis as above, we have (u − uh , φ)j = (u − Pk u, φ − Pk φ)j ≤ Ch2k+2 |u0 |k+1,Ij |φ|k+1,Ij . (4.8) The −(k + 1)-th order error estimate on Ω In this subsection, we proceed to prove (4.4). The proof mostly follows [34]. We begin by considering the solution to the dual problem: Find a function φ such that φ(·, t) satisfies φt + φx = 0, (x, t) ∈ Ω × (0, T ), (4.9) φ(x, T ) = Φ(x), x ∈ Ω. Assuming Φ is an arbitrary function in C0∞ (Ω), we have, following [34], Z T (u(T ) − uh (T ), Φ) = (u − Pk u, φ)(0) − [((uh )t , φ) + (uh , φt )]dt (4.10) 0 Z T N X = (u − Pk u, φ)(0) − [uh ](φ − Pk φ)+ |j− 1 (4.11) 2 0 j=1 Z à N !1/2 T X ≤ Chk+1/2 |Φ|k+1 + Chk+1/2 |Φ|k+1 [uh ]2j− 1 dt. 2 0 j=1 Using the Cauchy-Schwartz inequality and Lemma 2.4.1, we have Z à N !1/2 ÃZ N !1/2 T X T X [uh ]2j− 1 dt ≤ T 1/2 [uh ]2j− 1 dt 2 2 0 j=1 0 j=1 1/2 ≤T kuh (0)k ≤ CT 1/2 h−1/2 . 52 Combining the above, we can see (u(T ) − uh (T ), Φ) ku(T ) − uh (T )k−(k+1) = sup Φ∈C0∞ (Ω) kΦkk+1 Chk+1/2 |Φ|k+1 + CT 1/2 hk |Φ|k+1 ≤ sup Φ∈C0∞ (Ω) kΦkk+1 ≤ CT 1/2 hk + Chk+1/2 . Now, we consider the extension to higher dimensions. The proof of the following corollary is straightforward and is similar to the one-dimensional case, and is thus omitted. Corollary 4.2.1. Let Ω be an open set in Rd , and u be the exact solution of the following initial value problem Pd ut + j=1 uxj = 0, (x, t) ∈ Ω × (0, T ], u(x, 0) = δ(f (x)), x ∈ Ω, where f (x) : Rd → R is a smooth function. Denote Γh = {K} as a regular triangu- lation of Rd , whose elements K are open and have diameter hK less than or equal to h. In each K, denote ∂K− and ∂K+ as the inflow and outflow edges respectively. Let uh be the DG approximation which satisfies d X d X d X (uht , vh )K = (uh , (vh )xi )K + (u− + h , vh )∂K− − (u− − h , vh )∂K+ , vh ∈ Vh , i=1 i=1 i=1 where the finite element space Vh is made up of the piecewise polynomials of degree k ≥ 1. Suppose the total measure of the cells which contain δ-singularities initially is mh, then there holds the following estimate √ ku(T ) − uh (T )k−k−1 ≤ C mT 1/2 hk + Chk+d/2 , (4.12) 53 where the bounded constant C > 0 does not depend on h or T . The −(k + 2)-th order error estimate on Ω We shall prove (4.5). To do so, we apply P+ to estimate the term (uht , φ) + (uh , φt ). By using (2.4) and Lemma 2.4.3, we obtain (uht , φ) + (uh , φt ) = (uht , P⊥ + φ) + (uht , P+ φ) − (uh , φx ) = (uht , P⊥ + φ) + H(uh , P+ φ) − H(uh , φ) = (uht , P⊥ + φ). (4.13) Integrating in t, we obtain Z T Z T (uht , φ) + (uh , φt )dt = (uh , P⊥ + φ)(T ) − (uh , P⊥ + φ)(0) − (uh , P⊥ + φt )dt. (4.14) 0 0 Applying Lemma 2.4.1, we have Z T Z T ¡ ¢ (uht , φ) + (uh , φt )dt ≤ kuh (0)k k(P⊥ + φ)(0)k + k(P⊥ + φ)(T )k + kuh (0)k kP⊥ + φt (t)kdt 0 0 µ Z T ¶ k+1 k+1 ≤ Ckuh (0)k h |Φ|k+1 + h |Φ|k+2 dt 0 ≤ C(1 + T )hk+1 kuh (0)kkΦkk+2 . From the above we observe (u(T ) − uh (T ), Φ) ku(T ) − uh (T )k−(k+2) = sup Φ∈C0∞ (Ω) kΦkk+2 1 1 Chk+ 2 |Φ|k+1 + C(1 + T )hk+ 2 kΦkk+2 ≤ sup Φ∈C0∞ (Ω) kΦkk+2 1 ≤ C(1 + T )hk+ 2 . 54 4.2.3 The negative-order error estimate on Ω1 ⊂⊂ Ω\RT We proceed to prove (4.6). To estimate the negative-order norm of u − uh at time T on Ω1 , we need to assume Φ ∈ C0∞ (Ω1 ) instead of C0∞ (Ω). Moreover, we also assume the exact solution u ∈ C(Ω), this is because we are allowed to modify the exact solution in the cell which contains the δ-singularity, keeping the numerical solution uh untouched. More details of this assumption can be found in [117, 31] or Section 4.5.2. Therefore, in (4.11), we have N X N X + [uh ](φ − Pk φ) |j− 1 = [uh − Pk u + Pk u − u](φ − Pk φ)+ |j− 1 2 2 j=1 j=1 ≤ Chk (ku − Pk ukΩ1 + kPk u − uh kΩ1 ) |φ|k+1 ≤ Chk (ku − Pk ukΩ1 + ku − uh kΩ1 ) |φ|k+1 ≤ Ch2k+1 |φ|k+1 , where we use Lemmas 2.3.1 and 2.3.2 in the second inequality and Lemma 4.2.1 in the last one. Inserting this into (4.11), we have (u(T ) − uh (T ), Φ) ≤ (u − Pk u, φ)(0) + Ch2k+1 |φ|k+1 . (4.15) By using (4.8), we obtain the estimate we want (u(T ) − uh (T ), Φ) ku(T ) − uh (T )k−(k+1),Ω1 = sup Φ∈C0∞ (Ω1 ) kΦkk+1 Ch2k+1 kΦkk+1 + Ch2k+2 kΦkk+1 ≤ sup Φ∈C0∞ (Ω1 ) kΦkk+1 ≤ Ch2k+1 , where the constant C > 0 is independent of h. 55 4.3 Singular source term In this section, we briefly discuss a linear inhomogeneous evolution equation of a function u(x, t) : Ω × (0, ∞) → R of the form   u (x, t) + Lu(x, t) = f (x, t), (x, t) ∈ Ω × (0, ∞), t (4.16)  u(x, 0) = 0, x ∈ Ω, with L being a linear differential operator that does not involve time derivatives. If we multiply the above equation by a smooth function φ(x, t), then integrate over space and time, we obtain Z ∞ Z Z ∞ Z [φut + φLu] dxdt = f (x, t)φ(x, t)dxdt. 0 Ω 0 Ω Integrating by parts and assuming zero boundary condition, we have Z ∞ Z Z ∞ Z ∗ [uφt + uL φ] dxdt + f (x, t)φ(x, t)dxdt = 0, (4.17) 0 Ω 0 Ω where L∗ is the dual operator of L. Definition 4.3.1. The function u(x, t) is called a weak solution of the equation (4.16), if (4.17) holds for all functions φ ∈ C01 (Ω × R+ ). 4.3.1 Duhamel’s principle Now, we consider linear hyperbolic conservation laws with source terms. To deal with such problems we apply Duhamel’s principle, which is applicable to linear parabolic and hyperbolic PDE and yields an integral representation in terms of the solutions of more tractable PDEs. 56 Lemma 4.3.1 (Duhamel’s principle). The solution to (4.16) is Z t u(x, t) = (P s f )(x, t)ds, 0 where P s f is the solution of the problem   P (x, t) + LP (x, t) = 0, (x, t) ∈ Ω × (s, ∞), t (4.18)  P (x, s) = f (x, s), x ∈ Ω. Notice that P s f is the solution to the homogeneous PDE with the source term f serving as the initial condition at time t = s. To prove the lemma, we can simply check that the expression of u satisfies (4.16). More details and a proof can be found in [62], in which the PDE is a second order wave equation. The above lemma requires suitable regularity of u. However, Duhamel’s principle is also valid in the following weak sense. Lemma 4.3.2. Suppose u(x, t) is the weak solution of equation (4.16), then Z t u(x, t) = (P s f )(x, t)ds 0 in the sense of distribution, where P s f is the weak solution of equation (4.18). The proof directly follows from the definition of the weak solution and the proof of Duhamel’s principle, so we omit it here. Finally, we extend Duhamel’s principle to the DG schemes. For simplicity, we consider the following equation   u (x, t) + u (x, t) = δ(x), (x, t) ∈ Ω × (0, T ), t x (4.19)  u(x, 0) = u (x), x ∈ Ω, 0 with u0 = 0. For general smooth u0 (x), the same result can be obtained by superpo- sition. We define the finite element approximation uh : [0, T ] → Vh as the solution 57 to (uht , χ)j = Hj (uh , χ) + (δ(x), χ)j , ∀χ ∈ Vh , (4.20) uh (0) = 0, where Hj (·, ·) is the DG bilinear form defined in (2.4). Then the semi-discrete version of Duhamel’s principle is given in Lemma 4.3.3. Rt Lemma 4.3.3. The solution of (4.20) can be written in the form uh = 0 ps (x, t)ds where ps (x, t) is the solution of the following scheme: find p ∈ Vh such that (pt , χ)j = Hj (p, χ), ∀χ ∈ Vh , (4.21) p(s) = Pk δ(x). Rt The proof is straightforward, since uh in (4.20) and 0 ps (x, t)ds share the same initial condition and the same system of ODEs, noticing the fact that (Pk δ, χ) = (δ, χ). In what follows, we would like to rewrite the inhomogeneous equations (4.19) and (4.20) into homogeneous ones (4.26) and (4.21) by using Lemma 4.3.2 and Lemma 4.3.3 respectively. Then we apply the estimates of P s f − ps , which have been given in Theorem 4.2.1, to prove the main result, Theorem 4.3.1. 4.3.2 Error estimates We first state the main result Theorem 4.3.1 and then give the proof. Theorem 4.3.1. Suppose u is the exact solution of equation (4.19), and uh is the numerical solution which satisfies (4.20). Denote RT = Ii ∪ (T − C log(1/h)h1/2 , T + C log(1/h)h1/2 ), where Ii is the cell which contains the concentration of the δ-singularity 58 on the source term. Then we have the following estimates ku(T ) − uh (T )k−(k+1) ≤ Chk , (4.22) ku(T ) − uh (T )k−(k+2) ≤ Chk+1/2 , (4.23) ku − uh k−(k+1),Ω1 ≤ Ch2k+1 , (4.24) ku(T ) − Kh2k+2,k+1 ∗ uh (T )kΩ0 ≤ Ch2k+1 , (4.25) where Ω0 + 2supp(Kh2k+2,k+1 ) ⊂⊂ Ω1 ⊂⊂ R\RT . Here the mesh is assumed to be uniform for (4.25) but can be regular and non-uniform for the other three inequalities. Remark 4.3.1. As mentioned in Remark 4.2.1, (4.25), which requires uniform meshes, follows from (4.24). Moreover, we also skip the proofs of (4.22) and (4.23), since they follow easily from (4.4) and (4.5) in Theorem 4.2.1. Now we proceed to prove (4.24). Denote v s as the exact solution of the following equation ut + ux = 0, (x, t) ∈ Ω × (s, T ], (4.26) u(x, s) = δ(x), x ∈ Ω, and vhs as the solution of the numerical scheme (4.21). For convenience, if s = 0, the superscript will be omitted. We consider the dual problem defined the same way as (4.9). By Lemma 4.3.3 and Lemma 4.3.2, we have Z T (u − uh , Φ)(T ) = (v s − vhs , Φ)(T )ds. (4.27) 0 By using (4.10) and (4.13), and the fact that vh is the L2 -projection of v at t = 0, we obtain Z T s (v − vhs , Φ)(T ) = ((v − vh )(0), P⊥ + φ(s)) − (vht (t − s), P⊥ + φ(t))dt, s 59 which further yields (u − uh , Φ)(T ) = Π1 − Π2 , RT RT RT where Π1 = 0 ((v − vh )(0), P⊥ + φ(s))ds, and Π2 = 0 s (vht (t− s), P⊥ + φ(t))dtds. First consider the second term, Z T Z t Π2 = − ((vh )s (t − s), P⊥ + φ(t))dsdt 0 0 Z T = (vh (t) − vh (0), P⊥+ φ(t))dt. 0 Therefore, Z T Z T (u − uh , Φ)(T ) = (v(0), P⊥ + φ(s))ds − (vh (t), P⊥ + φ(t))dt 0 0 = G1 − G2 . R T +τ We claim G1 = 0. Actually, for any τ ∈ Ii , τ Φ(x)dx does not depend on τ , since Φ(x) vanishes in the neighborhood of x = 0 and x = T . Therefore, µ Z T ¶ µ Z x+T ¶ G1 = δ(x), P⊥ + φ(x, s)ds = δ(x), P⊥ + Φ(y)dy = 0. 0 x i Now, we only need to estimate G2 . Since (vh (t), P⊥ ⊥ + φ(t)) = (vh (t) − v(t) + v(t) − Pk−1 v, P+ φ(t)), by Lemma 4.2.1 and Lemma 2.3.2, we have G2 ≤ Ch2k+1 |φ|k+1 . Finally, we obtain ku − uh k−(k+1),Ω1 ≤ Ch2k+1 and complete the proof of Theorem 4.3.1. 60 4.4 Numerical examples In this section, we provide numerical experiments to demonstrate our theoretical results for the post-processor and to illustrate the performance of the DG schemes. We denote by d the distance between the singularities and the region under consid- eration. In all the figures, if not otherwise stated, the numerical solutions are plotted using six Gaussian points in each cell. 4.4.1 Singular initial condition Example 4.1. We solve the following problem ut + ux = 0, (x, t) ∈ [0, π] × (0, 1], (4.28) u(x, 0) = sin(2x) + δ(x − 0.5), x ∈ [0, π], with periodic boundary condition u(0, t) = u(π, t). Clearly, the exact solution is u(x, t) = sin(2x − 2t) + δ(x − t − 0.5). We use a ninth order SSP Runge-Kutta discretization in time [46] and take the time step ∆t = 0.1h. We test the example by using P k polynomials with k = 1, 2, 3 on uniform meshes, and compute the L2 -norm of the error after post-precessing in the region away from the singularity at t = 0.5. By taking d = 0.2, the region under consideration is [0, 0.8] ∪ [1.2, π]. In Table 4.1, we can observe at least (2k + 1)-th order convergence. Moreover, we observe that the rate of convergence settles to the dN 0.2×500 asymptotic value when the total number of cells is around π = π ≈ 30, no matter which degree of polynomials we use. The initial discretization is obtained by taking the L2 projection. Figure 4.1 shows the numerical solution with and without post-processing. We 61 Table 4.1: L2 -norm of the error between the numerical solution and the exact solution for (4.28) after post-processing in the region away from the singularity. P 1 polynomial P 2 polynomial P 3 polynomial N d error order error order error order 200 0.2 6.88E-05 - 8.40e-07 - 1.48E-09 - 300 0.2 1.41E-05 3.92 3.56e-10 19.2 3.98E-13 20.3 400 0.2 5.89E-06 3.02 1.98e-11 10.1 4.42E-16 23.7 500 0.2 3.01E-06 3.01 6.13e-12 5.25 7.49E-17 7.95 600 0.2 1.74E-06 3.00 2.37e-12 5.21 1.76E-17 7.94 use P 2 polynomials and take h = 0.01. From the figure we observe some localized oscillations near the discontinuity and that the post-processor does not significantly smear the singularity. 50 50 40 40 30 30 U U 20 20 10 10 0 0 -10 -10 0 0.5 1 1.5 2 2.5 3 0 1 2 3 x x Figure 4.1: Numerical solution for (4.28) at t = 0.5 with (right) and without (left) post-processing. Example 4.2. We consider the following two dimensional problem ut + ux + uy = 0, (x, y, t) ∈ [0, 2π] × [0, 2π] × (0, 1], (4.29) u(x, 0) = sin(x + y) + δ(x + y − 2π), (x, y) ∈ [0, 2π] × [0, 2π], 62 with periodic boundary condition. Clearly, the exact solution is u(x, t) = sin(x + y − 2t) + δ(x + y − 2t) + δ(x + y − 2t − 2π). We use Qk polynomial approximation spaces with k = 1 and 2, where Qk is the space of tensor product polynomials of degree at most k ≥ 0. We also apply the same time discretization as in Example 4.1 and compute the L2 -norm of the error after post-precessing in the region away from the singularity at t = 0.5. Moreover, we take d = 0.4. In Table 4.2, we can observe (2k + 1)-th order convergence. Table 4.2: L2 -norm of the error between the numerical solution and the exact solution for (4.29) after post-processing in the region away from the singularity. Q1 polynomial Q2 polynomial N d error order error order 400 0.4 2.60E-05 - 3.23e-08 - 500 0.4 1.24E-05 3.32 2.47e-10 20.0 600 0.4 7.16E-06 3.01 1.19e-11 16.6 700 0.4 4.50E-06 3.01 5.11e-12 5.47 800 0.4 3.01E-06 3.02 2.53e-12 5.29 It appears that similar results are valid in two dimensions. However, the technique of proof in this chapter, in particular the part related to the special projections in (2.10) and (2.11), does not seem to be easily extendable to two dimensions. Moreover, Figure 4.2 shows the numerical solution by plotting the numerical cell averages. We use Q2 polynomials and take N = 100. From the figure we can observe two lines of δ-singularities. Even though the theory in this chapter is given only for scalar linear equations for simplicity, it generalizes to linear systems in a straightforward way. 63 Z Y X 3 2.5 2 1.5 U 1 0.5 0 -0.5 -1 0 1 2 3 4 5 6 x,y Figure 4.2: Numerical solution (left) and the cut plot along x = y (right) for (4.29) at t = 0.5. Example 4.3. We solve the following linear system ut − vx = 0, (x, t) ∈ [0, 2] × (0, 0.4], vt − ux = 0, (x, t) ∈ [0, 2] × (0, 0.4], (4.30) u(x, 0) = δ(x − 1), v(x, 0) = 0, x ∈ [0, 2]. Clearly, the exact solution (the Green’s function) is 1 1 1 1 u(x, t) = δ(x − 1 − t) + δ(x − 1 + t), v(x, t) = δ(x − 1 + t) − δ(x − 1 − t). 2 2 2 2 We use a third order SSP Runge-Kutta discretization in time [46] and take the time step ∆t = 0.1h. Figure 4.3 shows the numerical solutions at t = 0.4 with P 3 polynomials and h = 0.01. We observe that the numerical solutions capture the profiles of the exact solutions quite well. Since we have not used any limiter, there are some localized oscillations near the singularities. 64 40 40 35 30 20 25 20 U V 0 15 10 -20 5 0 -40 0 0.5 1 1.5 2 0 0.5 1 1.5 2 X X Figure 4.3: Solutions of u (left) and v (right) for (4.30) at t = 0.4. 4.4.2 Singular source term Example 4.4. We solve the following problem ut + ux = δ(x − π), (x, t) ∈ [0, 2π] × (0, 1], u(x, 0) = sin(x), x ∈ [0, 2π], (4.31) u(0, t) = 0, t ∈ (0, 1]. Clearly, the exact solution is u(x, t) = sin(x − t) + χ[π,π+t] , where χ[a,b] denotes the indicator function of the interval [a, b]. We use the same time discretization as in the previous example, and use both P 1 and P 2 polynomials to approximate the exact solution on uniform meshes. We compute the L2 -norm of the error after post-precessing in the region away from the singularities at t = 0.5. In this example, we also take d = 0.2, and the region under consideration is [0, π − 0.2] ∪ [π + 0.2, π + 0.3] ∪ [π + 0.7, 2π]. In Table 4.3, we observe (2k + 1)-th order 65 convergence. The initial discretization is again obtained by taking the L2 projection. Table 4.3: L2 -norm of the error between the numerical solution and the exact solution for (4.31) after post-processing in the region away from the singularity. P 1 polynomial P 2 polynomial N d error order error order 401 0.2 1.74E-06 - 4.29E-08 - 801 0.2 5.92E-09 8.22 6.80E-13 15.9 1601 0.2 7.36E-10 3.03 1.34E-17 12.3 3201 0.2 9.19E-11 3.01 3.86E-18 5.13 6401 0.2 1.15E-11 3.01 1.16E-19 5.07 Moreover, Figure 4.4 shows the numerical solutions with and without post- processing. We use P 2 polynomials and take h = 0.01. We observe that the post-processor does not smear the singularity and that it effectively damps out the oscillations near the left singularity. 1.5 1.5 1 1 0.5 0.5 U U 0 0 -0.5 -0.5 -1 -1 0 1 2 3 4 5 6 0 1 2 3 4 5 6 X X Figure 4.4: Numerical solutions for (4.31) at t = 0.5 with (right) and without (left) post-processing. 66 Example 4.5. We solve the following problem ut + ((x + 1)u)x = δ(x − c), (x, t) ∈ [0, 1.5] × (0, 1], u(x, 0) = 0, x ∈ [0, 1.5], (4.32) u(0, t) = 0, t ∈ (0, 1]. The exact solution is 1 u(x, t) = [H(x − c) − H(x + 1 − (c + 1)et )], 1+x where H(x) is the Heaviside function defined as   0, x < 0, H(x) =  1, x ≥ 0. π We take c = 20 and compute the solution at t = 0.5 with P 1 and P 2 polynomials. Since (x + 1) is always positive, by using upwind fluxes, we always consider u− h to be the numerical flux at the cell interfaces. For time discretization, the classical fourth order Runge-Kutta method is used with ∆t = h2 . In this example, we take d = 0.1, π π π √ and the region under consideration is [0, 20 − 0.1] ∪ [ 20 + 0.1, ( 20 + 1) e − 1.1] ∪ π √ [( 20 + 1) e − 0.9, 1.5]. In Table 4.4, we can observe (k + 1)-th and (2k + 1)-th order convergence before and after post-processing respectively. Moreover, Figure 4.5 shows the numerical solution with P 2 polynomials and h = 0.01. We use the cell averages to plot the left panel of the figure. We observe that the numerical solution agrees well with the exact solution away from the singularities. Since we have not used any limiter, there are some localized oscillations near the singularity on the right. It is interesting to observe that there are very few numerical oscillations near the left singularity. In the middle panel of Figure 4.5, we use six Gaussian points to plot, and the detailed zoom for the left singularity is given in the right panel. Clearly, the numerical solution only oscillates in the cell [0.15,0.16]. No 67 Table 4.4: L2 -norm of the error between the numerical solution with and the exact solution for (4.32) before and after post-processing in the region away from the singularity. P 1 polynomial P 2 polynomial N d error order error order before post-processing 400 0.1 7.08E-07 - 1.10E-08 - 800 0.1 1.21E-07 2.56 4.37E-11 7.98 1600 0.1 3.02E-08 2.00 5.46E-12 3.00 3200 0.1 7.55E-09 2.00 6.83E-13 3.00 6400 0.1 1.89E-09 2.00 8.53E-14 3.00 after post-processing 400 0.1 4.92E-07 - 8.65E-09 - 800 0.1 7.49E-11 12.7 4.54E-14 17.5 1600 0.1 7.43E-12 3.33 5.13E-19 16.4 3200 0.1 9.31E-13 3.00 1.76E-20 4.86 6400 0.1 1.16E-13 3.00 5.75E-22 4.94 oscillation is observed in the left figure for cell averages, and only one undershoot can be observed in the middle and right panels for which six Gaussian points are plotted. This can be explained by the size of the pollution region. In Theorem 4.3.1, we have proved that, for such singularities, RT contains only one cell. This implies that the numerical solution will oscillate within that cell, which clearly agrees with our observation. 4.5 Proof of Lemma 4.2.1 In this section we prove Lemma 4.2.1. The main line of proof is based on ideas in [31, 117]. For simplicity, we only consider a δ-singularity in (4.1), hence u0 (x) = δ(x) + f (x), where f (x) is sufficiently smooth and has compact support on the computational domain Ω. 68 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 -0.2 -0.2 -0.2 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0.15 0.16 X X X Figure 4.5: Numerical solutions for (4.32) at t = 0.5 plotted for the cell averages (left), six Gaussian points (middle) and the detailed zoom (right). In the left panel, the solid line is the exact solution and the symbols are the cell averages of the numerical solution. 4.5.1 The weight function Let ϕ(x) be a positive bounded function, which can be taken as a weight function. For any function q ∈ Hh1 , we define the weighted L2 -norm as µZ ¶ 12 2 kqkϕ,D = q ϕdx D in the domain D. If ϕ = 1 or D = Ω, the corresponding subscript will be omitted. We will consider two weight functions ϕ1 (x, t) and ϕ−1 (x, t), respectively, to determine the left-hand and right-hand boundary of the region RT such that, outside this region, we can resume the (k +1)-th order accuracy in the L2 -norm. Both weight functions are related to the cut-off of the exponent function φ(r) ∈ C 1 : Ω → R,   2 − er , r < 0, φ(r) =  e−r , r > 0, 69 and they are defined as the solutions of the linear hyperbolic problem, ϕat + ϕax = 0, (4.33) µ ¶ a a(x − xc ) ϕ (x, 0) = φ , (4.34) γhσ where γ > 0, 0 < σ < 1 and xc are three parameters which will be chosen later. We always assume γhσ−1 ≥ 1 in this section. In [117], the authors have listed several properties about the two weight functions. Here, we state those that will be used. Proposition 4.5.1. For each weight function ϕa (x, t), the following properties hold 1 ≤ ϕa (x, t) ≤ 2, a(x − xc − t) ≤ 0, (4.35) 0 < ϕa (x, t) < hs , a(x − xc − t) > s log(1/h)γhσ . (4.36) Lemma 4.5.1. Let V be a Gauss-Radau projection, either P− or P+ . For any sufficiently smooth function p(x), there exists a positive constant C independent of h and p, such that kV⊥ pkϕ,D ≤ Chk+1 k∂xk+1 pkϕ,D , (4.37) kV⊥ (ϕvh )kϕ−1 ,D ≤ Cγ −1 h1−σ kvh kϕ,D , (4.38) kV(ϕvh )kϕ−1 ,D ≤ Ckvh kϕ,D . (4.39) where D is either the single cell Ij or the whole computational domain Ω. Lemma 4.5.2. For any function v ∈ Vh there holds the following identity 1X 1 H(v, ϕv) = − ϕj+ 1 [v]2j+ 1 + (v, ϕx v). (4.40) 2 j 2 2 2 70 4.5.2 The smooth solution We consider the following problem vt + vx = 0, (4.41) v(x, 0) = v0 (x), (4.42) where the initial condition v0 (x), is a sufficiently smooth function modified from the original initial condition u0 (x) = δ(x) + f (x) such that it agrees with u0 (x) for all x ∈ Ω\Ii , and satisfies |∂xα v0 (x)| ≤ Ch−α−1 , x ∈ Ii , where Ii is the cell containing x = 0. 4.5.3 Error representation and error equations Denote the error by e = v −uh , where uh approximates the solution to (4.1) or (4.41). Clearly, e also satisfies (2.5). We divide the error into the form e = η − ξ, where η = v − P− v = P⊥ − v, and ξ = uh − P− v. Following [117], we obtain dkξk2ϕ ¡ ¢ = 2 ξt , P⊥ + (ϕξ) + 2 (ξt , P+ (ϕξ)) − (ξ, ϕx ξ) dt ¡ ¢ = 2 ξt , P⊥ + (ϕξ) + 2 (ηt , P+ (ϕξ)) − 2 (et , P+ (ϕξ)) − (ξ, ϕx ξ) ¡ ¢ = 2 ξt , P⊥ + (ϕξ) + 2 (ηt , P+ (ϕξ)) − 2H(e, P+ (ϕξ)) − (ξ, ϕx ξ) ¡ ¢ = 2 ξt , P⊥ + (ϕξ) + 2 (ηt , P+ (ϕξ)) + 2H(ξ, P+ (ϕξ)) − (ξ, ϕx ξ) ¡ ¢ = 2 ξt , P⊥ + (ϕξ) + 2 (ηt , P+ (ϕξ)) + 2H(ξ, ϕξ) − (ξ, ϕx ξ) = 2Π1 + 2Π2 − Π3 , 71 where ¡ ¢ X Π1 = ξt , P⊥ + (ϕξ) , Π2 = (ηt , P+ (ϕξ)) , Π3 = ϕj+ 1 [ξ]2j+ 1 . 2 2 j First we estimate Π1 . Denote w = ξt − Pk−1 ξt . From the scheme (2.20), we have + (ξt , w)j = (ηt , w)j − (et , w)j = (ηt , w)j − [ξ]j− 1 wj− 1. 2 2 √ Inserting this into Π1 and defining ψ = ϕ, we obtain à ! (ξt , w)j (ξt , P⊥ + (ϕξ))j = 2 w, P⊥ + (ϕξ) kwkIj j à ! ³ w ´ + = (ηt , w)j − [ξ]j−1/2 wj− , P⊥ (ϕξ) kwk2Ij + 1 2 j C ³ ¯ ¯´ ° ° ¯ + ¯ °ψ −1 P⊥ ° ≤ |(ψηt , w)j | + ¯[ψξ]j−1/2 wj− 1¯ + (ϕξ) Ij kwkIj 2 Ch1−σ ³ ´ Ch1/2−σ ³ ´ ≤ kηt k2ϕ,Ij + kξk2ϕ,Ij + ϕj−1/2 [ξ]2j−1/2 + kξk2ϕ,Ij . γ γ Summing up with respect to j, we obtain à ! Ch1−σ ¡ ¢ Ch1/2−σ X (ξt , P⊥ + (ϕξ)) ≤ kη k2 t ϕ + kξk2 ϕ + ϕj−1/2 [ξ]2j−1/2 + kξk2ϕ . γ γ j For Π2 , it is not difficult to see that Π2 ≤ Ckηt kϕ kξkϕ ≤ C(kηt k2ϕ + kξk2ϕ ). Then if γ is large enough and σ = 12 , we have ¡ ¢ 2Π1 + 2Π2 − Π3 ≤ C kηt k2ϕ + kξk2ϕ . 72 By Gronwall’s inequality, Z T kξ(T )k2ϕ ≤C kηt k2ϕ dt + Ckξ(0)k2ϕ . (4.43) 0 4.5.4 The final estimate This part is almost the same as in [117]. We will only discuss the left-hand boundary of RT since the discussion for the right one is similar. Denote xL (t) = t + xc with xc = −2s log(1/h)γhσ , where s and γ are sufficiently large and σ = 1/2. As mentioned before, the δ- singularity in the initial datum is assumed to be contained in cell Ii . By proposition 4.5.1, we obtain 0 < φ(x) < hs for any x ∈ Ii . We choose v0 to satisfy Pk v0 = Pk u0 = uh (0). Then kξ(0)kϕ ≤ kξ(0)kϕ,L2 (R\Ii ) + kξ(0)kϕ,L2 (Ii ) ≤ Chk+1 kf kk+2 + Chs−1/2 . If s is large enough, then kξ(0)kϕ ≤ Chk+1 . Define the domain R+ T = (xL (T ), ∞), then kuh − vkR\R+ ≤ kuh − vkϕ,R\R+ ≤ kηkϕ,R\R+ + kξkϕ ≤ Chk+1 kf kk+1 + kξkϕ . T T T To estimate the second term on the right hand side, we use (4.43). Denote 1 w(t) = max{xj+ 1 : xj− 1 < t + xc , ∀j}, 2 2 2 and R1 (t) = (−∞, w(t)), R2 (t) = R\R1 (t) = (w(t), ∞). If γhσ−1 is large enough, R1 (t) stays away from the bad interval [t − h, t + h] where v(x, t) 6= u(x, t), then we 73 have kηt kϕ,R1 (t) ≤ Chk+1 kf kk+2 . Now we proceed to estimate kηt kϕ,R2 (t) . Since R2 contains the entire bad region, we will use the property of the weight function. By (4.36) we have ϕ ≤ hs in this zone. Then we obtain kηt kϕ,R2 (t) ≤ Chs/2 kηt kR2 (t) ≤ Chs/2+k+1 k∂xk+2 vkR2 (t) ≤ Ch(s−3)/2 +Chs/2+k+1 kf kk+2,R2 (t) . Similarly, we can estimate the right-hand side of the non-smooth region. If we take s large enough, we have kuh − u(x, T )kR\R+ = kuh − v(x, T )kR\R+ ≤ Chk+1 kf kk+2 + Ch(s−3)/2 ≤ Chk+1 . T T 4.6 Concluding remarks In this chapter, we use a DG method to solve linear hyperbolic conservation laws involving δ-singularities. We investigate the negative-order norm error estimates for the accuracy of the DG approximations to linear hyperbolic conservation laws with singular initial data or singular source terms, and obtain error estimates in the L2 -norm after post-processing in one space dimension. Numerical experiments demonstrate that the estimates are optimal. The results in this chapter offers evi- dence that the DG method is a good algorithm for problems involving δ-singularities in their solutions. Chapter 5 Applications to Krause consensus models and pressureless Euler equations In this chapter, we apply DG methods to solve hyperbolic conservation law ut + f (u)x = 0, (x, t) ∈ R × (0, T ], (5.1) u(x, 0) = u0 (x), x ∈ R, and its two dimensional version, where the exact solution u(x, t) contains δ-singularities. We extend the work in Chapter 4 and consider applications to two model equations: Krause’s consensus models and the pressureless Euler equations. These two models are used to describe the collision of particles, and the distributions can be identified as density functions. If the particles are places at a single point, the density func- tion is a δ-function and is difficult to approximate numerically. Recently, in [119], genuinely maximum-principle-satisfying high order DG schemes for scalar equations and two-dimensional incompressible flows in vorticity-streamfunction formulation have been constructed. Subsequently, positivity-preserving high order DG schemes for compressible Euler equations were given in [120]. We will extend the ideas in 74 75 [119, 120] to construct bound-preserving high order DG schemes for the Krause’s consensus models and pressureless Euler equations. For the first example, we discuss the following Krause’s consensus model equation ρt + Fx = 0, x ∈ [0, 1], t > 0, (5.2) ρ(x, 0) = ρ0 (x), t > 0, where ρ is the density function, which is always positive. The flux F is given by F (x, t) = v(x, t)ρ(x, t), and the velocity v is defined by Z 1 v(x, t) = (y − x)ξ(y − x)ρ(y, t)dy, 0 where 0 ≤ ξ(x) ≤ 1 is supported on a ball centered at zero with radius R. In [19], Canuto et al. investigated the discretized version of the PDE and proved that when the time t tends to infinity, the density function ρ will converge to some isolated δ- singularities, and the distance between any two of these δ-singularities cannot be less than R. Some computational results are shown in [19] based on a first order finite volume method. For two dimensions, if the initial density is rotationally invariant, the limit density should also be rotationally invariant, and hence can only be a single δ located at the center. However, direct computations on rectanglular meshes yield more than one δ singularity for sufficiently small R as a resultof the meshes not being invariant under rotation. In this chapter, following ideas in [23, 24], we construct a special mesh to obtain symmetry and convergence to physically relevant solutions. Computational results are given to demonstrate the advantages of high order DG schemes. 76 For the second example, we discuss the pressureless Euler equation wt + f (w)x = 0, t > 0, x ∈ R, (5.3)     ρ m w= , f (w) =  , m ρu2 with m = ρu, where ρ is the density function and u is the velocity. Pressureless Euler equations in one space dimension have been analyzed at the theoretical level intensively, e.g. [13, 14, 18, 22, 41]. Some numerical methods have also been studied by several authors [15, 10, 27, 16]. In [15], only first and second order numerical schemes were considered. Except for those in [15], no other methods seem to have been designed to solve (5.3) directly. In [16], the authors added an artificial vis- cosity and built a diffusive scheme. In [27], the authors applied the sticky particle methods to the equation and showed that the approximation satisfies the original system within a certain residual. In [10], the authors introduced a new variable and added one more equation to the system, leading to more computational cost. In this chapter, we consider high order DG scheme and approximate the equation without modification. Physically, the density ρ is positive and the velocity u satisfies a max- imum principle. We extend the idea in [120] and construct suitable limiters to fulfill these two requirements while maintaining high order accuracy. Moreover, numerical evidences demonstrate that the scheme is good for approximations in the presence of vacuum. Finally, our scheme works well in two dimensions. To the authors’ knowl- edge, not too much works in the literature focus on the two dimensional equations, and some theoretical results can be found in [96, 97]. However, complete existence and uniqueness results are not available. Therefore, our scheme offers a good tool to study two-dimensional pressureless Euler equations and other similar equations. 77 5.1 Preliminaries 5.1.1 Limiters In this subsection, we use forward Euler for time discretization and briefly discuss the construction of bound-preserving limiters [121]. We denote unj and unj to be the numerical solution and its cell average at time level n in cell Ij . If we consider generic numerical solution on the whole computational domain Ω, then the subscript j will be omitted. Suppose the exact solution of (5.1) is in some convex set G and we are interested in constructing numerical solutions which are also in G. The whole precess can be divided into three steps. In the first step, we consider a first order scheme ³ ¡ ¢ ¡ n n ¢´ un+1 = u n j + λ ˆ f u n j−1 , u n j − ˆf uj , uj+1 j 1³ n ¡ ¢´ 1 ³ n ¡ ¢´ = uj + 2λˆf unj−1 , unj + uj − 2λˆ f unj , unj+1 2 2 1 ¡ n ¢ 1 ¡ ¢ = H1 uj−1 , unj , 2λ + H2 unj , unj+1 , 2λ , (5.4) 2 2 where H1 (u, v, c) = v + c ˆ f (u, v), H2 (u, v, c) = u − c ˆ f (u, v). (5.5) ∆t Here unj = unj is a constant in each cell Ij , and λ = ∆x is the ratio of time and space mesh sizes. For many two-point first order numerical fluxes, we can prove the following property. Property 5.1.1. Suppose G is a convex set and u, v ∈ G, then there exists a positive constant C? , such that, for any 0 < c < C? , we have H1 (u, v, c) , H2 (u, v, c) ∈ G. Based on the above property, we can easily obtain that, under the CFL condition C? λ< 2 , un ∈ G implies un+1 ∈ G. Next, we study high order schemes and assume un ∈ G. Consider the equation 78 satisfied by the numerical cell averages ³ ´ ¯ n+1 u = ¯ nj u ˆ − + ˆ − + + λ f (uj− 1 , uj− 1 ) − f (uj+ 1 , uj+ 1 ) . (5.6) j 2 2 2 2 Let αi be the Legendre Gauss-Lobatto quadrature weights for the interval [− 12 , 12 ] P such that M i=0 αi = 1, with 2M − 3 ≥ k, and denote the corresponding Gauss- xji }. Then the Gauss-Lobatto quadrature yields Lobatto points in cell Ij as {ˇ M X ¯ nj = u xji ). αi unj (ˇ i=0 xj0 ) = u+ Clearly, unj (ˇ j− 1 xjM ) = u− and unj (ˇ j+ 1 . Therefore, 2 2 M X ³ ´ ¯ n+1 u = xji ) αi unj (ˇ ˆ − + ˆ − + + λ f (uj− 1 , uj− 1 ) − f (uj+ 1 , uj+ 1 ) j 2 2 2 2 i=0 M X −1 µ ¶ µ ¶ λ λ = xji ) αi unj (ˇ + α0 H1 u− j− 12 , u+ j− 12 , + αM H2 u− j+ 21 , u+ j+ 12 , . i=1 α0 αM ³ ´ If the numerical flux satisfies Property 5.1.1, we have H1 ∈ G u− j− 21 , u+ , λ j− 12 α0 ³ ´ and H2 u− j+ 1 , u+ , λ ∈ G, provided the suitable CFL condition λ < α0 C? is j+ 1 αM 2 2 xji ) ∈ G and G is a convex satisfied. Here, we use the fact that α0 = αM . Since unj (ˇ ¯ n+1 set, we have u j ∈ G. Finally, we can modify the numerical solution through the simple scaling limiter ¡ ¢ ˜ n+1 u j =u ¯ n+1 j + θ un+1 j −u¯ n+1 j ˜ n+1 . By taking suitable θ ∈ [0, 1], we have u j ∈ G, ˜ n+1 and u j is used as the numerical solution at time level n + 1. In many situations we can prove that this modification does not affect the high order accuracy of the original solution un+1 j [121]. 79 5.1.2 High order time discretizations We will use strong stability preserving (SSP) high order time discretizations to solve the ODE system ut = Lu. More details of these time discretizations can be found in [102, 101, 46]. In this chapter, we use the third order SSP Runge-Kutta method [102] u(1) = un + ∆tL(un ), 3 1 ¡ (1) ¢ u(2) = un + u + ∆tL(u(1) ) , (5.7) 4 4 1 2 ¡ (2) ¢ un+1 = un + u + ∆tL(u(2) ) , 3 3 and the third order SSP multi-step method [101] µ ¶ n+1 16 n 11 12 u = (u + 3∆tL(un )) + u n−3 n−3 + ∆tL(u ) . (5.8) 27 27 11 Since a SSP time discretization is a convex combination of the forward Euler, by using the limiter mentioned in Section 5.1.1, the numerical solution obtained from the full scheme is also in G. 5.2 Krause’s consensus models In this section we apply DG methods to Krause’s consensus models, extending the results in [19]. 5.2.1 Positivity-preserving high order schemes We consider (5.2) in more detail. For this model, we define G = {ρ : ρ > 0}. Clearly, G is a convex set. We start with the following first order scheme ³ ´ ρn+1 j = ρnj +λ vj− 1 h(ρnj−1 , ρnj ) − vj+ 1 h(ρnj , ρnj+1 ) , (5.9) 2 2 80 where h(·, ·) is a numerical flux, and ρnj = ρnj is the numerical approximation to the exact solution in cell Ij at time level n, with ρnj being its cell average. Moreover, vj− 1 is the numerical velocity at the interface xj− 1 , given by 2 2 N Z X vj− 1 = (y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy. 2 2 2 i=1 Ii For this model, we use an upwind flux, i.e.   vu v ≥ 0, vh(u, w) =  vw v < 0, and define H1 and H2 as H1 (u, w, c) = w + cvh(u, w), H2 (u, w, c) = u − cvh(u, w). Then the scheme (5.9) can be written as 1 1 ρn+1 j = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1 , 2λ). 2 2 Based on the above notations, we can prove the following lemma. Lemma 5.2.1. Suppose ρn > 0, then under the CFL condition 1 max |vj− 1 |λ < , j 2 2 we have ρn+1 > 0. Proof: We consider H1 (ρnj−1 , ρnj , 2λ) first. If vj− 1 < 0, then 2 H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = (1 + 2λvj− 1 )ρnj > 0. 2 2 81 On the other hand, if vj− 1 ≥ 0, then 2 H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = ρnj + 2λvj− 1 ρnj−1 > 0. 2 2 Similarly, we have H2 (ρnj , ρnj+1 , 2λ) = ρnj − 2λvj+ 1 h(ρnj , ρnj+1 ) > 0. 2 Therefore, 1 1 ρn+1 j = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1 , 2λ) > 0. 2 2 This completes the proof. Now, we consider high order schemes. The analysis in Section 5.1.1 implies the following theorem. Theorem 5.2.1. Suppose the DG solution ρn > 0, then under the CFL condition max |vj− 1 |λ < α0 , j 2 we have ρ¯n+1 > 0. Based on the above theorem, we can modify the density ρnj in the following steps. © ª • Set up a small number ε = min 10−13 , ρ¯nj . xji ), where {ˇ • Compute mj = mini ρnj (ˇ xji } are the Gauss-Lobatto points in cell Ij . • If mj < ε, then we take ρnj − ε θ= n , ρj − m j and use ρenj = ρnj + θ(ρnj − ρnj ) (5.10) 82 as the DG approximation in cell Ij at time level n. Based on the above steps, the numerical density is always positive. Therefore Z Z n n kρ kL1 (Ω) = ρ (x)dx = ρ0 (x)dx = kρ0 kL1 (Ω) , (5.11) Ω Ω where kukL1 (Ω) is the standard L1 -norm of u on Ω. Clearly, (5.11) implies the L1 stability for the DG scheme. Moreover, we can also derive a sufficient CFL condition which does not depend on the numerical velocity in Theorem 5.2.1. Actually, Note that N Z X vj− 1 = (y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy ≤ Rkρ0 kL1 (Ω) , 2 2 2 i=1 Ii such that the sufficient CFL condition is α0 λ≤ . (5.12) Rkρ0 kL1 (Ω) To summarize, we have the following theorem. Theorem 5.2.2. Under the CFL condition (5.12), the DG scheme with the positivity- preserving limiter for equation (5.2) is L1 stable and the density function is always positive. 5.2.2 Numerical experiments In this subsection, some numerical examples will be given to demonstrate the good performance of the DG scheme. Example 5.1. We consider the following problem ρt + (vρ)x = 0, x ∈ [0, 1], t > 0, (5.13) ρ(x, 0) = ρ0 (x), t > 0, 83 where the velocity v is defined by Z x+R v(x, t) = (y − x)ρ(y, t)dy. x−R We apply the positivity-preserving limiter and use P 0 and P 1 polynomials. Moreover, we use the third order SSP Runge-Kutta discretization in time [102] with ∆t = 0.1∆x. Figure 5.1 shows the numerical approximations of ρ(x) at t = 1000, with 30 30 25 25 20 20 15 15 10 10 5 5 0 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 X X Figure 5.1: Numerical density for (5.13) at t = 1000 with N = 400 when using P 0 (left) and P 1 (right) polynomials. N = 400, ρ0 = 1, and R = 0.02. We can observe 22 δ-singularities in each panel, and the distance between any two adjacent singularities is greater than R. The algorithm is quite stable in this simulation. Moreover, the P 1 solution in the right panel is more accurate than the P 0 one in the left panel, since the heights of the δ-singularities are almost doubled. Example 5.2. We consider the model problem in two dimensions. ρt + div(vρ) = 0, x ∈ [−1, 1]2 , t > 0, (5.14) ρ(x, 0) = ρ0 (x), t > 0, 84 where the velocity v is defined by Z v(x, t) = (y − x)ρ(y, t)dy. BR (x) In this example, we take R = 0.1 and   1 r < 0.5, ρ0 (x) =  0 r > 0.5, where r = kxk is the Euclidean norm of x. In [19], the authors demonstrated that the exact solution should be a single delta placed at the origin. However, by using rectangle meshes, we observe more than one delta singularity for R sufficiently small as a consequence of the meshes not being invariant under rotation. To address this problem, we follow [23, 24], and construct a special equal-angle-zoned mesh. The structure of the mesh is given in Figure 5.2. By using this mesh, the limit density given in Figure 5.3 is a single delta placed at the origin. Figure 5.2: Equal-angle-zoned mesh. 85 Figure 5.3: Numerical density ρ for (5.14) at t = 2000 with N = 200 when using P 0 polynomials. 5.3 Pressureless Euler equations In this section, we apply DG methods to the pressureless Euler equations. 5.3.1 Numerical schemes in one dimension We study (5.3) in more detail. Physically, the density is positive and the velocity satisfies the maximum principle. Therefore, we define      ρ  G = w =   : ρ > 0, aρ ≤ m ≤ bρ ,  m  where a = min u0 (x), b = max u0 (x), (5.15) 86 with u0 being the initial velocity. Clearly, G is a convex set. As mentioned in Section 5.1.1, we start with the following first order scheme, ¡ ¢ wjn+1 = wjn + λ h(wj−1 n , wjn ) − h(wjn , wj+1 n ) , (5.16) ¡ ¢T where h(·, ·) is a numerical flux and wjn = ρnj , mnj is the numerical approximation ¡ ¢T to the exact solution in cell Ij at time level n. Moreover, we define wnj = ρnj , mnj as its cell average. Clearly, for a first order scheme, wjn = wnj in (5.16). For simplicity, mn we use unj for ρn j as the numerical velocity throughout this section. In this problem, j we consider the Godunov flux [15]. Suppose that at the cell interface x = xj− 1 we 2 T T have two numerical approximations w` = (ρ` , m` ) and wr = (ρr , mr ) from the left and right respectively. Then the Godunov flux is given as    (m` , ρ` u2` ) u` > 0, ur > 0,       (0, 0) u` ≤ 0, ur > 0,      (m , ρ u2 ) u` ≤ 0, ur ≤ 0, r r r T d 2 ρuj− 1 , ρu j− 1 ) = (h (w` , wr )) = (c 2 2   (m` , ρ` u2` ) u` > 0, ur ≤ 0, v > 0,       (mr , ρr u2r ) u` > 0, ur ≤ 0, v < 0,      ( m` +mr , ρ u2 = ρ u2 ) u` > 0, ur ≤ 0, v = 0, 2 ` ` r r (5.17) where √ √ m` mr ρ` u` + ρr ur u` = , ur = , and v = √ √ . ρ` ρr ρ` + ρr For this problem, H1 and H2 are taken to be H1 (u, v, c) = v + c h (u, v) , H2 (u, v, c) = u − c h (u, v) . (5.18) Clearly, (5.16) can be written as 1 ¡ n ¢ 1 ¡ ¢ wjn+1 = H1 wj−1 , wjn , 2λ + H2 wjn , wj+1 n , 2λ . 2 2 87 Before proceeding to the theoretical results for the scheme, we would like to introduce the following lemma. The proof is trivial and is omitted. Lemma 5.3.1. Suppose {ˇ xi } are positive real numbers, and a ≤ yˇi ≤ b, ∀i, then Pn xˇi yˇi a ≤ Pi=1 n ≤ b. i=1 xˇi We will use Lemma 5.3.1 to prove the following lemma. Lemma 5.3.2. Suppose wn ∈ G, then under the CFL condition 1 λ< , 2 max(|a|, |b|) where a and b are defined in (5.15), we have wn+1 ∈ G. ¡ n ¢ ¡ ¢ Proof: We will only prove H1 wj−1 , wjn , 2λ ∈ G. The proof for H2 wjn , wj+1n , 2λ ∈ ¡ n ¢ G follows the same lines. Define H1 wj−1 , wjn , 2λ = (ˇ ˇ T , then the velocity de- ρ, m) rived from H1 is given as d2 1 mnj + 2λρu m ˇ j− 2 uˇ = = n . (5.19) ρˇ ρj + 2λc ρuj− 1 2 We will prove ρˇ > 0 and a ≤ uˇ ≤ b. To do so, we have to determine what {ˇ xi } and {ˇ yi } should be in Lemma 5.3.1, by testing the different choices for the numerical flux d2 1 ρu j− in (5.17). For simplicity, we define u b j− 12 = ρc 2 uj− 1 , if ρc uj− 1 6= 0. 2 2 uj− 1 = mnj−1 , then u • If ρc bj− 1 = unj−1 > 0. We take xˇ1 = ρnj , yˇ1 = unj and 2 2 xˇ2 = 2λmnj−1 , yˇ2 = unj−1 . uj− 1 = mnj , then u • If ρc bj− 1 = unj ≤ 0. We take xˇ1 = ρnj + 2λmnj , yˇ1 = unj . 2 2 d2 1 = (mn un + mn un )/2. We combine uj− 1 = (mnj−1 + mnj )/2, then ρu • If ρc 2 j− j−1 j−1 j j2 the two situations above, and take xˇ1 = ρnj + λmnj , yˇ1 = unj and xˇ2 = λmnj−1 , yˇ2 = unj−1 . 88 • If ρc d2 1 = 0. We take xˇ1 = ρn , yˇ1 = un . uj− 1 = 0, then ρu 2 j− 2 j j Clearly, in each case, xˇi > 0, therefore ρˇ > 0. Moreover, we observe a ≤ yˇi ≤ b, so that by Lemma 5.3.1, we have a ≤ uˇ ≤ b. Now, we consider high order schemes. By the same analysis as in Section 5.1.1, we have the following result. Theorem 5.3.1. Suppose wn ∈ G, then under the CFL condition α0 λ< , max(|a|, |b|) we have wn+1 ∈ G. Based on this, we can modify the numerical solution wjn while keeping the cell average untouched. Due to rounding error, we define      ρ m  Gε = w =   : ρ ≥ ε, a − ε ≤ ≤b+ε ,  m ρ       ρ m  ∂Gε = w =   : ρ ≥ ε, = a − ε or b + ε .  m ρ  Then the modification of wjn is given in the following steps. • Set up a small number ε = 10−13 . • If ρ¯nj > ε, then proceed to the following steps. Otherwise, ρnj is identified as the approximation to vacuum, and the velocity is undefined. Therefore, we take e jn = wnj as the numerical solution and skip the following steps. w xji ), where {ˇ • Modify the density first: Compute mj = mini ρnj (ˇ xji } are the Gauss- Lobatto points in cell Ij , and get ρenj by (5.10). Then use ρenj as the new numerical density ρnj . 89 • Modify the velocity: Define qji = wjn (ˇ xji ) in cell Ij . If qji ∈ Gε , then take θij = 1. Otherwise, take ° n ° °wj − sj ° i θij = ° n °, °wj − qj ° i where k·k is the Euclidean norm, and sji is the intersection point of the straight line s(t) = (1 − t)wnj + tqji , 0 ≤ t ≤ 1, and the surface ∂Gε . Define θj = mini=0,··· ,m θij , and use e jn = wnj + θj (wjn − wnj ), w as the DG approximation in cell Ij . 5.3.2 Numerical schemes in two dimensions We extend our work to two dimensions and consider the following equation wt + f (w)x + g(w)y = 0, t > 0, (x, y) ∈ R2 , (5.20)       ρ m n             w = m, f (w) =  ρu2  , g(w) =  ρuv  ,       n ρuv ρv 2 with m = ρu, n = ρv, where ρ is the density function and (u, v) is the velocity field. We define       ρ           2 2 2 2 G = w =  m  : ρ > 0, m + n ≤ S ρ ,           n 90 where ¡ ¢ S > 0, and S 2 = max u2 (x, y, 0) + v 2 (x, y, 0) x,y with (u, v)(x, y, 0) being the initial velocity field. Clearly, G is a convex set. For simplicity, we use uniform rectangular meshes. The cell is defined as Iij = h i h i xi− 1 , xi+ 1 × yj− 1 , yj+ 1 , and the mesh sizes in x and y directions are denoted as 2 2 2 2 ∆x and ∆y respectively. At time level n, we approximate the exact solution with a n vector of polynomials of degree k, wij = (ρnij , mnij , nnij )T , and define the cell average + − + − wnij = (ρnij , mnij , nnij )T . Moreover, we denote wi− 1 (y), w ,j i+ 1 ,j (y), wi,j− 1 (x), w i,j+ 1 (x) 2 2 2 2 as the traces of w on the four edges of cell Iij respectively. More details can be found mn nn in [120]. We use (unij , vijn ) for ( ρnij , ij ρn ) as the numerical velocity field in cell Iij at ij ij time level n, and define a1 = maxij |unij | and a2 = maxij |vijn |. For simplicity, if we consider a generic numerical solution on the whole computational domain at time level n, the subscript ij will be omitted. We only consider high order schemes, and the one satisfied by the cell averages can be written as Z y 1 ³ ´ ³ ´ ∆t j+ 2 wn+1 ij = wnij + h1 wi−− 1 (y), w + 1 (y) − h 1 w − 1 (y), w + 1 (y) dy ∆x∆y y 1 2 ,j i− 2 ,j i+ 2 ,j i+ 2 ,j j− 2 Z x 1 ³ ´ ³ ´ ∆t i+ 2 − + − + + h2 wi,j− 1 (x), w i,j− 12 (x) − h 2 wi,j+ 12 (x), w i,j+ 21 (x) dx,(5.21) ∆x∆y x 1 2 i− 2 where h1 (·, ·) and h2 (·, ·) are one-dimensional numerical fluxes. For this problem, we also use the Godunov flux. Suppose (x, y) = (xi− 1 , y0 ) is a point on the vertical 2 cell interface, at which we have two numerical approximations w` = (ρ` , m` , n` )T and wr = (ρr , mr , nr )T from left and right respectively. Then the Godunov flux 91 (h1 (w` , wr ))T can be written as    (m` , ρ` u2` , ρ` u` v` ) u` > 0, ur > 0,       (0, 0, 0) u` ≤ 0, ur > 0,    ³ ´  (m , ρ u2 , ρ u v ) r r r r r r u` ≤ 0, ur ≤ 0, ρc d 2 u, ρu , ρd uv =   (m` , ρ` u2` , ρ` u` v` ) u` > 0, ur ≤ 0, v > 0,       (mr , ρr u2r , ρr ur vr ) u` > 0, ur ≤ 0, v < 0,      1 (m + m , ρ u2 + ρ u2 , m v + m v ) u` > 0, ur ≤ 0, v = 0, 2 ` r ` ` r r ` ` r r where µ ¶ µ ¶ √ √ m` n` mr n` ρ` u` + ρr ur (u` , v` ) = , , (ur , vr ) = , , and v = √ √ . ρ` ρ` ρr ρ` ρ` + ρr ³ ´T The numerical flux h2 = ρc v, ρd c2 uv, ρv can be defined in a similar way on the horizontal cell interfaces. For accuracy, we use L-point Gauss quadratures with L ≥ k + 1 to approximate the integrals in (5.21). More details of this requirement can be found in [32]. The h i h i Gauss quadrature points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 are denoted by 2 2 2 2 n o n o β y β pˆxi = xˆi : β = 1, · · · , L and pˆj = yˆj : β = 1, · · · , L , respectively. Also, we denote wˆβ as the corresponding weights on the interval £ 1 1¤ − 2 , 2 . Following the same notation as in previous sections, we use © ª xαi : α = 0, · · · , M } and pˇyj = yˇjα : α = 0, · · · , M pˇxi = {ˇ h i h i as the Gauss-Lobatto points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 respectively. Also, we 2 2 2 £ 12 1 ¤ denote wˇα as the corresponding weights on the interval − 2 , 2 . 92 ∆t ∆t Let λ1 = ∆x and λ2 = ∆y , then the numerical scheme (5.21) becomes L X h ³ ´ ³ ´i wn+1 ij = wnij + λ1 − wˆβ h1 wi− 1 ,β , w + i− 1 ,β − h 1 w − i+ 1 ,β , w + i+ 1 ,β 2 2 2 2 β=1 L X h ³ ´ ³ ´i − + − + + λ2 wˆβ h2 wβ,j− 1,w β,j− 1 − h 2 wβ,j+ 1,w β,j+ 1 , (5.22) 2 2 2 2 β=1 − where wi− 1 ,β − = wi− y β ) is a point value in the Gauss quadrature. Likewise for 1 (ˆ ,j j 2 2 the other point values. As the general treatment, we rewrite the cell average on the right hand side as M X X L M X X L wnij = 1 wˇα wˆβ wαβ = 2 wˇα wˆβ wβα , α=0 β=1 α=0 β=1 1 where wαβ 2 and wβα n denote wij xαi , yˆjβ ) and wij (ˇ n xβi , yˇjα ) respectively. We extend the (ˆ definitions of H1 and H2 in (5.18) to two-dimensional problems and define H11 (u, v, c) = v + c h1 (u, v) , H12 (u, v, c) = u − c h1 (u, v) , H21 (u, v, c) = v + c h2 (u, v) , H22 (u, v, c) = u − c h2 (u, v) . Let µ = a1 λ1 + a2 λ2 , then scheme (5.22) can be written as ÃM −1 ! L X X ³ ´ ³ ´ wn+1 ij = C1 wˆβ 1 wˇα wαβ − + wˇ0 H11 wi− 1 ,β + , wi− 1 ,β − , µ1 + wˇM H12 wi+ 1 ,β + , wi+ 1 ,β , µ1 2 2 2 2 β=1 α=1 ÃM −1 ! L X X ³ ´ ³ ´ 2 − + − + + C2 wˆβ wˇα wβα + wˇ0 H21 wβ,j− 1,w β,j− 1 , µ2 + wˇM H22 wβ,i+ 1,w β,j+ 1 , µ2 , 2 2 2 2 β=1 α=1 where a1 λ 1 a2 λ 2 µ µ C1 = , C2 = , µ1 = , µ2 = . µ µ a1 wˇ0 a2 wˇ0 Now, we can state the main theorem. 93 Theorem 5.3.2. Suppose wn ∈ G, then under the CFL condition ∆t ∆t a1 + a2 ≤ wˆ0 , ∆x ∆y we have wn+1 ∈ G. ³ ´ − + Proof: For simplicity, we only prove H11 wi− 1 ,β , wi− 1 ,β , µ 1 ∈ G, ∀β, and define 2 2 ³ ´ m ˇ n ˇ − H11 wi− 1 , w + 1 , µ 1 = (ˇ ρ, m, ˇ nˇ )T , uˇ = , vˇ = . 2 ,β i− 2 ,β ρˇ ρˇ Following the same analysis as in Lemma 5.3.2, we have ρˇ > 0. Therefore, we need only prove uˇ2 + vˇ2 ≤ S 2 . By the assumption, we have ³ ´T ³ ´T − − − − + + + + wi− 1 ,β = ρ i− 1 ,β , m i− 1 ,β , n i− 1 ,β ∈ G and wi− 1 ,β = ρi− 1 ,β , m i− 1 ,β , n i− 1 ,β ∈ G. 2 2 2 2 2 2 2 2 ³ ´ ³ ´T − Denote h1 wi− , w + = ρ c u 1 , d2 1 , ρd ρu uv 1 as the corresponding 1 ,β 1 i− ,β i− ,β i− ,β i− ,β 2 2 2 2 2 numerical flux, and for any unit vector n = (n1 , n2 )T , define w ˇ = uˇn1 + vˇn2 . Then ³ ´ m+ n + n + n + µ d ρu 2 1 n +ρ d uv 1 n i− 1 ,β 1 i− 1 ,β 2 1 i− ,β 1 i− ,β 2 2 2 wˇ = 2 2 ρ+ i− 21 ,β + µ1 ρc ui− 1 ,β 2 ρ+ w+ i− 21 ,β i− 21 ,β + µ1 ρc ui− 1 ,β w bi− 1 ,β 2 2 = , ρ+ i− 21 ,β + µ1 ρc ui− 1 ,β 2 where m+ n + n+ i− 1 ,β 1 n i− 1 ,β 2 d2 1 n1 + ρd ρu i− ,β uv i− 1 ,β n2 + wi− 1 = 2 2 , w bi− 1 ,β = 2 2 . 2 ,β ρ+ i− 1 ,β 2 ρc ui− 1 ,β 2 2 + We can easily show that |wi− 1 ,β | ≤ S and |w bi− 1 ,β | ≤ S. Following the same lines as 2 2 is the proof of Lemma 5.3.2, we have |w| ˇ ≤ S. Choosing n to be parallel with (ˇ u, vˇ), we have uˇ2 + vˇ2 ≤ S 2 , completing the proof. 94 Remark 5.3.1. Since a1 ≤ S and a2 ≤ S, another sufficient CFL condition in ∆t ∆t w ˆ0 Theorem 5.3.2 is ∆x + ∆y ≤ S . n Based on the above theorem, we can modify the numerical solution wij keeping the cell average untouched. Due to the rounding error, we define       ρ         ε   2 2 2 2 G = w =  m  : ρ ≥ ε, m + n ≤ (S + ε) ρ ,           n       ρ         ε   2 2 2 2 ∂G = w =  m  : ρ ≥ ε, m + n = (S + ε) ρ .           n n Then the modification of wij is given in the following steps. • Set up a small number ε = 10−13 . • If ρnij > ε, then proceed to the following steps. Otherwise, ρnij is identified as the approximation to vacuum, and the velocity is undefined. Therefore, we n e ij take w = wnij as the numerical solution and skip the following steps. n o xαi , yˆjβ ), ρnij (ˆ • Modify the density first: Compute mij = minαβ ρnij (ˇ xβi , yˇjα ) . If mij < ε, then take ρenij as ¡ ¢ ρenij = ρnij + θij ρnij − ρnij , with ρnij − ε θij = , ρnij − mij and use ρenij as the new numerical density ρnij . 1 2 1 • Modify the velocity: Consider wαβ and wβα in cell Iij . If wαβ ∈ Gε , then take 95 1 θαβ = 1. Otherwise, take ° n ° °wij − s1 ° 1 αβ θαβ =° ° °wnij − w1 ° , αβ where k · k is the Euclidean norm, and s1αβ is the intersection point of the straight line 1 s1 (t) = (1 − t)wnij + twαβ , 0 ≤ t ≤ 1, and the surface ∂Gε . Similarly, we can define θβα 2 2 in the same way for wβα . Finally, we use n © 1 2 ª e ij w = wnij + θ(wij n − wnij ), θ = min θαβ , θβα , α,β as the DG approximation in cell Iij . 5.3.3 Numerical experiments Let us provide numerical experiments to demonstrate the good performance of the DG scheme for solving pressureless Euler equations. In all numerical simulations, if not otherwise stated, we use third order schemes and take N = 100. One space dimension We consider the problem in one space dimension and solve (5.3) with different initial conditions. Example 5.3. We consider the following initial data ρ0 (x) = sin(x) + 2, u0 (x) = sin(x) + 2, (5.23) with periodic boundary condition. 96 Clearly, the exact solution is ρ0 (x0 ) u(x, t) = u0 (x0 ), ρ(x, t) = , 1 + u00 (x0 ) where x0 is given implicitly by x0 + tu0 (x0 ) = x. We use the third order SSP multi-step method in time [101] with ∆t = 0.01∆x2 , and test the example by using P k polynomials with k = 1, 2, 3 on uniform meshes. Table 5.1 shows the L2 -norm of the error at t = 0.1. We observe (k + 0.5)-th order convergence. Table 5.1: L2 -norm of the error between the numerical density and the exact density for initial condition (5.23). k=1 k=2 k=3 N error order error order error order 20 1.41E-02 - 6.84E-04 - 3.40e-5 - 40 4.18E-03 1.76 1.04E-04 2.72 2.82e-6 3.59 80 1.30E-03 1.68 1.55E-05 2.74 2.26e-7 3.64 160 4.24E-04 1.62 2.41E-06 2.69 1.83e-8 3.62 320 1.51E-04 1.49 3.80E-07 2.67 1.49e-9 3.63 Example 5.4. We consider the following initial condition   1 x < 0,  1 x < 0, ρ0 (x) = u0 (x) = (5.24)  0.25 x > 0,  0 x > 0. Clearly, the exact solution is   (1, 1) x < 2t/3, (ρ(x, t), u(x, t)) =  (0.25, 0) x > 2t/3, 97 2t and at x = 3 , the density should be a δ-function. Figure 5.4 shows the numerical 1 15 0.8 0.6 10 u 0.4 5 0.2 0 0 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 x x Figure 5.4: Numerical density (left) and velocity (right) at t = 0.5 with P 1 polyno- mials for initial condition (5.24). density and velocity at t = 0.5 using P 1 polynomials. From the figure, we observe the numerical solution capture the profile of the exact solution quite well. Example 5.5. We consider the following initial condition    −0.5 x < −0.5,      0.4 −0.5 < x < 0, ρ0 (x) = 0.5, u0 (x) = (5.25)   0.4 − x 0 < x < 0.8,      −0.4 x > 0.8. The exact solution for t < 1 is    (0.5, −0.5) x < −0.5 − 0.5t,       (0, undefined) −0.5 − 0.5t < x < −0.5 + 0.4t,   (ρ(x, t), u(x, t)) = (0.5, 0.4) −0.5 + 0.4t < x < 0.4t,     0.5 0.4−x   ( 1−t , 1−t ) 0.4t < x < 0.8 − 0.4t,     (0.5, −0.4) x > 0.8 − 0.4t. 98 Figure 5.5 shows the numerical density and velocity at t = 0.5. From the figure, 0.4 1 0.2 0.8 0.6 0 u 0.4 -0.2 0.2 -0.4 0 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x x Figure 5.5: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (5.25). The solid line shows the exact solution while the symbols show the numerical solution. we can observe some local oscillations near the singularities. This is not surprising as we have not used any limiters other than the bound-preserving ones for the DG scheme. Example 5.6. We consider the following initial condition   −0.5 x < 0, ρ0 (x) = 0.5, u0 (x) = (5.26)  0.4 x > 0. The exact solution is    (0.5, −0.5) x < −0.5t,   (ρ(x, t), u(x, t)) = (0, undefined) −0.5t < x < 0.4t,     (0.5, 0.4) x > 0.4t. Figure 5.6 shows the numerical density and velocity at t = 0.5. From the figure, we can observe some local oscillations near the singularities. 99 0.5 0.4 0.4 0.2 0.3 0 u 0.2 -0.2 0.1 -0.4 0 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 x x Figure 5.6: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (5.26). Two dimensions We consider the problem in two dimensions and solve (5.20) with different initial conditions. Example 5.7. We consider the following initial condition ρ(x, y, 0) = ρ0 (x + y) = exp(sin(x + y)), u(x, y, 0) = u0 (x + y) = 31 (cos(x + y) + 2), (5.27) v(x, y, 0) = v0 (x + y) = 13 (sin(x + y) + 2). The exact solution is ρ0 (z0 ) u(x, y, t) = u0 (z0 ), v(x, y, t) = v0 (z0 ), ρ(x, y, t) = , 1 + u00 (z0 ) + v00 (z0 ) where z0 is given implicitly by z0 + t(u0 (z0 ) + v0 (z0 )) = x + y. 100 We use the third order SSP multi-step method in time [101] with ∆t = 0.01∆x3/2 , and test the example by using P k polynomials with k = 1, 2, 3. Table 5.2 shows the L2 -norm of the error at t = 0.1. From the table, we again observe about (k + 0.5)-th order convergence. Table 5.2: L2 -norm of the error between the numerical density and the exact density for initial condition (5.27). k=1 k=2 k=3 N error order error order error order 10 0.512 - 0.107 - 3.42E-02 - 20 0.176 1.54 3.12E-02 1.78 3.57E-03 3.26 40 6.48E-02 1.44 8.52E-03 1.87 4.86E-04 2.88 80 2.32E-02 1.48 1.39E-03 2.62 3.97E-05 3.61 160 9.08E-03 1.35 1.92E-04 2.86 3.65E-06 3.45 Example 5.8. We consider the following initial condition 1 1 1 ρ(x, y, 0) = , (u, v)(x, y, 0) = (− cos θ, − sin θ), (5.28) 100 10 10 where θ is the polar angle. Since all the particles are moving towards the origin, the density function at t > 0 should be a single delta at the origin. Different from Example 5.2 in Section 5.2.2, we can observe only one delta located at the origin by using rectangle mesh as shown in Figure 5.7. Example 5.9. We consider the following initial condition    (−0.25, −0.25) x > 0, y > 0,     1  (0.25, −0.25) x < 0, y > 0, ρ(x, y, 0) = , (u, v)(x, y, 0) = (5.29) 10   (0.25, 0.25) x < 0, y < 0,      (−0.25, 0.25) x > 0, y < 0. 101 0.4 0.2 0 y -0.2 -0.4 -0.4 -0.2 0 0.2 0.4 x Figure 5.7: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.28). Figure 5.8 shows the numerical density and velocity field at t = 0.5. From the figure, we can observe δ-singularities located at the origin and two axes. Example 5.10. We consider the following initial condition  1  (cos θ, sin θ) r < 0.3, ρ(x, y, 0) = , (u, v)(x, y, 0) = (5.30) 100  (− 1 cos θ, − 1 sin θ) r > 0.3, 2 2 p where r = x2 + y 2 and θ is the polar angle. Figure 5.9 shows the numerical density (contour plot) and velocity field at t = 0.5. From the figure, we can observe δ-shocks located on a circle and vacuum inside. Example 5.11. We consider the following initial condition    (0.3, 0.4) x > 0, y > 0,      (−0.4, 0.3) x < 0, y > 0, ρ(x, y, 0) = 0.5, (u, v)(x, y, 0) = (5.31)   (−0.3, −0.4) x < 0, y < 0,      (0.4, −0.3) x > 0, y < 0. 102 1 0.8 0.6 y 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 x Figure 5.8: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.29). 0.4 0.2 0 y -0.2 -0.4 -0.4 -0.2 0 0.2 0.4 x Figure 5.9: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (5.30). 103 1 0.8 0.6 y 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 x Figure 5.10: Numerical density (left) and velocity field (right) at t = 0.4 with N = 50 for initial condition (5.31). Figure 5.10 shows the numerical density (contour plot) and velocity field with N = 50 at t = 0.4. From the figure, we can observe that the numerical solution approximates the vacuum quite well. 5.4 Concluding remarks In this chapter, we developed DG methods to solve hyperbolic conservation laws in- volving δ-singularities. We study Krause’s consensus models and pressureless Euler equations to demonstrate the stability and high resolution of the DG approximations. Moreover, numerical experiments show that the scheme is also good for approxima- tions in the presence of vacuum. In future work we will extend DG methods to other equations involving δ-singularities to wider areas of applications. Part II Numerical cosmology 104 105 In this second part, we will study the Lyα photons transfer in an optically thick medium. Lyα photons have been widely applied to study the physics of luminous objects at various epochs of the universe, such as Lyα emitters, Lyα blob, damped Lyα system, Lyα forest, fluorescent Lyα emission, star-forming galaxies, quasars at high redshifts as well as optical afterglow of gamma ray bursts [49, 43, 38, 69]. The resonant scattering of Lyα photons with neutral hydrogen atoms has a profound effect on the time, space and frequency dependencies of Lyα photons transfer in an optically thick medium. Lyα photons emerging from an optically thick medium would carry rich information of the photon sources and halo surrounding the source of the Lyα photon. The profiles of the emission and absorption of the Lyα radiation are powerful tools to constrain the mass density, velocity, temperature and the fraction of neutral hydrogen of the optically thick medium. Radiation transfer of Lyα photons in an optically thick medium is fundamentally important. The radiative transfer of Lyα photons in a medium consisting of neutral hydrogen atoms has been extensively studied either analytically or numerically for more than half a century. In [3], the focus was on the numerical approximation of the redis- tribution function of resonant scattering, and no solution of the integro-differential equation of the radiative transfer has been found. Before that, Field [44] gave the first analytical solution of the integro-differential equation for the case of both medium and source uniformly distributed in the whole space. Analytical solutions of the fre- quency profile of photons emergent from optically thick halo are found based on the Fokker-Planck (F-P) approximation of the integro-differential equation [51, 76, 37]. Besides the F-P approximation, Monte Carlo (MC) simulations are also popular in solving the transfer of resonant photons (e.g. [75, 122, 107, 110, 70, 82, 114, 115]. However, many important topics cannot be seen with the above-mentioned solu- tions. Besides the Field’s analytical solution, all others are time-independent, and therefore miss the detailed balance relationship of resonant scattering [93] and can- not be used to describe the formation and evolution of the Wouthuysen-Field (W-F) 106 local thermalization of the Lyα photon frequency distribution [113, 44, 45], which is important for the emission and absorption of the hydrogen 21 cm line (e.g. [42]. The rich features of the Lyα photon transfer referring to the W-F local thermalization are fully missed. The F-P equation is based on the Eddington approximation, which assumes that the radiation intensity is a linear function of angular (direction) vari- able. Therefore, the solutions of the F-P equation do not provide the information of the evolution of the angular distribution of Lyα photons. Recently, a state-of-the-art numerical method has been introduced to solve the integro-differential equation of the radiative transfer with resonant scattering [83, 84, 85, 87]. The solver is based on the weighted essentially non-oscillatory (WENO) scheme [61]. With the WENO solver, many physical features of the transfer of Lyα photons in an optically thick medium [88, 89, 90], missed in the F-P equation ap- proximations, have been revealed. For instance, the WENO solution shows that the time scale of the formation of the W-F local thermal equilibrium is only about a few hundred times that of the resonant scattering. It also shows that the double peaked frequency profile of the Lyα photon, emerging from an optically thick medium, does not follow the time-independent solutions of the F-P equation [42, 87, 88, 89, 90]. These results indicate the needs of re-visiting problems which have been studied only via the F-P time-independent approximation. In this part, we will use a WENO solver to study the effect of dust and angular distribution of Lyα resonant photons transfer in an optically thick halo. Chapter 6 WENO solver of transfer equations of resonant photons 6.1 Basic theory 6.1.1 Radiative transfer equation of a dusty halo The radiative transfer equation of Lyα photons in a spherical halo with dust is given by ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, µ, x0 , µ0 ; a)I(η, r, x0 , µ0 )dx0 dµ0 /2 Z −κ(x)I + Aκ(x) Rd (x, x0 ; µ, µ0 ; a)I(η, r, x0 , µ0 )dx0 dµ0 + S, (6.1) where S is the source and I(t, rp , x, µ) is the specific intensity, which is a function of time t, radial coordinate rp , frequency x and the direction angle, µ = cos θ, with respect to the radial vector r. In (6.1), we use the dimensionless time η defined as η = cnHI σ0 t and the dimen- sionless radial coordinate r defined as r = nHI σ0 rp , where nHI is the numer density 107 108 of HI, and σ0 /π 1/2 is the cross section of HI resonent scatering of Lyα photons at resonant frequency ν0 = 2.46×1015 s−1 . That is, η and r are, respectively, in the units of mean free flight-time and mean free path of photon ν0 with respect to the resonant scattering without dust scattering and absorption. Without resonant scattering, a signal propagates in the radial direction with the speed of light and the orbit of the signal is r = η + const. φ(x, a) is the normalized Voigt profile [57] given as Z ∞ 2 a e−y φ(x, a) = dy . (6.2) π 3/2 −∞ (x − y)2 + a2 As usual, the photon frequency ν in (6.1) is described by the dimensionless frequency x ≡ (ν − ν0 )/∆νD , and ∆νD = ν0 (vT /c) = 1.06 × 1011 (T /104 )1/2 Hz is the Doppler p broadening by the thermal motion vT = 2kB T /m, T being the gas temperature of the halo. The parameter a in (6.2) is the ratio of the natural to the Doppler broadening. For the Lyα line, a = 4.7 × 10−4 (T /104 )−1/2 . The optical depth of Lyα photons with respect to HI resonant scattering is τs (x) = nHI σ(x)drp , where σ(x) = σ0 φ(x, a) is the cross section of scattering at ν, and therefore, the dimensionless size of the halo R is equal to the optical depth τ0 = nHI σ0 R. The re-distribution function R(x, µ, x0 , µ0 ; a) of (6.2), the derivation of which is given in Section 6.1.4, gives the probability of a photon absorbed at the frequency x0 direction µ0 , and re-emitted at the frequency x direction µ. It depends on the details of the scattering [54, 56, 58]. If we consider coherent scattering without recoil, the re-distribution function with the Voigt profile is Z Z " µ ¶ ¶2 #−1 µ 2π ∞ 0 0 a −u2 2 x + x0 (x − x0 )2 R(x, µ, x , µ ; a) = e a + exp − − αu dudφ, 0 −∞ 4π 3 β 4β 2 2 (6.3) p p q q where H = 1 − µ2 1 − µ02 cos φ + µµ0 , α = 1+H 2 , and β = 1−H 2 . In the case of a = 0, i.e., considering only the Doppler broadening, the re-distribution function 109 is Z · 2 ¸ 2π 0 0 1 2 − 2xx0 H + x02 R(x, µ, x , µ ) = √ exp − dφ, (6.4) 0 2π 2 1 − H 2 1 − H2 where H is exactly the same as in (6.3). The redistribution function of (6.4) is normalized as Z 1 Z ∞ 1 2 R(x, µ, x0 , µ0 )dx0 dµ0 = φ(x, 0) = π −1/2 e−x . 2 −1 −∞ With this normalization, the total number of photons is conserved in the evolution described by (6.1). That is, the destruction processes of Lyα photons, such as the two-photon process [105, 80], are ignored in (6.3). The recoil of atoms is not considered in (6.3) or (6.4). The absorption and scattering of dust are described by the term κ(x)I in (6.1), where κ(x) = σd /σ0 , which is of the order of 10−8 (T /104 )1/2 [40, 39]. The term with A in (6.1) describes albedo, i.e. A ≡ σs /σd , where σs is the cross section of dust scattering and σd is the effective cross-section per hydrogen atom, which describes the absorption and scattering of dust. Generally, A lies approximately between 0.3 and 0.4 [81, 112]. With dust, the optical depth is given by τ (x) = τ0 φ(x, a) + τd (x) (6.5) where the dust optical depth τd (x) = nHI σd (x)R. This is equal to assuming that dust is uniformly distributed in IGM. The effects of inhomogeneous density distributions of dust [77, 48] will not be studied. Since dust generally is much heavier than a single atoms, the recoil of dust par- ticles can be neglected when colliding with a photon. Under this “heavy dust” ap- proximation, photons do not change their frequency during the collision with dust. The redistribution function of dust Rd is independent of x and x0 , and is simply 110 given by a phase function as Z 2π ∞ X (2l + 1) d 0 1 0 1 − g2 R (µ, µ ) = dφ = g l Pl (µ)Pl (µ0 ), (6.6) 4π 0 (1 + g 2 − 2g µ ¯)3/2 l=0 2 p ¯ = µµ0 + where µ (1 − µ2 )(1 − µ02 )cosφ0 and Pl is the Legendre function. The factor g in (6.6) is the asymmetry parameter. For isotropic scattering, g = 0. The cases of g = +1 and -1 correspond to complete forward and backward scattering, respectively. Generally, the factor g is a function of the wavelength. For the Lyα photon, we will take g = 0.73 for realistic dust scattering [73]. The integral of (6.6) is performed in Section 6.1.5. In (6.1) we neglect the effect of collision transition from H(2p) state to H(2s) state, which can significantly affect the escape of Lyα photons when the HI column density is higher than 1021 cm−2 and dust absorption is very small [76]. This generally is out of the parameter range used below. We are also not considering the effects of bulk motion of the medium of halos (e.g. [104, 114]). In (6.1), the term with the parameter γ is due to the expansion of the universe. If −1 nH is equal to the mean of the number density of cosmic hydrogen, we have γ = τGP , and τGP is the Gunn-Peterson optical depth. Since the Gunn-Peterson optical depth is of the order of 106 at high redshift (e.g. [89]), the parameter γ is of the order of 10−5 − 10−6 . Therefore, if the optical depth of halos is less than or equal to 106 , the term with γ in (6.1) can be ignored. 6.1.2 Integrated redistribution function In general, it is difficult to solve the transfer equation for noncoherent scatter- ing. Therefore, the scattering is assumed to be isotropic and we need to integrate R(x, µ, x0 , µ0 ; a) over the angular direction. Denote R(x, x0 ; a) as the angular aver- aged re-distribution function. It gives the probability of a photon absorbed at the frequency x0 , and re-emitted at the frequency x. If we consider coherent scattering 111 without recoil, the angular-averaged re-distribution function with the Voigt profile can be written as, R(x, x0 ; a) = (6.7) Z ∞ · µ ¶ µ ¶¸ 1 −u2 −1 xmin + u −1 xmax − u e tan − tan du π 3/2 |x−x0 |/2 a a where xmin = min(x, x0 ) and xmax = max(x, x0 ). In the case of a = 0, i.e. including only the Doppler broadening, the re-distribution function is 1 R(x, x0 ) = erfc[max(|x|, |x0 |)]. (6.8) 2 R∞ The angular-averaged re-distribution function of (6.8) is normalized as −∞ R(x, x0 )dx0 = 2 φ(x, 0) = π −1/2 e−x . With this normalization, the total number of photons is also conserved in the evolution described by (6.1). 6.1.3 Eddington approximation Equation (6.6) indicates that the transfer equation (6.1) can be solved with the P Legendre expansion I(η, r, x, µ) = l Il (η, r, x)Pl (µ). If we take only the first two terms, l = 0 and 1, it is the Eddington approximation as I(η, r, x, µ) ' J(η, r, x) + 3µF (η, r, x) (6.9) where Z +1 Z +1 1 1 J(η, r, x) = I(η, r, x, µ)dµ, F (η, r, x) = µI(η, r, x, µ)dµ. (6.10) 2 −1 2 −1 112 They are, respectively, the angularly averaged specific intensity and flux. Defining j = r2 J and f = r2 F , (6.1) yields the equations of j and f as Z ∂j ∂f ∂j + = −(1 − A)κj − φ(x; a)j + R(x, x0 ; a)jdx0 + γ + r2 S, (6.11) ∂η ∂r ∂x ∂f 1 ∂j ∂f 2j + = −(1 − Ag)κf + γ − φ(x; a)f + . (6.12) ∂η 3 ∂r ∂x 3r The mean intensity j(η, r, x) describes the x photons trapped in the position r at time η by the resonant scattering, while the flux f (η, r, x) describes the photons in transit. For spherical halo with a central source, the term S of (6.1) can be replaced by a boundary condition of I(η, r, x, µ) at r = 0. If the angular distribution of photons is independent of photon’s frequency, we have r2 I(η, r, x, µ)|r→0 = S0 T (η)Θ(µ)φ(x). (6.13) where the functions T (η), Θ(µ), and φ(x) describe, respectively, the time-dependence, angular- and frequency-distributions of photons of the source. In this case, the source of (6.13) can be replaced by a boundary condition at r = 0 as Z 1 f (η, 0, x) = S0 T (η)φ(x) µΘ(µ)dµ = S0 φs (x). (6.14) 2 For example, if we take the boundary condition at r = 0 to be   6µπ −1/2 e−x2 , µ > 0, r2 I(η, r, x, µ)|r→0 = (6.15)  0, µ < 0. With (6.15) one can find I from (6.1), and then find j and f via (6.10). Therefore, the corresponded boundary condition of (6.14) is 2 f (η, r = 0, x) = π −1/2 e−x . (6.16) 113 We will use these two boundary conditions in Chapter 8. In (6.13) and (6.14), S0 is the intensity. Since (6.1), (6.11) and (6.12) are linear, the solutions of j(x), f (x) and I(x) for given S0 = S are equal to Sj1 (x), Sf1 (x) and SI1 (x), where j1 (x), f1 (x) and I1 (x) are the solutions of S0 = 1. On the outside of the halo, r > R, no photons propagate in the direction µ < 0. The boundary condition at r = R of (6.1) should be I(η, R, x, µ) = 0, µ < 0. (6.17) R −1 For (6.11), we have 0 µI(η, R, x, µ)dµ = 0 [108], and the boundary condition is then j(η, R, x) = 2f (η, R, x). (6.18) If the source becomes to emit photon at t = 0, the initial condition should be I(0, r, x, µ) = 0, (6.19) for (6.1), and j(0, r, x) = f (0, r, x) = 0, (6.20) for (6.11). 6.1.4 Re-distribution function In this section, we proceed to the re-distribution function R(x, µ, x0 , µ0 ; a) in (6.1). We consider isotropic scattering and the case of a = 0. The re-distribution function · 2 ¸ 0 0 1 x − 2xx0 cos α + x02 R(x, n, x , n ) = exp − π sin α sin2 α gives the probability that a photon with frequency x0 and direction n0 within an element of solid angle dω 0 is absorbed and re-emitted with frequency x and direction 114 n in dω [56], where α is the angle between n and n0 . Choose x-axis such that n0 lies in the xz-plane, then n0 = (sin θ0 , 0, cos θ0 ) and n = (sin θ cos φ, sin θ sin φ, cos θ), where φ is the azimuthal angle, µ = cos θ and µ0 = cos θ0 . With the above notation, 1 dω = 4π dφdµ and cos α = n0 · n = sin θ sin θ0 cos φ + cos θ cos θ0 . Integrating over φ, we obtain Z 2π · 2 ¸ 0 0 1 2 − 2xx0 H + x02 R(x, µ, x , µ ) = √ exp − dφ, 0 2π 2 1 − H 2 1 − H2 p p where H = 1 − µ2 1 − µ02 cos φ + µµ0 . If we consider a 6= 0, and follow the same line above, we obtain Z Z " µ ¶2 #−1 µ ¶ 2π ∞ 0 0 a −u2 2 x + x0 (x − x0 )2 R(x, µ, x , µ ; a) = e a + − αu exp − dudφ, 0 4π 3 β −∞ 2 4β 2 p p q q 1+H where H = 1− µ2 02 0 1 − µ cos φ + µµ , α = 2 , and β = 1−H 2 . We can verify numerically that the angular averaged re-distribution function is exactly the same as the one obtained by Hummer [56], i.e. Z 1 Z 1 1 0 01 0 R(x, µ, x , µ ; a)dµ = R(x, µ, x0 , µ0 ; a)dµ = R(x, x0 ; a). 2 −1 2 −1 115 6.1.5 Integral of the phase function Equation (6.6) can be rewritten as Z 2π d 0 1 1 − g2 R (µ, µ ) = dφ0 3 (6.21) 4π 0 |I − gI0 | 2 where I and I0 are unit vector on the direction of polar angle θ and θ0 , and azimuth angle φ and φ0 , respectively. That is I · I = I0 · I0 = 1 and I · I0 = cos γ = cos θ cos θ0 + sin θ sin θ0 cos(φ − φ0 ), and µ = cos θ, µ0 = cos θ. We have d 1 1 − g2 1 0 1/2 = 0 3/2 − , dg |I − gI | 2g|I − gI | 2g|I − gI0 |1/2 and therefore, 1 − g2 d 1 1 0 3/2 = 2g 0 1/2 + . (6.22) |I − gI | dg |I − gI | |I − gI0 |1/2 The expansion with Legendre functions Pl (cos γ) gives X ∞ 1 = g l Pl (cos γ), (6.23) |I − gI0 |1/2 l=0 and then X ∞ X ∞ 1 − g2 l = 2lg P l (cos γ) + g l Pl (cos γ). (6.24) |I − gI0 |3/2 l=1 l=0 Since cos γ = cos θ cos θ0 + sin θ sin θ0 cos(φ − φ0 ), we have the following identity for the Legendre function Pl (cos γ) as m=l X (l − m)! m Pl (cos γ) = Pl (cos θ)Pl (cos θ0 ) + 2 Pl (cos θ)Plm (cos θ0 ) cos[m(φ − φ0 )]. m=1 (l + m)! (6.25) 116 The integral of φ0 in (6.21) kills the second term of (6.25), we have "∞ ∞ # 1 X X Rd (µ, µ0 ) = 2π 2lg l Pl (cos θ)Pl (cos θ0 ) + g l Pl (cos θ)Pl (cos θ0 ) (6.26) 4π " ∞ l=1 ∞ l=0 # 1 X l X = 2lg Pl (µ)Pl (µ0 ) + g l Pl (µ)Pl (µ0 ) . 2 l=1 l=0 R1 2 Using the orthogonal relation −1 Pl (µ)Pl0 (µ)dµ = δ 0, 2l+1 l,l we have Z 1 Z 1 1 R0 (g) = dµ dµ0 Rd (µ, µ0 ) = 1, (6.27) 2 −1 −1 for which only the term l = 0 in (6.26) has contribution. Similarly, Z 1 Z 1 Z 1 Z 1 1 0 d 1 0 R1 (g) = dµ dµ µR (µ, µ ) = dµ dµ0 µ0 Rd (µ, µ0 ) = 0, (6.28) 2 −1 −1 2 −1 −1 Z 1 Z 1 1 g R2 (g) = dµ dµ0 µµ0 Rd (µ, µ0 ) = . (6.29) 2 −1 −1 3 These results are used in deriving (6.11) and (6.12). 6.2 Numerical algorithm for equation (6.1) We solve (6.1) with initial and boundary conditions (6.13), (6.17) and (6.19). For simplicity, we ignore the effect of dust (i.e. κ(x) = 0). Our computational domain is (r, x, µ) ∈ [0, rmax ] × [xleft , xright ] × [−1, 1], where rmax , xleft and xright are chosen such that the solution vanishes to zero outside the boundaries. We choose mesh sizes with grid refinement tests to ensure proper numerical resolution. In the following, we describe numerical techniques involved in our algorithm, including approximations to spatial derivatives, numerical integration, numerical boundary condition and time evolution. 117 6.2.1 Conservation law To use the WENO algorithm, we first rewrite (6.1) into the form of a conservation law. Noticing the boundary condition (6.13), we define I 0 = r2 I, so(6.1) becomes ∂I 0 ∂I 0 1 ∂(1 − µ2 )I 0 ∂I 0 +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, µ, x0 , µ0 ; a)I 0 (η, r, x0 , µ0 )dx0 dµ0 /2 + r2 S. 0 (6.30) For simplicity, we drop the prime, and use I(η, r, x, µ) for I 0 (η, r, x, µ) below. 6.2.2 The WENO algorithm: approximations to the spatial derivatives The spatial derivative terms in (6.30) are approximated by a fifth-order finite differ- ence WENO scheme. ∂I We first give the WENO reconstruction procedure for approximating ∂x , ∂I(η n , ri , xj , µk ) 1 ˆ ˆ j−1/2 ) ≈ (hj+1/2 − h ∂x ∆x ˆ j+1/2 is obtained by the with fixed η = η n , r = ri and µ = µk . The numerical flux h fifth-order WENO approximation in an upwind fashion, because the wind direction is fixed (negative). Denote hj = I(η n , ri , xj , µk ), j = −2, −1, · · · , Nx + 3 with fixed n, i and k. The numerical flux from the WENO procedure is obtained by ˆ j+1/2 = ω1 h h ˆ (3) , ˆ (2) + ω3 h ˆ (1) + ω2 h (6.31) j+1/2 j+1/2 j+1/2 118 ˆ (m) are the three third-order fluxes on three different stencils given by where h j+1/2 ˆ (1) = − 1 hj−1 + 5 hj + 1 hj+1 , h j+1/2 6 6 3 ˆ (2) 1 5 1 h j+1/2 = hj + hj+1 − hj+2 , 3 6 6 ˆ (3) 11 7 1 h j+1/2 = hj+1 − hj+2 + hj+3 , 6 6 3 and the nonlinear weights ωm are given by ω ˘m γl ωm = P3 , ω ˘l = , l=1 ω ˘l (² + βl )2 where ² is a parameter to avoid the denominator to become zero and is taken as ² = 10−8 . The linear weights γl are given by 3 3 1 γ1 = , γ2 = , γ3 = , 10 5 10 and the smoothness indicators βl are given by 13 1 β1 = (hj−1 − 2hj + hj+1 )2 + (hj−1 − 4hj + 3hj+1 )2 , 12 4 13 1 β2 = (hj − 2hj+1 + hj+2 ) + (hj − hj+2 )2 , 2 12 4 13 1 β3 = (hj+1 − 2hj+2 + hj+3 ) + (3hj+1 − 4hj+2 + hj+3 )2 . 2 12 4 ∂(1−µ2 )I Similarly, we give the WENO procedure in approximating ∂µ , ∂(1 − µ2j )I(η n , ri , xk , µj ) 1 ˆ ≈ ˆ j−1/2 ) (hj+1/2 − h ∂µ ∆µ ˆ j+1/2 is also obtained with fixed η = η n , r = ri and x = xk . The numerical flux h by the fifth-order WENO approximation in an upwind fashion, however the wind 119 ∂I direction here is positive, opposite from that of ∂x . Denote hj = (1 − µ2j )I(η n , ri , xk , µj ), j = −3, −2, · · · , Nµ + 2 with fixed n, i and k. The numerical flux from the WENO procedure is obtained by ˆ j+1/2 = ω1 h h ˆ (1) + ω2 h ˆ (2) + ω3 h ˆ (3) , (6.32) j+1/2 j+1/2 j+1/2 ˆ (m) are the three third-order fluxes on three different stencils given by where h j+1/2 ˆ (1) = − 1 hj+2 + 5 hj+1 + 1 hj , h j+1/2 6 6 3 1 5 1 ˆ (2) = hj+1 + hj − hj−1 , h j+1/2 3 6 6 ˆ (3) = 11 7 1 h j+1/2 hj − hj−1 + hj−2 , 6 6 3 and the nonlinear weights ωm are given as ω ˘m γl ωm = P3 , ω ˘l = , l=1 ω ˘l (² + βl )2 where ² is taken as ² = 10−8 . The linear weights γl are also given by 3 3 1 γ1 = , γ2 = , γ3 = , 10 5 10 and the smoothness indicators βl are given by 13 1 β1 = (hj+2 − 2hj+1 + hj )2 + (hj+2 − 4hj+1 + 3hj )2 , 12 4 13 1 β2 = (hj+1 − 2hj + hj−1 )2 + (hj+1 − hj−1 )2 , 12 4 13 1 β3 = (hj − 2hj−1 + hj−2 )2 + (3hj − 4hj−1 + hj−2 )2 . 12 4 120 In the end, we approximate the r-derivative in (6.30), following the reconstruction procedures mentioned above. However, we need to check the wind direction at the r-boundary of each cell. When µ > 0, the wind direction is positive, and we use (6.32) to approximate the numerical flux, while when µ < 0, we use equation (6.31). 6.2.3 High order numerical integration The integration of the resonance scattering term is calculated by a fifth order quadra- ture [100] Z µright Nµ X f (µ)dµ = ∆µ ωk f (µk ) + O(∆µ5 ), µlef t k=1 where µk = µlef t + (k − 12 )dµ and the weights are defined as, 6463 1457 741 5537 ω1 = , ω2 = , ω3 = , ω4 = , 5760 1920 640 5760 5537 741 1457 6463 ωNµ −3 = , ωNµ −2 = , ωNµ −1 = , ω Nµ = , 5760 640 1920 5760 and ωk = 1 otherwise. 6.2.4 Numerical boundary condition Following Carrillo et al. [20], at µ = −1 and µ = 1, we take the boundary conditions as, for µ > 0, I(η, r, x, −1 − µ) = I(η, r, x, −1 + µ), I(η, r, x, 1 + µ) = I(η, r, x, 1 − µ), motivated by the physical meaning of µ as the cosine of the angle to the z−axis. ˆ1 = h We also explicitly impose h ˆ Nµ + 1 = 0 for the first and last numerical fluxes in 2 2 order to enforce conservation of mass. 121 6.2.5 Time Evolution To evolve in time, we use the third-order TVD Runge-Kutta time discretization [102]. For system of ODEs ut = L(u), the third order Runge-Kutta method is u(1) = un + ∆τ L(un , τ n ), 3 1 u(2) = un + (u(1) + ∆τ L(u(1) , τ n + ∆τ )), 4 4 1 2 1 un+1 = un + (u(2) + ∆τ L(u(2) , τ n + ∆τ )). 3 3 2 6.3 Numerical algorithm for equations (6.11) and (6.12) We will solve (6.11) and (6.12) with boundary and initial conditions (6.14), (6.18) and (6.20) by using the WENO solver. To solve (6.11) and (6.12) as a system, our computational domain is (r, x) ∈ [0, rmax ] × [xleft , xright ]. As mentioned in the Section 6.2, rmax , xleft and xright are cho- sen such that the solution vanishes outside the boundaries. In the following, we describe numerical techniques involved in our algorithm, including the characteris- tic decomposition, and numerical boundary condition. All other algorithms can be found in Section 6.2. 6.3.1 Characteristic decomposition We consider the WENO reconstruction procedure for approximating the r-derivatives only. We need to perform the WENO procedure based on a characteristic decompo- 122 sition. To accomplish this, we write the left-hand side of (6.11) and (6.12) as ut + Aur , where u = (j, f )T and   01 A=  1 3 0 is a constant matrix. To perform the characteristic decomposition, we first compute the eigenvalues, the right eigenvectors and the left eigenvectors of A and denote them by Λ, R and R−1 . We then project u to the local characteristic fields v with v = R−1 u. Now ut + Aur of the original system is decoupled as two independent equations as vt + Λvr . We approximate the derivative vr component by component, each with the correct upwind direction, with the WENO reconstruction procedure similar to the procedure described in Section 6.2. In the end, we transform vr back to the physical space by ur = Rvr . We refer the readers to [29] for more implementation details. 6.3.2 Numerical Boundary Condition To implement the boundary condition (6.18), we also need to perform a characteristic decomposition as discussed above. Using the same notation as before, we project u to the local characteristic fields v with v = R−1 u. Denote v = (v1 , v2 )T , now ut + Aur of the original system is decoupled to two independent scalar operators given by ∂v1 ∂v1 ∂v2 ∂v2 + λ1 ; + λ2 ∂t ∂r ∂t ∂r √ √ 3 3 where λ1 = 3 and λ2 = − 3 . The characteristic line starting from the boundary r = rmax for the first equation is pointing outside the computational domain while the one for the second equation is pointing inside. For well-posedness of our system, 123 we need to impose the boundary condition there as v2 = αv1 + β with constants α and β. We can calculate the values of α and β based on equation (6.18) and the left and right eigenvectors of A. For example, if we take √ √  3 3 R= 2 2 , 1 2 − 21 √ we recover that α = 7 − 4 3 and β = 0. We use extrapolation to obtain the value of v1 and then compute the value v2 . In the end, we transfer v back to the physical space by u = Rv. Chapter 7 Effect of dust on Lyα photon transfer in optically thick halo We will investigate, in this chapter, the effects of the dust on the Lyα photons transfer in an optically thick medium. Dust can be produced at epochs of low and moderate redshifts, and even at redshift as high as 6 [106]. Absorption and scattering of dust have been used to explain the observations on Lyα emission and absorption [59], such as the escaping fraction of Lyα photons [52, 53, 11]; the redshift-dependence of the ratio between Lyα emitters and Lyman Break galaxies [109]; and the “evolution” of the double-peaked profile [71]. However, it is still unclear whether the time scale of a photon escaping from optically thick halo will be increasing (or decreasing) when the halo is dusty. It is also unclear whether the effects of dust absorption can be estimated by the random walk picture [50]. As for the dust effect on the double-peaked profile, the current results given by different studies seem to be contradictory: some claims that the dust absorption leads to the narrowing of the double-peaked profile [71], while others conclude that the width between the two peaks apparently should be increasing due to the dust absorption [110]. We will focus on these basic problems, and examine them with the solution of the integro-differential equation of radiative transfer. 124 125 7.1 Basic theory We study the radiative transfer (6.11) and (6.12) with boundary and initial condition (6.14), (6.18) and (6.20). 7.1.1 Dust models We consider three models of the dust as follows: I. pure scattering, A = 1, g = 0.73: dust causes only anisotropic scattering, but no absorption; II. scattering and absorption, A = 0.32, g = 0.73: dust causes both absorption and anisotropic scattering. III. pure absorption, A = 0: dust causes only absorption, but no scattering; Models I and III do not occur in reality. They are, however, helpful to reveal the effects of pure scattering and absorption on the radiative transfer. Since κ(x) is on the order of 10−8 , its effect will be significant only for halos with optical depth τ0 ≥ 106 , and ignorable for τ0 ≤ 105 . To illustrate the dust effect, we use halos of R = τ0 ≤ 104 , and take larger κ to be ' 10−4 − 10−2 . We consider below only the case of grey dust, i.e. κ is independent of frequency x. This certainly is not realistic dust. Yet, the frequency range given in the solution below mostly are in the range |x| < 4. Therefore, the approximation of grey dust would be proper if the cross section of dust is not strongly frequency dependent in the range |x| < 4. 7.1.2 Numerical example: Wouthuysen-Field thermalization As the first example of numerical solutions, we show the Wouthuysen-Field (W-F) effect, which requires that the distribution of Lyα photons in the frequency space should be thermalized near the resonant frequency ν0 . The W-F effect illustrates the 126 difference between the analytical solutions of the Fokker-Planck approximation and that of equations (6.11) and (6.12). The former can not show the local thermaliza- tion [76], while the latter can [88]. All problems related to the W-F local thermal equilibrium should be studied with the integro-differential equation (6.1). 100 100 κ=0 100 κ=0 κ=10 -4 κ=10 -4 j(η,r,x) j(η,r,x) j(ηr,x) 10-1 10-1 10-1 κ=10 -3 κ=10 -3 κ=10 -2 κ=10 -2 -2 -2 -2 10 10 10 -5 0 5 -5 0 5 -5 0 5 x x x Figure 7.1: The mean intensity j(η, r, x) at η = 500 and r = 100 for dust models I (left panel), √ II −x (middle panel) and III (right panel). The source is S0 = 1 and 2 φs (x) = (1/ π)e . The parameter a = 10−3 . In each panel, κ is taken to be 0, 10−4 , 10−3 and 10−2 . Figure 7.1 presents a solution of the mean intensity j(η, r, x) at time radial η = 500 coordinate r = 102 for halo with size R À r = 102 . The three panels correspond to dust models I (left panel), II (middle panel) and III (right panel). The source is √ 2 taken to have a Gaussian profile φs (x) = (1/ π)e−x and unit intensity S0 = 1. The solutions of Figure 7.1 actually are independent of R, if R À 102 . The intensity of j is decreasing from left to right in Figure 7.1, because the absorption is increasing with the models from I to III. A remarkable feature shown in Figure 7.1 is that all j(η, r, x) have a flat plateau in the range |x| ≤ 2. This gives the frequency range of the W-F local thermalization [88, 89]. The range of the flat plateau |x| ≤ 2 is almost dust-independent, either for model I or for models II and III. This is expected, as neither the absorption nor scattering given by the κ term of (6.1) changes the frequency distribution of photons. The redistribution function (6.6) also does not change the frequency distribution of photons. This point can also be seen from (6.11) and (6.12), in which the κ terms 127 are frequency-independent. The evolution of the frequency distribution of photons is due only to the resonant scattering. Since thermalization will erase all frequency features within the range |x| ≤ 2, the double-peaked structure does not retain information of the photon frequency distribution within |x| < 2 at the source. That is, the results in Figure 7.1 will hold for any source S0 φs (x) with arbitrary φs (x) which is non-zero within |x| < 2 [88, 89]. This property can also be used as a test of the simulation code. It requires that simulation result in a flat plateau, regardless of the whether source is monochromatic or with a finite width around ν0 . 7.2 Dust effects on photon escape 7.2.1 Model I: scattering of dust To study the effects of dust scattering on the Lyα photon escape, we show in Figure 7.2 the flux f (η, r, x) of Lyα photons emerging from halos at the boundary r = R = 102 for Model I. The three panels of Figure 7.2 correspond to κ = 10−4 , 10−3 , and 10−2 from left to right, respectively. The source starts to emit photons at η = 0 with √ 2 a stable luminosity S0 = 1, and with a Gaussian profile φs (x) = (1/ π)e−x . Figure 7.2 clearly shows that the time-evolution of f (η, r, x) is κ-independent. Although the cross section of dust scattering increases about 100 times from κ = 10−4 to κ = 10−2 , the curves of the left and right panels in Figure 7.2 are almost identical. According to the scenario of “single longest excursion”, photon escape is not a process of Brownian random walk in the spatial space, but a transfer in the frequency space [80, 9, 1, 2, 51, 12]. A photon will escape, once its frequency is transferred from |x| < 2 to |x| > 2, on which the medium is transparent. On the other hand, dust scattering given by the redistribution function equation (6.6) does not change photon frequency. Dust scattering has no effect on the transfer in the frequency space. 128 0.6 0.6 0.6 0.4 0.4 0.4 η=2000,3000 η=2000,3000 η=2000,3000 0.2 0.2 0.2 f(η,r,x) f(η,r,x) f(η,r,x) η=500 η=500 η=1000 η=500 η=1000 η=1000 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 x x x Figure 7.2: Flux f (η, r, x) of Lyα photons emergent from halos at the boundary R = 102 , and for the dust model I A = 1, g = 0.73. The parameter κ is taken to√be 10−4 2 (left), 10−3 (middle) and 10−2 (right). The source is S0 = 1 and φs (x) = (1/ π)e−x . The parameter a = 10−3 . Moreover, photons with frequency |x| < 2 are quickly thermalized after a few hundred resonant scattering. In the local thermal equilibrium state, the angular distribution of photons is isotropic. Thus, even if the dust scattering is anisotropic g 6= 0 with respect to the direction of the incident particle, the angular distribution will keep isotropic undergoing a g 6= 0 scattering. Hence, dust scattering also has no effect on the angular distribution. 7.2.2 Model III: absorption of dust Similar to Figure 7.2, we present in Figure 7.3 the flux of Model III, i.e. dust causes only absorption without scattering. All other parameters of Figure 7.3 are the same as in Figure 7.2. In the left panel of Figure 7.3, the curves at the time η = 2000 and 3000 are the same. It means the flux f (η, R, x) at the boundary R is already stable, or saturated at the time η ≥ 2000. The small difference between the curves of η = 1000 and η ≥ 2000 of the left panel indicates that the flux is still not yet completely saturated at the time η = 1000. However, comparing the middle and right panels of Figure 7.3, we see that for κ = 10−3 , the flux has already saturated at η = 1600, while it has saturated at η = 800 for κ = 10−2 . That is, the stronger the 129 dust absorption, the shorter the saturation time scale. The time scales of escape or saturation do not increase by dust absorption, and even decrease with respect to the medium without dust. Stronger absorption leads to shorter time scale of saturation. 0.4 0.04 0.4 η=2000,3000 η=1600,2500 η=800,1500 0.02 0.2 0.2 f(η,r,x) f(η,r,x) f(η,r,x) η=1000 η=800 η=500 η=500 η=200 η=300 -2 0 2 -2 0 2 -2 0 2 x x x Figure 7.3: Flux f (η, r, x) of Lyα photons emergent from halos at the boundary r = R = 102 . The parameters of the dust are A = 0 and κ = 10−4 (left), 10−3 (middle) and 10−2 (right). Other parameters are the same as in Figure 7.2. Obviously, dust absorption does not help in producing photons for the “single longest excursion”. Therefore, dust absorption cannot make the time scale of pro- ducing photons for “single longest excursion” smaller. However, dust absorptions are effective in reducing the number of photons trapped in the state of local thermalized equilibrium |x| < 2 (see also Section 7.3.2). This indicates that the higher the value of κ, shorter the time scale of saturation. 7.2.3 Effective absorption optical depth Since Lyα photons underwent a large number of resonant scattering before escap- ing from the halo with optical depth τ0 À 1, it is generally believed that a small absorption of dust will lead to a significant decrease of the flux. However, it is still unclear what the exact relationship between the dust absorption and the resonant scattering is. This problem should be measured by the effective optical depth of dust absorption of Lyα photons in R = τ0 À 1 halos. To calculate the effective optical depth, we first give the total flux of Lyα photons 130 R emerging from a halo of radius R, which is defined as F (η) = f (η, R, x)dx. Figure 7.4 plots F (η) as a function of time η for halo with sizes R = τ0 = 102 and 104 . The curves typically are growing and then saturating. The three panels correspond to the dust models I, II and III from left to right. The upper panels are of R = 102 , and lower panels for R = 104 . In each panel of R = 102 , we have three curves corresponding to κ = 10−4 , 10−3 and 10−2 , respectively. In cases of R = 104 , we take κ = 10−4 and 10−3 . κ=10 κ=10 -4 -4 100 100 100 κ=10 -3 κ=10 -3 F(η) F(η) F(η) 10-1 10-1 10-1 κ=10 -2 κ=10 -2 10-2 10-2 10-2 0 1000 2000 3000 0 1000 2000 3000 0 1000 2000 3000 η η η 0 10 100 100 κ=10 -4 10-2 10-2 κ=10 -4 10 -2 κ=10 -4 κ=10 -3 10-4 10-4 10-4 F(η) F(η) F(η) 10-6 10 -6 10-6 κ=10 -3 κ=10 -3 10-8 10-8 10-8 10-10 10-10 10-10 0 200000 0 200000 0 200000 η η η Figure 7.4: The time evolution of the total flux F (η) at the boundary of halos with R = τ0 = 102 (upper √ panels), and R = τ0 = 104 (lower Panels). The source of S0 = 1 2 and φs (x) = (1/ π)e−x starts to emit photons at time η = 0. The parameters of dust are (A = 1, g = 0.73) (left); (A = 0.32, g = 0.73) (middle) and A = 0 (right). In each panel of R = 102 , κ is taken to be 10−4 , 10−3 and 10−2 . In the cases of R = 104 , κ is taken to be 10−4 , 10−3 . The left panel of Figure 7.4 shows that the three curves of κ = 10−4 , 10−3 and 10−2 are almost the same. This is consistent with Figure 7.2 that for Model I, the time-evolution of f are κ-independent for the pure scattering dust. For the pure 131 absorption dust (the right panel of Figure 7.4), the saturated flux is smaller for larger κ. We can also see from Figure 7.4 that the time scale of approaching saturation is smaller for larger κ. The result of model II is in between that for models I and III. With the saturated flux of Figure 7.4, one can define the effective absorption optical depth by τeffect ≡ −(1/κ) ln fS . The results are shown in Table 7.1, in which τa is the dust absorption depth. It is interesting to see that the effective absorption optical depth is always equal to a few times of the optical depth of resonant scattering τ0 , regardless of whether τa is less than 1. Namely, the effective absorption depth τeffect of dust is roughly proportional to τ0 . Table 7.1: Effective absorption optical depth τeffect Model II Model III R = τ0 κ τa fS τeffect τa fS τeffect 102 10−4 0.0068 0.978 2.2 × 102 0.01 0.963 3.8 × 102 102 10−3 0.068 0.760 2.7 × 102 0.10 0.670 4.0 × 102 102 10−2 0.68 0.116 2.2 × 102 1.00 0.057 2.9 × 102 104 10−4 0.68 6.28 × 10−2 2.8 × 104 1.00 3.02 × 10−2 3.5 × 104 104 10−3 6.8 4.07 × 10−7 1.5 × 104 10.0 2.87 × 10−9 1.97 × 104 According to the random walk scenario, if a medium has optical depths of ab- sorption τa and scattering τs , the effective absorption optical depth should be equal p to τeffect = τa (τa + τs ) [95]. However, the results of the last line of Table 7.1 show that the random walk scenario does not work for the dust effect on resonant photon transfer. This result is consistent with Figures 7.2 and 7.3. When the optical depth of dust is lower than the optical depth of resonant scattering τ0 , the time scale of photon escaping is not affected by the dust, but is proportional to τ0 , and therefore, the absorption is also proportional to τ0 . 132 7.2.4 Escape coefficient With the total flux, we can define the escaping coefficient of Lyα photon as fesc (η, τ0 ) ≡ F (η)/F0 , where F0 is the flux of the center source. Figure 7.5 shows fesc (η, τ0 ) at three times η = 5 × 103 , 104 and 3.2×104 for Model II and κ = 10−3 . At η = 5 × 103 , the flux of halos with τ0 ≤ 103 is saturated. At η = 104 , halos with τ0 ≤ 3 × 103 are saturated, and all halos of τ0 ≤ 104 are saturated at η = 3.2 × 104 . 100 10-1 -2 10 fesc(η,τ0) 10-3 10-4 10-5 10-6 2 3 4 10 10 10 τ0 Figure 7.5: Escaping coefficient fesc (η) as a function of the optical depth τ0 of halo at time η = 5 × 103 , 104 , and 3.2×104 from bottom to up. Dust is modeled by II, A = 0.32, g = 0.73, and κ = 10−3 . 7.3 Dust effects on double-peaked profile 7.3.1 Dust and the frequency of double peaks A remarkable feature of Lyα photon emerging from an optically thick medium is the double-peaked profile. Figures 7.1, 7.2 and 7.3 have shown that the double peak frequencies x+ = |x− | are almost independent of either the scattering or the absorption of dust. In this section, we consider halos of size R or τ0 larger than 102 . Figure 7.6 presents the double peak frequency |x± | as a function of aτ0 , where the 133 parameter a is taken to be 10−2 (left) and 5 × 10−3 (right). Comparing the curves with dust and without dust in Figure 7.6 we conclude that the dust effect on |x± | is very small woth aτ0 = aR = 102 . 5 4.5 4 4 3.5 -3 A= 0 0 3 -3 A= 0 κ= 0 3 κ= 1 κ= 0 0 κ= 1 2.5 X+ X+ 2 2 1.5 1 0 1 10 10 1 10 2 20 40 60 aτ0 aτ0 Figure 7.6: The two-peak frequencies x+ = |x− | as a function of aτ0 . The parameter a is taken to be 10−2 (left) and 5 × 10−3 (right). Dust model III (pure absorption) is used, and κ is taken to be 10−3 . The dashed straight line gives log x± -log aτ with slope 1/3, which is to show the (aτ )1/3 -law of x± . In the range aτ0 < 20, the |x± |-τ0 relation is almost flat with |x± | ' 2. It is because the double-peaked profile is given by the frequency range of the locally thermal equilibrium. The positions of the two peaks, x+ and x− , essentially are at the maximum and minimum frequencies of the local thermalization. The frequency range of the local thermal equilibrium state is determined mainly by the Doppler broadening, and weakly dependent on τ0 . Thus, we always have x± ' ±2. When the optical depth is larger, aτ0 ∼ 102 , more and more photons of the flux are attributed to the resonant scattering by the Lorentzian wing of the Voigt profile. In this phase, |x± | will increase with τ0 . Figure 7.6 also shows a line x± = ±(aτ0 )1/3 , which is given by the analytical solution of the Fokker-Planck approximation, in which the Doppler broadening core in the Voigt profile being ignored [51, 76, 37]. The numerical solutions of (6.1) or (6.11) and (6.12) deviate from the (aτ0 )1/3 -law at all parameter range of Figure 134 7.6. The deviation at aτ0 < 20 is caused by the Doppler broadening core in the Voigt profile is ignored in the Fokker-Planck approximation, so no locally thermal equilibrium can be reached. Therefore, in the range aτ0 < 20, |x± | of the WENO solution is larger than the (aτ0 )1/3 -law. In the range of aτ0 > 20, the Fokker-Planck approximation yields a faster diffusion of photons in the frequency space. This point can be seen in the comparison between a Fokker-Planck solution with Field’s analytical solution (Figure 1 in [94]). In this range, the numerical results of |x± | is less than the (aτ0 )1/3 -law. 7.3.2 No narrowing and no widening The dust effect has been used to explain the narrowing of the width between the two peaks [71]. conversely, it is also used to explain the widening of the width between the two peaks [110]. However, Figures 7.1, 7.2, 7.3 and 7.6 show that the width between the two peaks of the profile is very weakly dependent on dust scattering and absorption. This result supports, at least in the parameter range considered in Figures 7.1, 7.2, 7.3, neither the narrowing nor the widening of the two peaks. If dust absorption can cause narrowing, the absorption should be weaker at |x| ∼ 0, and stronger at |x| ≥ 2. Similarly, if dust absorption can cause widening, the absorption should be weaker at |x| ∼ 2, and stronger at |x| ∼ 0. To test these assumptions, Figure 7.7 plots ln[f (η, r, x, κ = 0)/f (η, r, x, κ)] as a function of x. It measures the x(frequency)-dependence of the flux ratio with and without dust absorption. We take large η, and then the fluxes in Figure 7.7 are saturated. Figure 7.7 shows that the absorption in the range |x| < 2 is much stronger than for |x| > 2, and therefore, the assumption of the narrowing is ruled out. Figure 7.7 shows also that the curves of ln[f (η, r, x, κ = 0)/f (η, r, x, κ = 10−3 )] are almost flat in the range |x| < 2. Therefore, the assumption of widening of the two peaks can also be ruled out. Since the cross sections of dust absorption and scattering are assumed to be 135 2.1 ln[f(η,r,x,κ=0)/f(η,r,x,κ=10 -3)] ln[f(η,r,x,κ=0)/f(η,r,x,κ=10 -2)] 0.22 2 0.21 1.9 0.2 1.8 0.19 1.7 0.18 1.6 -4 -2 0 2 4 -4 -2 0 2 4 x x ln[f(η,r,x,κ=0)/f(η,r,x,κ=10 -3)] ln[f(η,r,x,κ=0)/f(η,r,x,κ=10 -2)] 0.32 2.8 0.31 0.3 2.6 0.29 0.28 2.4 0.27 0.26 2.2 -4 -2 0 2 4 -4 -2 0 2 4 x x Figure 7.7: ln[f (η, r, x, κ = 0)/f (η, r, x, κ)] as function of x for model II (up), and III (bottom), and κ = 10−3 (left) and 10−2 (right). Other parameters are the same as in Figure 7.2. frequency-independent. (6.11) and (6.12) do not contain any frequency scales other than that from resonant scattering. However, either narrowing or widening would require to have frequency scales different from that of resonant scattering. This is not possible if the dust is gray. 7.3.3 Profile of absorption spectrum If the radiation from the sources has a continuous spectrum, the effect of a neutral hydrogen halos is to produce an absorption line at ν = ν0 . The profile of the absorption line can also be found by solving (6.11) and (6.12), but replacing the boundary equation (6.14) by f (η, 0, x) = S0 . (7.1) 136 That is, we assume that the original spectrum is flat in the frequency space. The spectrum of the flux emerging from the halo of R = 102 and 104 with central source for (7.1) for dust models I, II and III are shown in Figure 7.8. All curves are for large η, i.e. they are saturated. κ=0 κ=0 0 0 0 10 10 10 f(η,r,x) f(η,r,x) f(η,r,x) κ=10 κ=10 -4 -4 -1 -1 -1 κ=10 κ=10 10 10 -3 10 -3 κ=10 -2 κ=10 -2 10-2 10-2 10-2 -5 0 5 -5 0 5 -5 0 5 x x x 100 100 κ=0 100 κ=0 -2 -2 -2 10 10 10 κ=10 -4 f(η,r,x) f(η,r,x) f(η,r,x) 10 -4 10 -4 10-4 κ=10 -4 10-6 10-6 κ=10 -3 10-6 κ=10 -8 -8 -8 -3 10 10 10 10-10 10-10 10-10 -5 0 5 -5 0 5 -5 0 5 x x x Figure 7.8: The spectrum of the flux emergent from halo of R = 102 (upper panels) and 104 (lower panels) with central source of equation (7.1) for the dust model I (left), II (middle) and III (right). Other parameters are the same as in Figure 7.2. The optical depths at the frequency |x| > 4 are small, and therefore, the Ed- dington approximation might no longer be valid. However, those photons do not strongly involve the resonant scattering, and hence they do not significantly affect the solution around x = 0. The solutions of Figure 7.8 is still useful to study the profiles of f around x = 0. The flux profile of Figure 7.8 are absorption lines with the width given by the double peaks similar to the double peaked structure of the emission line. The flux at the double peaks is even higher than for the flat wing. It is because more photons 137 are stored in the frequency range |x| < 2. According to the redistribution function equation (6.7), the probability of transferring a x0 photon to a |x| < |x0 | photon is larger than that from |x0 | to |x| > |x0 |. Therefore, if the original spectrum is flat, the net effect of resonant scattering is to bring photons with frequency |x| > 2 to |x| < 2. Photons stored |x| < 2 are thermalized, and therefore, in the range |x| < 2, the profile will be the same as the emission line, and the double peaks can be higher than the wing. It makes the shoulder at |x| ∼ 2. As expected, for model I (left panels of Figure 7.8), the double profile is com- pletely κ-independent. Dusty scattering does not change the flux and its profile. For models II and III, the higher the κ, the lower the flux of the wing, because the dust absorption is assumed to be frequency-independent. The positions of the double peaks, in the absorption spectrum are also κ-independent. This once again shows that dust absorption and scattering causes neither narrowing nor widening of the double-peaked profile. However, for higher κ the flux of the peaks is lower. When the absorption is very strong, the double-peaked structure might disappear, but will not be narrowed or widened. 7.4 Discussions and conclusions The study of dust effects on radiative transfer has had a long history related to extinction. However, dust effects on radiative transfer of resonant photons actually have not been carefully investigated. Existing works are mostly based on the solu- tions of the Fokker-Planck approximation, or Monte Carlo simulation. These results are important. We revisited these problems with the WENO solver of the integro- differential equation of the resonant radiative transfer, and have found some features which have not been addressed in previous works. These features are summarized as follows. First, the random walk picture in the physical space will no longer be available 138 for estimating the effective optical depth of dust absorption. For a medium with the optical depth of absorption and resonant scattering to be τa À 1, τ (ν0 ) À 1 and τs (ν0 ) À τa , the effective absorption optical depth is found to be almost independent of τa , and to be equal to about a few times of τs (ν0 ). Second, dust absorption will, of course, yield the decrease of the flux of Lyα pho- tons emergent from optical thick medium. However, if the absorption cross-section of dust is frequency independent, the double-peaked structure of the frequency pro- file is basically dust-independent. The double-peaked structure does not narrow or widen by the absorption and scattering of dust. Third, the time scales of Lyα photon transfer basically are independent of dust scattering and absorption. Since these time scales are mainly determined by the ki- netics in the frequency space. However, dust does not affect the behavior of the trans- fer in the frequency space if the cross section of the dust is wavelength-independent. The local thermal equilibrium makes the anisotropic scattering ineffective on the angular distribution of photons. Dust absorption and scattering do not lead to the increase or decrease of the time of storing Lyα photons in the halos. The differences between the time-independent solutions of the Fokker-Planck approximation, or Monte Carlo simulation and the WENO solution of (6.1) is mainly related to the W-F effect. Therefore, all above-mentioned features can already be clearly seen with halos of τ0 ∼ 102 , in which the W-F local thermal equilibrium has been well established. In this context, most calculation in this chapter is on holes with τ0 < 105 . This range of τ0 certainly is unable to describe halos with column number density of HI larger than 1017 cm−2 (e.g. [90]). Nevertheless, the result of τ0 < 105 would already be useful for studying the 21 cm region around high-redshift sources, of which the optical depth typically is [74, 89]. µ ¶−1/2 µ ¶3 µ ¶µ ¶ 5 T 1+z Ωb h2 Rph τ0 = 3.9 × 10 fHI , (7.2) 104 K 10 0.022 10kpc 139 where fHI is the fraction of HI. All other parameters in (7.2) is taken from the concordance ΛCDM mode. For these objects the relation between dimensionless η and physical time t is given by µ ¶1/2 µ ¶−3 µ ¶−1 −1 T 1+z Ωb h2 t = 5.4 × 10−2 fHI η, yr. (7.3) 104 K 10 0.022 The 21 cm emission rely on the W-F effect. On the other hand, the time-scale of the evolution of the 21 region is short. The effect of dust on the time-scales of Lyα evolution should be considered. We have not considered the Lyα photons produced by the recombination in the ionized halo. If the halo is optical thick, photons from the recombination will also be thermalized. The information of where the photon comes from will be forgotten during the thermalization. Therefore, photons from recombination should not show any difference from those emitted from central sources. Only the photons formed very close to the boundary of the halo will not be thermalized, and may yield different behavior. Chapter 8 Angular distribution of Lyα resonant photons emergent from optically thick medium This chapter will study the angular distribution of Lyα photon transferring in an op- tically thick medium. Previous methods are based on the Eddington approximation and the evolution of the angular distribution is completely ignored. However, the evolution of angular distribution actually is significant. In a thermalized or statistical equilibrium state, the angular distribution of photons should be isotropic, regardless of the initial angular distribution. Therefore, one can expect that the angular distri- bution of Lyα photons with resonant frequency ν0 should be isotropic. On the other hand, the angular distribution of photons with frequency different from ν0 might be anisotropic, as those photons are not involved in the evolution of thermalization or statistical equilibrium. Consequently, the angular distributions of Lyα resonant photons from optically thick medium should be frequency-dependent. It definitely cannot be described by the Eddington approximation. The evolution of the angular distribution of resonant photons is not trivial. We still use the WENO solver, and solve the photon transfer in both frequency and angular spaces. 140 141 8.1 Transfer equations of resonant photons with- out dust We solve the radiative transfer equation of Lyα resonant photon in a spherically symmetric medium containing neutral HI. For simplicity, we ignore the effect of dust (i.e. κ(x) = 0), and (6.1) is ∂I ∂I (1 − µ2 ) ∂I ∂I +µ + −γ = ∂η ∂r r ∂µ ∂x Z −φ(x; a)I + R(x, µ, x0 , µ0 ; a)I(η, r, x0 , µ0 )dx0 dµ0 /2 + S. (8.1) The boundary and initial conditions are still given in (6.13), (6.17) and (6.19). 8.1.1 Test with Field’s analytical solution We first test the WENO solver with analytical solutions. Assuming that the specific intensity and source S are homogeneous in the r and µ space, i.e. I(η, r, x, µ) is independent of variables r and µ. (8.1) becomes Z ∂J ∂J −γ = −φ(x)J + R(x, x0 )J(η, x0 )dx0 + S, (8.2) ∂η ∂x where Z 1 J(η, x) = I(η, r, x, µ)dµ. (8.3) 2 2 Take γ = 0, Voigt parameter a = 0, the source S = φ(x) = π −1/2 e−x , and the initial radiative field I(x, η = 0) = 0. The time-dependent solution of(8.2) is [44, 94]. 2 J(x, η) = π −1/2 [1 − exp(−ηe−x )] (8.4) Z ∞ 2 2 2 + ew [1 − (1 + ηe−w ) exp(−ηe−w )]erf(w)dw. x 142 Our solver seeks a solution I from (8.1). One can then give J via (8.3). It is interesting to see whether the solution to (8.4) can be reproduced, if we also assume that the source S in (8.1) is spatially homogeneous, but µ-dependent, i.e. S = 2 Θ(µ)φ(x) = Θ(µ)π −1/2 e−x where Θ(µ) describes the angular distribution of photons from the source. We consider both an isotropic source 2 S = π −1/2 e−x , −1 ≤ µ ≤ 1, (8.5) and an anisotropic source as follows,   2(n + 1)µn π −1/2 e−x2 , 0 < µ ≤ 1, S= (8.6)  0, −1 ≤ µ ≤ 0, where n is taken to be a positive integer.Obviously, the larger the n, the stronger the emission in the direction µ = 1. The factor 2(n + 1) is for normalization: R 1 1 2 0 2(n + 1)µn dµ = 1. The numerical results with sources (8.5) and (8.6) with n = 4 and 6 are shown in Figure 8.1. It is expected that the numerical solution with source (8.5) (the left =10000 =10000 =10000 103 103 103 =1000 =1000 =1000 102 102 102 =100 =100 =100 10 1 =10 101 =10 101 =10 J(X) J(X) J(X) 100 =1 100 =1 100 =1 -1 -1 -1 10 =0.1 10 =0.1 10 =0.1 -2 -2 10 -2 =0.01 10 =0.01 10 =0.01 -3 -3 10-3 -5 0 5 10 -5 0 5 10 -5 0 5 X X X Figure 8.1: The WENO numerical solutions (solid lines) of equation (8.1) assuming 2 2 the sources S is a) S = π −1/2 e−x for all µ (left); b) S = 10µ4 π −1/2 e−x (middle); and 2 c) S = 14µ6 π −1/2 e−x (right) for µ > 0 and S=0 for µ < 0. The Field’s analytical solution is shown with dot lines. 143 panel of Figure 8.1) should follow the analytical solution (8.4) well, as the isotropic source is the same as that used to find the analytical solution. It is interesting to see that the WENO solutions of n = 4 (middle panel) and n = 6 (right panel) also follow the analytical solution to (8.4) well. It seems to indicate that the evolution of the frequency space is independent of the µ-space. 8.1.2 Time scale of the statistical equilibrium of the angular distribution A remarkable feature of the solutions of Figure 8.1 is a flat plateau in the range |x| < 2 at time η > 100. The flat plateau is caused by the Wouthuysen-Field local thermalization of frequency distribution of resonant photon [113, 44, 45]. The flat plateau actually is the Boltzmann statistical equilibrium distribution around x = 0 when the atomic mass is infinite. If the mass is finite, i.e. considering the recoil in the re-distribution functions (6.3) or (6.4), the flat plateau will become e−2bx , where b = hν0 /mvT c, which is the local Boltzmann distribution required by the Wouthuysen-Field effect [88]. The resonant scattering between photons and HI atoms leads to the Boltzmann distribution of the photon frequency distribution around x = 0 with the temperature equal to that of HI atoms. When resonant photons undergo the local thermalization in the frequency space, the angular distribution should be approaching statistical equilibrium. The anisotropic µ-distributions have to evolve to be isotropic (statistical equilibrium). We calculate all the µ-distributions at the time η corresponding to the three panels of Figure 8.1. The result is plotted in Figure 8.2. The µ-distribution of left panel is always isotropic. This is expected as the source is isotropic, which is already in the state of statistical equilibrium. The middle and right panels of Figure 8.2 show the evolution of an anisotropic µ- distribution equation (8.6). The time scale for approaching an isotropic distribution seems to be independent of the anisotropy of sources. It is always equal to about 144 3 =10000 3 =10000 3 =10000 10 10 10 =1000 =1000 =1000 102 =100 102 102 =100 =100 1 1 10 =10 10 =10 101 =10 I( ,X, ) I( ,X, ) I( ,X, ) 10 0 =1 10 0 10 0 10 -1 =0.1 10 -1 =1 10 -1 =1 10 -2 =0.01 10 -2 10 -2 -3 -3 =0.1 -3 =0.1 10 10 10 -4 -4 -4 10 10 10 =0.01 =0.01 10-5 -0.5 0 0.5 10-5 -0.5 0 0.5 10-5 -0.5 0 0.5 Figure 8.2: The WENO numerical solutions of angular distributions from equation 2 (8.1) at x = 0, assuming the sources S are a) S = π −1/2 e−x for all µ (left) and b) 2 2 S = 10µ4 π −1/2 e−x (middle); and c) S = 14µ6 π −1/2 e−x (right) for µ > 0 and S=0 for µ < 0. η ∼ 100 for both n = 4 and n = 6, i.e. the µ-distribution will become isotropic after 100 times of resonant scattering. This time scale is about the same as that of the W-F thermalization (Figure 8.1). Therefore, the thermalization in the frequency space and the isotropic distribution in the µ-space are realized at about the same time. 8.2 Precision of the Eddington approximation 8.2.1 Equations of the Eddington approximation Following the same analysis as in Section 6.1.3, (8.1) yields the equations of j and f as Z ∂j ∂f ∂j + = −φ(x; a)j + R(x, x0 ; a)jdx0 + γ + r2 S, (8.7) ∂η ∂r ∂x ∂f 1 ∂j 2 j ∂f + − = −φ(x; a)f + γ , (8.8) ∂η 3 ∂r 3 r ∂x with initial and boundary conditions (6.14), (6.18) and (6.20). 145 8.2.2 Profiles of j and f Figure 8.3 does show small differences between the solutions with and without the Eddington approximation, even though both solutions are given by the same source. The difference comes from the contribution of the terms of l > 2 in the Legendre expansion. The difference between the profiles with and without the Eddington approximation becomes smaller when the time η is larger. It is because larger η corresponds to larger optical depth. The Eddington approximation generally is good for optically thick medium. 0 0 10 10 =500 =300 =2000,3000 -1 -1 10 10 f( ,r,x) j( ,r,x) =200 =1000 10-2 10-2 =500 -3 -3 10 10 -5 0 5 -2 -1 0 1 2 x x Figure 8.3: A comparison of the solutions j and f with the Eddington approximation (dashed curves) and the solutions of f without the Eddington approximation (solid curves). Relevant parameters are r = R = 102 , and a = 10−3 . Now we consider different sources. We re-do the solutions of j and f with equation (8.1) by taking S = δ(µ − 21 ) and S = δ(µ − 1). We use polynomials of degree 6 to approximate the delta sources. The results are given in Figure 8.4, which shows the same shape of the profiles. That is, the profiles of j and f are not affected by the angular distribution of photons from the source. It is probably because the µ-distribution quickly evolves into the statistical equilibrium state, the initial 146 anisotropy of the µ distribution is forgotten. 100 0.6 =500 0.4 =300 0.2 10-1 f( ,r,x) j( ,r,x) =200 =2000,3000 -2 10 =1000 =500 10-3 -5 0 5 -2 0 2 x x Figure 8.4: A comparison of the solution s j and f with S=δ(µ − 1/2) and δ(µ − 1). Relative parameters are r = R = 100, a = 10−3 . 8.3 Angular distributions 8.3.1 Frequency dependence Although the Eddington approximation is acceptable when calculating the profile of Lyα photons in the frequency space, it fails in the µ-space. The result in Figure 8.2 shows that the µ-distribution is isotropic at frequency ν0 . On the other hand, the µ-distribution will no longer be isotropic at frequency |x| ≥ 2, because photons of |x| ≥ 2 have not undergone scattering. Consequently, the angular distribution of photons emerging from optically thick halo should be frequency(energy)-dependent. We calculate the µ-distribution of photons from a halo with R = 500 with the central source given by (6.15), i.e. photons from the source can be described by the Eddington approximation (6.9). The result is shown in Figure 8.5. The µ distributions at frequencies x = 0 and 0.8 are basically straight lines in the whole 147 8 x=2.4 6 I( ,r,x, ) 4 x=1.6 2 x=0 x=0.8 -0.5 0 0.5 Figure 8.5: The µ distribution of photons emergent from a halo with radius R = 500. The frequencies are x = 0.0, 0.8, 1.6 and 2.4. The relevant parameters of the calculation are η = 1.2 × 104 , γ = 0, and a = 10−3 . 148 range −1 ≤ µ ≤ 1. That is, I can be described by the Eddington approximation equation (6.9). At x = 1.6, the µ-distribution begins to deviate from a straight line, i.e. deviating from an Eddington approximation. At x = 2.4, the µ-distribution shows a very sharp spike at µ = 1. That is, the angular distribution of photons with frequency at the two peaks (Figure 8.3) is significantly different from isotropic, but is dominated by photons of µ = 1. This result is consistent with the “single shot picture” [2, 12], in which photons with frequency |x| < 2 mainly undergo a diffusion in the frequency space; once a photon diffuses to |x| ≥ 2, it will take “single longest excursion” to leave for outside of the halo. Therefore, the two peaks of the flux f at frequency x± ' ±(2 − 3) are dominated by photons from “single longest excursion” photons, of which µ ∼ 1. 8.3.2 Dependence of the initial anisotropy The source of Figure 8.5 given by (6.15) has Θ(µ) = 6µ (µ > 0), which is linear in µ. We now consider sources with higher anisotropy with Θ(µ) given by   2(n + 2)µn , 0 < µ ≤ 1, Θ(µ) = (8.9)  0, µ < 0. When the integer n is large, Θ(µ) is similar to a δ function δ(µ−1), i.e. most photons are in the direction µ = 1. We repeat the calculation of Figure 8.5, but using the source equation (8.9) with n = 1, 2, 4, 6 and 8. The result is plotted in Figure 8.6. It is interesting to see that the µ-distributions are independent of n, but depend on x. It is easy to explain the n-independence of the two top panels of Figure 8.6, both of which have frequency |x| ≤ 2. In this frequency range, the evolution of the specific intensity I is governed by the local thermalization of x-space and entropy increasing of µ- space. These processes lead to the Boltzmann distribution in the energy space and 149 4 x=1.6 3 x=2.4 I( ,r,x, ) x=0.8 2 x=0 1 -0.5 0 0.5 Figure 8.6: The µ-distribution of halo with radius R =100 and source equation (8.9) with n=1, 2, 4, 6, and 8. For each x, the curves for different n overlap with each other. The frequencies are taken to be x = 0, 0.8, 1.6, 2.4. The parameters of the halos are η = 3500, γ = 0, and a = 10−3 . isotropic distribution in the angular space, regardless of the initial distributions in either frequency- or angular spaces. In other words, the initial distribution is forgotten during the local thermalization and approaching statistical equilibrium. However, the mechanism of the local thermalization and approaching isotropic distribution seems to be unable to explain why the two curves at x=1.6 and x=2.4 of Figure 8.6 also show n-independence. The µ-distribution of these two curves of Figure 8.6 are highly anisotropic. Therefore, they do not have to be the result of the local thermalization and approaching statistical equilibrium. Why do they also show the behavior of forgetting the initial angular distributions? The reason is as follows. In the first phase of resonant photon evolution, Lyα photons are trapped in the range of |x| ≤ 2 within the time scales of a few tens or hundred scattering [89]. The trapped photons have already forgotten their initial state. On the other hand, photons with |x| ≥ 2 mostly come from the diffusion of trapped photons from 150 |x| ≤ 2 to |x| ≥ 2 [88]. Thus, all photons of |x| ≥ 2 emerging from the optically thick halo essentially have the same initial condition, given by the |x| space diffusion of trapped photons. Therefore, the initial distributions before they are trapped have been forgotten. This property can also be seen in Figure 8.2, in which, although the sources of the middle and right panels are different from each other, the behaviors of the time-evolution of the µ-distribution are about the same. This result also implies that it is impossible to recover the information of the distribution of photons emitted by the central source. 8.3.3 Collimation of photons of the double peaks A common feature of Figure 8.6 is to show a very sharp spike at µ ∼ 1 when |x| = 1.6 and |x| = 2.4, corresponding to the double peaks of Figures 8.3 and 8.4. Therefore, the spiky distribution of µ indicates that the photons with frequency at the double peaks have formed a forward beam. In order to measure the angular size of the µ = 1 spikes, we fit the µ-distributions of Figure 8.5 at x = 1.6 and x = 2.4 with polynomials of µ. We find that both curves can be well fitted with polynomials of µ having leading terms Aµ16 + Bµ15 + ..., A, B being fitting coefficients. The terms of either µ16 or µ15 are much sharper than the central source equation (8.9) µn with n ≤ 6. Therefore, the radiative transfer at the double peaks of frequency space plays the role of forward collimator. It made the photons form forward beams. If we define the spread angle β of the forward beam as the angle of half intensity, this number can be estimated by cos16 β = 1/2, and therefore, β ∼ 0.29rad. This result is again consistent with the “single shot picture”. The double peaks mainly consist of photons from a single shot, which moves in the forward direction. 151 8.3.4 Large halo We calculate the µ-distribution of Lyα photons in a halo with large radius R = 1000, and the central source is given by equation (8.9) and n = 6. The results are given in Figure 8.7, which shows the dependence of the µ distribution on the radial variable r in the halo. 101 101 101 x=0 x=0.8 x=1.6 x=0 x=0.8 x=1.6 x=0 x=0.8 x=1.6 I( ,r,x, ) I( ,r,x, ) I( ,r,x, ) 100 100 100 x=2.4 x=2.4 x=2.4 10-1 10-1 10-1 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 Figure 8.7: The µ-distributions at radial positions r = 100 (left), 300 (middle), 500 (right) of a halo with radius R = 1000. The source is given by equation (8.9) and n = 6. The frequencies are taken to be x = 0, 0.8, 1.6, 2.4. The dotted curves in the middle and right panels are µ-distributions at r = 100 x = 2.4. Other relevant parameters are η = 1.2 × 104 , γ = 0 and a = 10−3 . Although the photons from the source of n = 6 are highly anisotropic, all the µ-distributions of x = 0.0 and 0.8 at r = 100, 300, 500 are straight lines. That is, the specific intensity I can be well approximated by the Eddington approximation equation (6.9). This result is consistent with Section 8.3.2. The r-dependence of the µ-distribution of |x| = 2.4 photons is also consistent with the result of section 8.3.3: the larger the r, the sharper the µ-distribution. The r transfer leads to the collimation. The behavior of r-dependence of the µ-distribution at x = 1.6 is very different from that of x = 0.0, 0.8, and 2.4. The µ distribution at r = 100 is about the same as Figure 8.6, i.e. it has undergone an evolution of forward collimation, having a sharp spike at µ = 1. However, the µ distribution will no longer show a spike at 152 r = 300 and 500. When r is equal to or less than about 200, the r-dependence of the µ distribution is similar to the |x| = 2.4 photons. However, when r ≥ 200, the r-dependence of the µ distribution is similar to the |x| = 0 and 0.8 photons. This is because the optical depth at |x| = 1.6 is larger than |x| = 2.4, the “single shot picture” is working well at r ∼ 100 for both |x| = 1.6 and 2.4. However, at r ≥ 200, the single shot picture can still work well for |x| = 2.4, but not so well for |x| = 1.6. (x=2.0)=0.930 r=300 (x=2.0)=0.930 r=300 8 8 (x=2.4)=0.558 6 r=500 (x=2.4)=3.139 6 r=300 I( ,r,x, ) 4 r=500 (x=2.0)=5.232 I( ,r,x, ) 4 r=200 (x=2.4)=0.372 r=300 (x=1.6)=13.138 r=100 (x=2.4)=0.186 2 2 r=200 (x=0.0)=112.71 0 -0.5 0 0.5 0 -0.5 0 0.5 Figure 8.8: The µ-distributions with respect to the effective optical depth. Relevant parameters are η = 1.2 × 104 , γ = 0 and a = 10−3 . We consider the evolution of the angular distribution with respect to effective radial optical depth τr (x) = τ0 φ(x, a). The result is plotted in Figure 8.8. The µ-distribution is isotropic if the effective radial optical depth is large and it will no longer be isotropic when the depth becomes small. If we define the transition between isotropic and anisotropic µ-distribution occurs when I(r, x, 1) = 2I(r, x, −1), then at the transition, the critical optical depth is τcrit ≈ 5, as can be found from table 8.1. 153 Table 8.1: critical effective optical depth τcrit r 100 200 300 400 500 x 1.6 1.8 1.8 2.0 2.0 τr (x) 4.41 4.46 6.69 4.17 5.21 8.3.5 Effect of anisotropic scattering All calculations in the previous sections are based on the re-distribution function equation (6.3), which consider only isotropic scattering. If we consider dipole scatter- ing, it seems to introduce a new factor leading to anisotropy and yield new anisotropic behavior. However, HI atoms are in thermal equilibrium, and their distribution is isotropic. The dipole scattering, as average, does not contain any parameter of spe- cific direction. It will not add any anisotropic behavior. Therefore, all conclusions in the previous sections should still hold. 8.4 Conclusion The transfer of Lyα resonant photons from a central source in a halo consisting of HI generally is considered as a problem of radiative transfer in an optically thick medium. However, the “optically thick medium” assumption is true only when the frequency of Lyα photons lies in a narrow range |x| ≤ 2. The cross section of resonant scattering is very sensitive to the photon frequency. It quickly becomes small when the frequency of Lyα photons has only a small deviation from the range |x| ≤ 2. For those photons, the halo is optically moderate thick, or even thin. Therefore, in order to understand the transfer of Lyα photons with frequency around the resonant peak, we need to find the solutions of the integro-differential equation (8.1) in optically thick as well as moderate thick and even thin medium. That is, although the halo is optically thick for resonant photons, one should not treat (8.1) by using the condition 154 of optical thick. To find solution of (8.1) with desired precision in frequency ranges of optically thick as well as moderately thick, the algorithm must handle the extremely flat distribution (|x| < 2) and its sharp boundary (|x| ∼ 2) of I. These features can be properly captured by the state-of-the-art numerical method, WENO scheme, as it has high order of accuracy and good convergence in capturing discontinuities as well as being superior to piecewise smooth solutions containing discontinuities. The WENO solver has been shown to be powerful to solve the integro-differential equation of radiative transfer of resonant photons Lyα. In this chapter, we develop the WENO algorithm to be able to solve the integral-differential equation (8.1) in frequency and angular space simultaneously. We have first shown that the Eddington approximation can yield reasonable results of the frequency profile of photons emergent from optically thick halos. Since the Eddington approximation assumes I is linearly depends on µ, all the physics of the angular distribution of Lyα photons are missing. A cost of the Fokker-Planck equations is also to ignore all the effects of the evolutions of angular distribution. The physics of the evolution of the angular distribution is rich. As has been known, resonant scattering couples the transfers of resonant photons in the phys- ical space and the frequency space. We show, in this chapter, that the resonant scattering leads to the coupling between the evolutions of resonant photons in the frequency space and the angular space as well. The evolution of the resonant photon distribution in the µ space is strongly dependent on the frequency. Photons with frequency |x| ≤ 2 undergo the procedure of approaching statistical equilibrium, and their angular distribution is isotropic after a few tens or hundred scattering, regard- less of whether the initial angular distribution is isotropic. On the other hand, the angular distribution of photons with |x| ≥ 2 is strongly anisotropic, even if the initial angular distribution is isotropic. An interesting feature is that the anisotropic angular distributions at frequency 155 |x| ∼ 2 are independent of the initial angular distributions. Different initial angular distributions yield the same anisotropic angular distributions after a few tens or hundred scattering. This is because photons at frequency |x| ∼ 2 are not directly from the source, but come from the trapped photons within |x| ≤ 2, for which the initial distributions have been forgotten. Therefore, it seems to be impossible to recover the property of the source with the observed µ-distribution of Lyα photons either in the range of |x| ≤ 2 or in |x| ≥ 2. Another interesting feature of an optically thick halo is the collimation of photons with frequencies of the double peaks. This is also because photons trapped in |x| ≤ 2 are thermal. When the trapped photons diffuse to |x| ≥ 2, they have two possible fates. One is to get out of the holes by a single shot if photons move forward. If a photon has not taken a single shot, the resonant scattering will lead it back to the region of |x| ≤ 2. Therefore, photon transfer in optically thick medium is a collimator. Although photons stored in an optically thick halo are thermal, the µ- distribution is isotropic and the double peak only picks up photons of a single shot, i.e. moving forward. Bibliography [1] T. Adams, The escape of resonance-line radiation from extremely opaque media, The Astrophysical Journal, 174 (1972), 439-448. [2] T. Adams, The mean photon path length in extremely opaque media, The As- trophysical Journal, 201 (1975), 350-351. [3] T. Adams, D.G. Hummer, G.B. Rybicki Numerical evaluation of the redistribu- tion function RII−A (x, xy) and of the associated scattering integral, Journal of Quantitative Spectroscopy and Radiative Transfer, 11 (1971), 1365-1376. [4] S. Adjerid and M. Baccouch, Asymptotically exact a posteriori error estimates for a one-dimensional linear hyperbolic problem, Applied Numerical Mathemat- ics 60 (2010), 903-914. [5] S. Adjerid, K. Devine, J. Flaherty and L. Krivodonova, A posteriori error esti- mation for discontinuous Galerkin solutions of hyperbolic problems, Computa- tional Methods in Applied Mechanics and Engineering 191 (2002), 1097-1112. [6] S. Adjerid and T. Massey, Superconvergence of discontinuous Galerkin solu- tions for a nonlinear scalar hyperbolic problem, Computer Methods in Applied Mechanics and Engineering 195 (2006), 3331-3346. [7] S. Adjerid and T. Weinhart, Discontinuous Galerkin error estimation for linear symmetric hyperbolic systems, Computer Methods in Applied Mechanics and Engineering 198 (2009), 3113-3129. 156 157 [8] S. Adjerid and T. Weinhart, Discontinuous Galerkin error estimation for lin- ear symmetrizable hyperbolic systems, Mathematics of Computations 80 (2011), 1335-1367. [9] L.W. Avery and L. House, An investigation of resonance-line scattering by the Monte Carlo technique, The Astrophysical Journal, 152 (1968), 493-508. [10] C. Berthon, M. Breuss and M.-O. Titeux, A relaxation scheme for the approxi- mation of the pressureless Euler equations, Numerical Methods for Partial Dif- ferential Equations, 22, (2006), 484-505. [11] G. Blanc et. al., The HETDEX pilot survey. II. the evolution of the Lyα escape fraction form the ultraviolet slope and luminosity function of 1.9¡z¡3.8 LAEs, The Astrophysical Journal, 736 (2010), 31(21). [12] J.R.M. Bonilha et. al., Monte Carlo calculation for tesonance scattering with absorption or differential expansion, The Astrophysical Journal, 233 (1979), 649-660. [13] F. Bouchut, On zero pressure gas dynamics, Advances in Kinetic Theorey and Computing, Would Scientific, River Edge, NJ, (1994), 171-190. [14] F. Bouchut and F. James Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness, Communications in Partial Differen- tial Equations, 24, (1999), 2173-2189. [15] F. Bouchut, S. Jin and X. Li, Numerical approximations of pressureless and isothermal gas dynamics, SIAM Journal on Numerical Analysis, 41, (2003), 135-158. [16] L. Boudin and J. Mathiaud, A numerical scheme for the one-dimensional pres- sureless gases system, Numerical Methods for Partial Differential Equations, 28, (2006), 1729-1746. 158 [17] J.H. Bramble and A.H. Schatz, High order local accuracy by averaging in the finite element method, Mathematics of Computation, 31 (1977), 94-111. [18] Y. Brenier and E. Grenier, Sticky particles and scalar conservation laws, SIAM Journal on Numerical Analysis, 35, (1998), 2317-2328. [19] C. Canuto, F. Fagnani and P. Tilli, An Eulerian approach to the analysis of Krause’s consensus models, SIAM Journal on Control and Optimization, 50, (2012), 243-265. [20] J. Carrillo, 2D semiconductor device simulations by WENO-Boltzman schemes:Efficiency, boundary conditions and comparison to Monte Carlo meth- ods, Journal of Computational Physics, 214 (2006), 55-80. [21] A. Chalaby, On convergence of numerical schemes for hyperbolic conservation laws with stiff source terms, Mathematics of Computation, 66 (1997), 527-545. [22] G.-Q. Chen and H. Liu, Formation of δ-shocks and vacuum states in the vanish- ing pressure limit of solutions to the Euler equations for isentropic fluids, SIAM Journal on Mathematical Analysis, 34, (2003), 925-938. [23] J. Cheng and C.-W. Shu, A cell-centered Lagrangian scheme with the preser- vation of symmetry and conservation properties for compressible fluid flows in two-dimensional cylindrical geometry, Journal of Computational Physics, 229, (2010), 7191-7206. [24] J. Cheng and C.-W. Shu, Improvement on spherical symmetry in two- dimensional cylindrical coordinates for a class of control volume Lagrangian schemes, Communications in Computational Physics, 11, (2012), 1144-1168. [25] Y. Cheng and C.-W. Shu, Superconvergence and time evolution of discontinuous Galerkin finite element solutions, Journal of Computational Physics, 227 (2008), 9612-9627. 159 [26] 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 (2010), 4044-4072. [27] A. Chertock, A. Kurganov and Y. Rykov, A new sticky particle method for pressureless gas dynamics, SIAM Journal on Numerical Analysis, 45, (2007), 2408-2441. [28] P.G. Ciarlet. Finite Element Method for Elliptic Problems, North-Holland, Am- sterdam, 1978. [29] B. Cockburn et. al., In advanced numerical approximation of nonlinear hyper- bolic equations, ed. A. Quarteroni (Lecture Notes in Mathematics, Vol. 1697 Berlin: Springer), 450 [30] B. Cockburn, A Introduction to the discontinuous Galerkin methods for convec- tion dominated problems, High-Order Method for Computational Physics (T. Barth and H. Deconink, eds.), Lecture Notes in Computational Science and Engineering, vol. 9, Springer-Verlag Berlin Heidelberg, 1999, 69-224. [31] B. Cockburn and J. Guzm´an, Error estimate for the Runge-Kutta discontinuous Galerkin method for the transport equation with discontinuous initial data, SIAM Journal on Numerical Analysis, 46 (2008), 1364-1398. [32] B. Cockburn, S. Hou and C.-W. Shu, The Runge-Kutta local projection discon- tinuous Galerkin finite element method for conservation laws, IV: the multidi- mensional case, Mathematics of Computation, 54 (1990), 545-581. [33] 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 (1989), 90-113. 160 [34] 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 (2003), 577-606. [35] 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 (1989), 411-435. [36] B. Cockburn and C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws, V: multidimensional system, Journal of Computational Physics, 35 (1998), 199-224. [37] M. Dijkstra, Z. Haiman, M. Spaans, Lyα radiation from collapsing protogalaxies, The Astrophysical Journal, 649 (2006), 14-36. [38] M. Dijkstra, A. Loeb, Lyα blobs as an observational signature of cold accretion streams into galaxies, Monthly Notices of the Royal Astronomical Society, 400 (2009), 1109-1120 [39] B.T. Draine, Scattering by interstellar dust grains. I. optical and ultraviolet, The Astrophysical Journal, 598 (2003), 1017-1025. [40] B.T. Draine and H.M. Lee Optical properties of interstellat graphite and silicate grains, The Astrophysical Journal, 285 (1984), 89-108. [41] W. E., Yu. G. Rykov and Ya. G. Sinai, Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics, Communications in Mathematical Physics, 177 (1996), 349-380. [42] L.Z. Fang The zeroth law of thermodynamics of the photon-hydrogen system and 21cm cosmology, International Journal of Modern Physics D, 18 (2009), 1943-1954. 161 [43] M. Fardal Cooling radiation and the Lyα luminosity of forming galaxies, The Astrophysical Journal, 562 (2001), 605-617. [44] G.B. Field Excitation of the hydrogen 21-cm line, Proceedings of the IRE, 46 (1958), 240-250. [45] G.B. Field The time relaxation fo a resonance-line profile, The Astrophysical Journal, 129 (1959), 551-564. [46] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review, 43 (2001), 89-112. [47] J.M. Greenberg, A.Y. Leroux, R. Baraille and A. Noussair, Analysis and ap- proximation of conservation laws with source terms, SIAM Journal on Numerical Analysis, 34 (1997), 1980-2007. [48] Z. Haiman and M. Spaans, Models for dusty Lyα emitters at high redshift, The Astrophysical Journal, 518 (1999), 138-144. [49] Z. Haiman, M. Spaans and E. Quataert, Lyα cooling radiation from high-redshift halos, The Astrophysical Journal, 537 (2000), L5-L8. [50] M. Hansen and S.P. Oh, Lyα radiative transfer in a multiphase medium, Monthly Notices of the Royal Astronomical Society, 367 (2006), 979-1002. [51] J.P. Harrington, The scattering of resonance-line radiation in the limit of large optical depth, Monthly Notices of the Royal Astronomical Society, 162 (1973), 43-52. [52] M. Hayes, Escape of about five per cent of Lyα photons froom high-redshift star- forming galaxies, Nature, 464 (2010), 562-565. [53] M. Hayes, On the redshift evolution of the Lyα escape fraction and the dust content of galaxies, The Astrophysical Journal, 730 (2011), 8(13). 162 [54] L.G. Henyey and J.L. Greenstein, Diffuse radiation in the galaxy, The Astro- physical Journal, 93 (1941), 70-83. [55] S. Hou and X.-D. Liu, Solutions of multi-dimensional hyperbolic systems of conservation laws by square entropy condition satisfying discontinuous Galerkin method, Journal of Scientific Computing, 31 (2007), 127-151. [56] D.G. Hummer, Non-coherent scattering: I. The redistribution function with Doppler broadening, Monthly Notices of the Royal Astronomical Society, 125 (1962), 21-37. [57] D.G. Hummer, The Voigt function: An eight-significant-figure table and gener- ating procedure, Memoirs of the Royal Astronomical Society, 70 (1965), 1-32. [58] D.G. Hummer, Non-coherent scattering: VI. Solution of the transfer problem with a frequency-dependent source function, Monthly Notices of the Royal As- tronomical Society, 145 (1969), 95-120. [59] D.G. Hummer and P.B. Kunasz, Energy loss by resonance line photons in an absorbing medium, The Astrophysical Jornal, 236 (1980), 609-618. [60] G.-S. Jiang and C.-W. Shu, On a cell entropy inequality for discontinuous Galerkin methods, Mathematics of Computation, 62 (1994), 531-538. [61] G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics, 126 (1996), 202-228. [62] F. John, Partial different equations, Springer-Verlag, New York, 1971 [63] C. Johnson, U. N¨avert and J. Pitk¨aranta, Finite element methods for linear hyperbolic problems, Computer Methods in Applied Mechanics and Engineering, 45 (1984), 285-312. 163 [64] C. Johnson and J. Pitk¨aranta, An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation, Mathematics of Computation, 46 (1986), 1-26. [65] C. Johnson, A. Schatz and L. Walhbin, Crosswind smear and pointwise errors in streamline diffusion finite element methods, Mathematics of Computation, 49 (1987), 25-38. [66] B. Koren, A robust upwind discretization method for advection, diffusion and source terms, Notes on Numerical Fluid Mechanics, 45, Vieweg, Braunschweig (1993), 117-138. [67] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon and J.E. Flaherty, Shock detection and limiting with discontinuous Galerkin methods for hyperbolic con- servation laws, Applied Numerical Mathematics, 48 (2004), 323-338. [68] P. Kurasov, Distribution theory for discontinuous test functions and differential operators with generalized coefficients, Journal of Mathematical Analysis and Applications, 201 (1996), 297-323. [69] M. Latif, et. al., Lyα emission from the first galaxies: signatures of accretion and infall in the presence of line trapping, Monthly Notices of the Royal Astro- nomical Society: Letters, 413 (2011), L33-L37. [70] P. Laursen and J. Sommer-Larsen, Lyα resonsant scattering in young galax- ies: Predictions from cosmological simulations, The Astrophysical Journal, 657 (2007), L69-L72. [71] P. Laursen, J. Sommer-Larsen and A. Andersen, Lyα radiative transfer with dust: escape fractions from simulated high-redshift galaxies, The Astrophysical Jornal, 704 (2009), 1640-1656. 164 [72] R.J. LeVeque and H.C. Yee, A study of numerical methods for hyperbolic con- servation laws with stiff source terms, Journal of Computational Physics, 86 (1990), 187-210. [73] A. Li, B.T. Draine, On ultrasmall silicate grains in the diffuse interstellar medium, The Astrophysical Jornal, 550 (2001), L213-L217. [74] J. Liu, et. al., 21cm signals from early ionizing sources, The Astrophysical Jor- nal, 663 (2007), 1-9. [75] A. Loeb and G.B. Rybicki, Scattered Lyα radiation around sources before cos- mological reionization, The Astrophysical Journal, 524 (1999), 527-535. [76] D. Neufeld, The transfer of resonance-line radiation in static astrophysical me- dia, The Astrophysical Jornal, 350 (1990), 216-241. [77] D. Neufeld, The escape of Lyα radiation from a multiphase interstellar medium, The Astrophysical Jornal, 370 (1991), L85-L88. [78] A. Noussair, Analysis of nonlinear resonance in conservation laws with point sources and well-balanced scheme, Studies in Applied Mathematics, 104 (2000), 313-352. [79] P. de Oliveira and J. Santos, On a class of high resolution methods for solving hyperbolic conservation laws with source terms, Applied Nonlinear Analysis, 432- 445, (Ad´ elia Sequeira, Hugo Beir˜ ao da Veiga, and Juha Hans Videman, eds.), Kluwer academic publishers, New York, Boston, Dordrecht, London, Moscow, 1999. [80] D.E. Osterbrock, The escape of resonance-line radiation from an optically thick nebula, The Astrophysical Jornal, 135 (1962), 195-216. [81] Y. Pei, Interstellar dust from the milky way to the magellanic clouds, The As- trophysical Jornal, 395 (1992), 130-139. 165 [82] M. Pierleoni, A. Maselli and B. Ciardi, Crashα: coupling continuum and line ra- diative transfer, Monthly Notice of the Royal Astronomical Society, 393 (2009), 872-884. [83] J. Qiu, et. al., A WENO algorithm for the radiative transfer and ionized sphere at reionization, New Astronomy, 12 (2006), 1-10. [84] J. Qiu, et. al., A WENO algorithm of the temperature and ionization profiles around a point source, New Astronomy, 12 (2007), 398-409. [85] J. Qiu, et. al., A WENO algorithm for the growth of ionized regions at the reionization epoch, New Astronomy, 13 (2008), 1-11. [86] W.H. Reed and T.R. Hill, Triangular mesh for the Neutron transport equa- tion, Los Alamos Scientific Laboratory Report LA-UR-73-479, Los Alamos, NM, 1973. [87] I. Roy, et. al., A WENO algorithm for radiative transfer with resonant scattering and the Wouthuysen-Field coupling, New Astronomy, 14 (2009), 513-520. [88] I. Roy, et. al., Time evolution of Wouthuysen-Field coupling, The Astrophysical Journal, 694 (2009), 1121-1130. [89] I. Roy, et. al., Wouthuysen-Field coupling in the 21 cm region around high- redshift sources, The Astrophysical Journal, 704 (2009), 1992-2003. [90] I. Roy, C.-W. Shu and L.-Z. Fang, Resonant scattering and Lyα radiation emer- gent from nertral hydrogen halos, The Astrophysical Journal, 716 (2010), 604- 614. [91] J.K. Ryan and B. Cockburn, Local derivative post-processing for the discontinu- ous Galerkin method, Journal of Computational Physics 228 (2009), 8642-8664. 166 [92] J.K. Ryan and C.-W. Shu, On a one-sided post-processing technique for the dis- continuous Galerkin methods, Methods and Applications of Analysis, 10 (2003), 295-308. [93] G.B. Rybicki, Improved Fokker-Planck equation for resonance-line scattering, The Astrophysical Journal, 647 (2006), 709-718. [94] G.B. Rybicki and I. dell’Antonio, The time development of a resonance line in the expanding universe, The Astrophysical Journal, 427 (1994), 603-617. [95] G.B. Rybicki and A.P. Lightman, Radiative Processes in Astrophysics, (New York: Wiley). [96] Yu. G. Rykov, Propagation of shock wave type singularities in equations of two- dimensional zero-pressure gas dynamics, Mathematical notes, 66 (1999), 628- 635. [97] Yu. G. Rykov, On the nonhamiltonian character of shocks in 2-D pressureless gas, Bollettino della Unione Matematica Italiana. Serie VIII. Sezione B. Articoli di Ricerca Matematica, 5 (2002), 55-78. [98] J. Santos and P. de Oliveira, A converging finite volume scheme for hyperbolic conservation laws with source terms, Journal of Computational and Applied Mathematics, 111 (1999), 239-251. [99] H.J. Schroll and R. Winther, Finite-difference Schemes for scalar conservation laws with source terms, IMA Journal of Numerical Analysis, 16 (1996), 201-215. [100] J. Shen, C.-W. Shu, M. Zhang, A high order WENO scheme for a hierarchical size-structured population, Journal of Scientific Computing, 33 (2007), 279-291. [101] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing, 9 (1988), 1073-1084. 167 [102] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics, 77 (1988), 439-471. [103] P.V. Slingerland and J.K. Ryan and C. Vuik Position-dependent smooth- increasing accuracy-conserving (SIAC) filtering for improving discontinuous Galerkin solutions, SIAM Journal on Scientific Computing, 33 (2011), 802-825. [104] M. Spaans and J. Silk, Pregalactic black hole formation with an atomic hydro- gen equation of state, The Astrophysical Journal, 652 (2006), 902-906. [105] L. Spitzer and J. Greenstein, Continuous emission from planetary nebulae, The Astrophysical Journal, 114 (1951), 407-420. [106] G. Stratta, et. al., Dust properties at z=6.3 in the host galaxy of GRB 050904, The Astrophysical Journal, 661 (2007), L9-L12. [107] A. Tasitsiomi, Lyα radiative transfer in cosmological simulations and applica- tion to a z 8 Lyα emitter, The Astrophysical Journal, 645 (2006), 792-813. [108] W. Unno, Theoretical line contour of the Lyα radiation of ionized helium and the excitation of Bowen lines in planetary nebulae, Publications of the astro- nomical society of Japan, 7 (1955), 81-103. [109] A. Verhamme, 3D Lyα radiation transfer. III. constraints on gas and stellar properties of z 3 Lyα break galaxies (LBG) and implications for high-z LBGs and Lyα emmitters, The Astrophysical Journal, 491 (2008), 89-111. [110] A. Verhamme, D. Schaerer and A. Maselli, 3D Lyα radiation transfer. I. Under- standing Lyα line profile morphologies, The Astrophysical Journal, 460 (2006), 397-413. [111] D. Walfisch, J.K. Ryan, R.M. Kirby and R. Haimes, One-sided smoothness- increasing accuracy-conserving filtering for enhanced streamline integration 168 through discontinuous fields, Journal of Scientific Computing, 38 (2009), 164- 184. [112] J.C. Weingartner and B.T. Draine, Dust grain-size distribution and extinction in the milky way, large magellanic cloud, and small magellaanic cloud, The Astrophysical Journal, 548 (2001), 296-309. [113] S.A. Wouthuysen, On the excitation mechanism of the 21-cm (radio-frequency) interstellar hydrogen emission line, The Astrophysical Journal, 57 (1952), 31-32. [114] W. Xu and X.-P. Wu, On the formation of Lyα emission from resonantly scattered continnum photons of gamma-ray burst’s afterglow, The Astrophysical Journal, 710 (2010), 1432-1443. [115] W. Xu and X.-P. Wu and L.-Z. Fang, Time-dependent behaviour of Lyα photon transfer in a high-redshift optically thick medium, Monthly Notice of the Royal Astronomical Society, 418 (2011), 853-862. [116] M. Zhang and C.-W. Shu, An analysis of and a comparison between the discon- tinuous Galerkin and the spectral finite volume methods, Computers and Fluids, 34 (2005), 581-592. [117] Q. Zhang and C.-W. Shu, Error estimates for the third order explicit Runge- Kutta discontinuous Galerkin method for linear hyperbolic equation in one- dimension with discontinuous initial data, submitted to Numerische Mathe- matik. [118] 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 (2004), 641-666. 169 [119] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, Journal of Computational Physics, 229 (2010), 3091-3120. [120] X. Zhang and C.-W. Shu, On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, Jour- nal of Computational Physics, 229 (2010), 8918-8934. [121] X. Zhang and C.-W. Shu, Maximum-principle-satisfying and positivity- preserving high order schemes for conservation laws: Survey and new devel- opments, Proceedings of the Royal Society A, 467 (2011), 2752-2776. [122] Z. Zheng and J. Miralda-Escude, Monte Carlo simulation of Lyα scattering and application to damped Lyα systems, The Astrophysical Journal, 578 (2002), 33-42. [123] 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, 200 (2011), 2814-2827.