Finite Difference Computing with PDEs

Contents:

• Finite Difference Computing with PDEs
• Preface
  • Why finite differences?
  • Simplify, understand, generalize
  • Constructive mathematics
  • All nuts and bolts
  • Python as programming language
  • Program verification
  • Vectorized code
  • Analysis via exact solutions of discrete equations
  • Code-inspired mathematical notation
  • Limited scope
  • Focus on wave phenomena
  • Independent chapters
  • Supplementary materials
  • Acknowledgments
• Vibration ODEs
  • Finite difference discretization
    • A basic model for vibrations
    • A centered finite difference scheme
  • Implementation
    • Making a solver function
    • Verification
    • Scaled model
  • Visualization of long time simulations
    • Using a moving plot window
    • Making animations
    • Using Bokeh to compare graphs
    • Using a line-by-line ascii plotter
    • Empirical analysis of the solution
  • Analysis of the numerical scheme
    • Deriving a solution of the numerical scheme
    • The error in the numerical frequency
    • Empirical convergence rates and adjusted \(\omega\)
    • Exact discrete solution
    • Convergence
    • The global error
    • Stability
    • About the accuracy at the stability limit
  • Alternative schemes based on 1st-order equations
    • The Forward Euler scheme
    • The Backward Euler scheme
    • The Crank-Nicolson scheme
    • Comparison of schemes
    • Runge-Kutta methods
    • Analysis of the Forward Euler scheme
  • Energy considerations
    • Derivation of the energy expression
    • An error measure based on energy
  • The Euler-Cromer method
    • Forward-backward discretization
    • Equivalence with the scheme for the second-order ODE
    • Implementation
    • The Stoermer-Verlet algorithm
  • Staggered mesh
    • The Euler-Cromer scheme on a staggered mesh
    • Implementation of the scheme on a staggered mesh
  • Exercises and Problems
    • Problem 1.1: Use linear/quadratic functions for verification
    • Exercise 1.2: Show linear growth of the phase with time
    • Exercise 1.3: Improve the accuracy by adjusting the frequency
    • Exercise 1.4: See if adaptive methods improve the phase error
    • Exercise 1.5: Use a Taylor polynomial to compute \(u^1\)
    • Problem 1.6: Derive and investigate the velocity Verlet method
    • Problem 1.7: Find the minimal resolution of an oscillatory function
    • Exercise 1.8: Visualize the accuracy of finite differences for a cosine function
    • Exercise 1.9: Verify convergence rates of the error in energy
    • Exercise 1.10: Use linear/quadratic functions for verification
    • Exercise 1.11: Use an exact discrete solution for verification
    • Exercise 1.12: Use analytical solution for convergence rate tests
    • Exercise 1.13: Investigate the amplitude errors of many solvers
    • Problem 1.14: Minimize memory usage of a simple vibration solver
    • Problem 1.15: Minimize memory usage of a general vibration solver
    • Exercise 1.16: Implement the Euler-Cromer scheme for the generalized model
    • Problem 1.17: Interpret \([D_tD_t u]^n\) as a forward-backward difference
    • Exercise 1.18: Analysis of the Euler-Cromer scheme
• Generalization: damping, nonlinearities, and excitation
  • A centered scheme for linear damping
  • A centered scheme for quadratic damping
  • A forward-backward discretization of the quadratic damping term
  • Implementation
  • Verification
    • Constant solution
    • Linear solution
    • Quadratic solution
    • Catching bugs
  • Visualization
  • User interface
  • The Euler-Cromer scheme for the generalized model
  • The Stoermer-Verlet algorithm for the generalized model
  • A staggered Euler-Cromer scheme for a generalized model
  • The PEFRL 4th-order accurate algorithm
• Exercises and Problems
  • Exercise 1.19: Implement the solver via classes
  • Problem 1.20: Use a backward difference for the damping term
  • Exercise 1.21: Use the forward-backward scheme with quadratic damping
• Applications of vibration models
  • Oscillating mass attached to a spring
    • Scaling
    • The physics
  • General mechanical vibrating system
    • Scaling
  • A sliding mass attached to a spring
  • A jumping washing machine
  • Motion of a pendulum
    • Simple pendulum
    • Physical pendulum
  • Dynamic free body diagram during pendulum motion
    • Writing the solver
    • Drawing the free body diagram
    • Making the animated free body diagram
  • Motion of an elastic pendulum
    • Remarks about an elastic vs a non-elastic pendulum
    • Initial conditions
    • The complete ODE problem
    • Scaling
    • Remark on the non-elastic limit
  • Vehicle on a bumpy road
  • Bouncing ball
  • Two-body gravitational problem
    • The governing equations
    • Scaling
    • Solution in a special case: planet orbiting a star
  • Electric circuits
• Exercises
  • Exercise 1.22: Simulate resonance
  • Exercise 1.23: Simulate oscillations of a sliding box
  • Exercise 1.24: Simulate a bouncing ball
  • Exercise 1.25: Simulate a simple pendulum
  • Exercise 1.26: Simulate an elastic pendulum
  • Exercise 1.27: Simulate an elastic pendulum with air resistance
    • Remarks
  • Exercise 1.28: Implement the PEFRL algorithm
    • Remarks
• Wave equations
  • Simulation of waves on a string
    • Discretizing the domain
    • The discrete solution
    • Fulfilling the equation at the mesh points
    • Replacing derivatives by finite differences
    • Formulating a recursive algorithm
    • Sketch of an implementation
  • Verification
    • A slightly generalized model problem
    • Using an analytical solution of physical significance
    • Manufactured solution and estimation of convergence rates
    • Constructing an exact solution of the discrete equations
  • Implementation
    • Callback function for user-specific actions
    • The solver function
    • Verification: exact quadratic solution
    • Verification: convergence rates
    • Visualization: animating the solution
    • Running a case
    • Working with a scaled PDE model
  • Vectorization
    • Operations on slices of arrays
    • Finite difference schemes expressed as slices
    • Verification
    • Efficiency measurements
    • Remark on the updating of arrays
  • Exercises
    • Exercise 2.1: Simulate a standing wave
    • Exercise 2.2: Add storage of solution in a user action function
    • Exercise 2.3: Use a class for the user action function
    • Exercise 2.4: Compare several Courant numbers in one movie
    • Exercise 2.5: Implementing the solver function as a generator
    • Project 2.6: Calculus with 1D mesh functions
• Generalization: reflecting boundaries
  • Neumann boundary condition
  • Discretization of derivatives at the boundary
  • Implementation of Neumann conditions
  • Index set notation
  • Verifying the implementation of Neumann conditions
  • Alternative implementation via ghost cells
    • Idea
    • Implementation
• Generalization: variable wave velocity
  • The model PDE with a variable coefficient
  • Discretizing the variable coefficient
  • Computing the coefficient between mesh points
  • How a variable coefficient affects the stability
  • Neumann condition and a variable coefficient
  • Implementation of variable coefficients
  • A more general PDE model with variable coefficients
  • Generalization: damping
• Building a general 1D wave equation solver
  • User action function as a class
    • The code
    • Dissection
  • Pulse propagation in two media
• Exercises
  • Exercise 2.7: Find the analytical solution to a damped wave equation
  • Problem 2.8: Explore symmetry boundary conditions
  • Exercise 2.9: Send pulse waves through a layered medium
  • Exercise 2.10: Explain why numerical noise occurs
  • Exercise 2.11: Investigate harmonic averaging in a 1D model
  • Problem 2.12: Implement open boundary conditions
    • Remarks
  • Exercise 2.13: Implement periodic boundary conditions
  • Exercise 2.14: Compare discretizations of a Neumann condition
  • Exercise 2.15: Verification by a cubic polynomial in space
• Analysis of the difference equations
  • Properties of the solution of the wave equation
  • More precise definition of Fourier representations
  • Stability
    • Preliminary results
    • Numerical wave propagation
  • Numerical dispersion relation
  • Extending the analysis to 2D and 3D
• Finite difference methods for 2D and 3D wave equations
  • Multi-dimensional wave equations
  • Mesh
  • Discretization
    • Discretizing the PDEs
    • Handling boundary conditions where \(u\) is known
    • Discretizing the Neumann condition
• Implementation
  • Scalar computations
    • Domain and mesh
    • Solution arrays
    • Index sets
    • Computing the solution
  • Vectorized computations
  • Verification
    • Testing a quadratic solution
  • Visualization
    • Matplotlib
    • Gnuplot
    • Mayavi
• Exercises
  • Exercise 2.16: Check that a solution fulfills the discrete model
  • Project 2.17: Calculus with 2D mesh functions
  • Exercise 2.18: Implement Neumann conditions in 2D
  • Exercise 2.19: Test the efficiency of compiled loops in 3D
• Applications of wave equations
  • Waves on a string
    • Damping
    • External forcing
    • Modeling the tension via springs
  • Elastic waves in a rod
  • Waves on a membrane
  • The acoustic model for seismic waves
    • Anisotropy
  • Sound waves in liquids and gases
  • Spherical waves
  • The linear shallow water equations
    • Wind drag on the surface
    • Bottom drag
    • Effect of the Earth’s rotation
  • Waves in blood vessels
  • Electromagnetic waves
• Exercises
  • Exercise 2.20: Simulate waves on a non-homogeneous string
  • Exercise 2.21: Simulate damped waves on a string
  • Exercise 2.22: Simulate elastic waves in a rod
  • Exercise 2.23: Simulate spherical waves
  • Problem 2.24: Earthquake-generated tsunami over a subsea hill
  • Problem 2.25: Earthquake-generated tsunami over a 3D hill
  • Problem 2.26: Investigate Mayavi for visualization
  • Problem 2.27: Investigate visualization packages
  • Problem 2.28: Implement loops in compiled languages
  • Exercise 2.29: Simulate seismic waves in 2D
  • Project 2.30: Model 3D acoustic waves in a room
  • Project 2.31: Solve a 1D transport equation
  • Problem 2.32: General analytical solution of a 1D damped wave equation
  • Problem 2.33: General analytical solution of a 2D damped wave equation
• Diffusion equations
  • An explicit method for the 1D diffusion equation
    • The initial-boundary value problem for 1D diffusion
    • Forward Euler scheme
    • Implementation
    • Verification
    • Numerical experiments
  • Implicit methods for the 1D diffusion equation
    • Backward Euler scheme
    • Sparse matrix implementation
    • Crank-Nicolson scheme
    • The unifying \(\theta\) rule
    • Experiments
    • The Laplace and Poisson equation
  • Analysis of schemes for the diffusion equation
    • Properties of the solution
    • Analysis of discrete equations
    • Analysis of the finite difference schemes
    • Analysis of the Forward Euler scheme
    • Analysis of the Backward Euler scheme
    • Analysis of the Crank-Nicolson scheme
    • Analysis of the Leapfrog scheme
    • Summary of accuracy of amplification factors
    • Analysis of the 2D diffusion equation
    • Explanation of numerical artifacts
  • Exercises
    • Exercise 3.1: Explore symmetry in a 1D problem
    • Exercise 3.2: Investigate approximation errors from a \(u_x=0\) boundary condition
    • Exercise 3.3: Experiment with open boundary conditions in 1D
    • Exercise 3.4: Simulate a diffused Gaussian peak in 2D/3D
    • Exercise 3.5: Examine stability of a diffusion model with a source term
  • Diffusion in heterogeneous media
    • Discretization
    • Implementation
    • Stationary solution
    • Piecewise constant medium
    • Implementation of diffusion in a piecewise constant medium
    • Axi-symmetric diffusion
    • Spherically-symmetric diffusion
  • Diffusion in 2D
    • Discretization
    • Numbering of mesh points versus equations and unknowns
    • Algorithm for setting up the coefficient matrix
    • Implementation with a dense coefficient matrix
    • Verification: exact numerical solution
    • Verification: convergence rates
    • Implementation with a sparse coefficient matrix
    • The Jacobi iterative method
    • Implementation of the Jacobi method
    • Test problem: diffusion of a sine hill
    • The relaxed Jacobi method and its relation to the Forward Euler method
    • The Gauss-Seidel and SOR methods
    • Scalar implementation of the SOR method
    • Vectorized implementation of the SOR method
    • Direct versus iterative methods
    • The Conjugate gradient method
    • What is the recommended method for solving linear systems?
  • Random walk
    • Random walk in 1D
    • Statistical considerations
    • Playing around with some code
    • Equivalence with diffusion
    • Implementation of multiple walks
    • Demonstration of multiple walks
    • Ascii visualization of 1D random walk
    • Random walk as a stochastic equation
    • Random walk in 2D
    • Random walk in any number of space dimensions
    • Multiple random walks in any number of space dimensions
  • Applications
    • Diffusion of a substance
    • Heat conduction
    • Porous media flow
    • Potential fluid flow
    • Streamlines for 2D fluid flow
    • The potential of an electric field
    • Development of flow between two flat plates
    • Flow in a straight tube
    • Tribology: thin film fluid flow
    • Propagation of electrical signals in the brain
  • Exercises
    • Exercise 3.6: Stabilizing the Crank-Nicolson method by Rannacher time stepping
    • Project 3.7: Energy estimates for diffusion problems
    • Exercise 3.8: Splitting methods and preconditioning
    • Problem 3.9: Oscillating surface temperature of the earth
    • Problem 3.10: Oscillating and pulsating flow in tubes
    • Problem 3.11: Scaling a welding problem
    • Exercise 3.12: Implement a Forward Euler scheme for axi-symmetric diffusion
• Advection-dominated equations
  • One-dimensional time-dependent advection equations
    • Simplest scheme: forward in time, centered in space
    • Analysis of the scheme
    • Leapfrog in time, centered differences in space
    • Upwind differences in space
    • Periodic boundary conditions
    • Implementation
    • A Crank-Nicolson discretization in time and centered differences in space
    • The Lax-Wendroff method
    • Analysis of dispersion relations
  • One-dimensional stationary advection-diffusion equation
    • A simple model problem
    • A centered finite difference scheme
    • Remedy: upwind finite difference scheme
  • Time-dependent convection-diffusion equations
    • Forward in time, centered in space scheme
    • Forward in time, upwind in space scheme
  • Two-dimensional advection-diffusion equations
  • Applications of advection equations
    • Transport of a substance
    • Transport of a heat
  • Exercises
    • Exercise 4.1: Analyze 1D stationary convection-diffusion problem
    • Exercise 4.2: Interpret upwind difference as artificial diffusion
• Nonlinear problems
  • Introduction of basic concepts
    • Linear versus nonlinear equations
    • A simple model problem
    • Linearization by explicit time discretization
    • Exact solution of nonlinear algebraic equations
    • Linearization
    • Picard iteration
    • Linearization by a geometric mean
    • Newton’s method
    • Relaxation
    • Implementation and experiments
    • Generalization to a general nonlinear ODE
    • Systems of ODEs
• Systems of nonlinear algebraic equations
  • Picard iteration
  • Newton’s method
  • Stopping criteria
  • Example: A nonlinear ODE model from epidemiology
    • Implicit time discretization
    • A Picard iteration
    • Newton’s method
• Linearization at the differential equation level
  • Explicit time integration
  • Backward Euler scheme and Picard iteration
  • Backward Euler scheme and Newton’s method
    • Linearization via Taylor expansions
    • Similarity with Picard iteration
    • Implementation
    • Derivation with alternative notation
  • Crank-Nicolson discretization
• 1D stationary nonlinear differential equations
  • Finite difference discretization
  • Solution of algebraic equations
    • The structure of the equation system
    • Picard iteration
    • Mesh with two cells
    • Newton’s method
• Multi-dimensional nonlinear PDE problems
  • Finite difference discretization
    • Picard iteration
    • Newton’s method
  • Continuation methods
• Operator splitting methods
  • Ordinary operator splitting for ODEs
  • Strang splitting for ODEs
  • Example: Logistic growth
    • Splitting techniques
    • Verbose implementation
    • Compact implementation
    • Results
  • Reaction-diffusion equation
  • Example: Reaction-Diffusion with linear reaction term
  • Analysis of the splitting method
• Exercises
  • Problem 5.1: Determine if equations are nonlinear or not
  • Problem 5.2: Derive and investigate a generalized logistic model
  • Problem 5.3: Experience the behavior of Newton’s method
  • Exercise 5.4: Compute the Jacobian of a \(2\times 2\) system
  • Problem 5.5: Solve nonlinear equations arising from a vibration ODE
  • Exercise 5.6: Find the truncation error of arithmetic mean of products
  • Problem 5.7: Newton’s method for linear problems
  • Problem 5.8: Discretize a 1D problem with a nonlinear coefficient
  • Problem 5.9: Linearize a 1D problem with a nonlinear coefficient
  • Problem 5.10: Finite differences for the 1D Bratu problem
  • Problem 5.11: Discretize a nonlinear 1D heat conduction PDE by finite differences
  • Problem 5.12: Differentiate a highly nonlinear term
  • Exercise 5.13: Crank-Nicolson for a nonlinear 3D diffusion equation
  • Problem 5.14: Find the sparsity of the Jacobian
  • Problem 5.15: Investigate a 1D problem with a continuation method
• Appendix: Useful formulas
  • Finite difference operator notation
  • Truncation errors of finite difference approximations
  • Finite differences of exponential functions
    • Complex exponentials
    • Real exponentials
  • Finite differences of \(t^n\)
• Appendix: Truncation error analysis
  • Overview of truncation error analysis
    • Abstract problem setting
    • Error measures
  • Truncation errors in finite difference formulas
    • Example: The backward difference for \(u'(t)\)
    • Example: The forward difference for \(u'(t)\)
    • Example: The central difference for \(u'(t)\)
    • Overview of leading-order error terms in finite difference formulas
    • Software for computing truncation errors
• Exponential decay ODEs
  • Forward Euler scheme
  • Crank-Nicolson scheme
  • The \(\theta\)-rule
  • Using symbolic software
  • Empirical verification of the truncation error
  • Increasing the accuracy by adding correction terms
  • Extension to variable coefficients
  • Exact solutions of the finite difference equations
  • Computing truncation errors in nonlinear problems
• Vibration ODEs
  • Linear model without damping
    • The truncation error of a centered finite difference scheme
    • The truncation error of approximating \(u'(0)\)
    • Truncation error of the equation for the first step
    • Computing correction terms
  • Model with damping and nonlinearity
  • Extension to quadratic damping
  • The general model formulated as first-order ODEs
    • A centered scheme on a staggered mesh
• Wave equations
  • Linear wave equation in 1D
  • Finding correction terms
  • Extension to variable coefficients
  • 1D wave equation on a staggered mesh
  • Linear wave equation in 2D/3D
• Diffusion equations
  • Linear diffusion equation in 1D
    • The Forward Euler scheme in time
    • The Crank-Nicolson scheme in time
  • Nonlinear diffusion equation in 1D
• Exercises
  • Exercise B.1: Truncation error of a weighted mean
  • Exercise B.2: Simulate the error of a weighted mean
  • Exercise B.3: Verify a truncation error formula
  • Problem B.4: Truncation error of the Backward Euler scheme
  • Exercise B.5: Empirical estimation of truncation errors
  • Exercise B.6: Correction term for a Backward Euler scheme
  • Problem B.7: Verify the effect of correction terms
  • Problem B.8: Truncation error of the Crank-Nicolson scheme
  • Problem B.9: Truncation error of \(u'=f(u,t)\)
  • Exercise B.10: Truncation error of \([D_t D_tu]^n\)
  • Exercise B.11: Investigate the impact of approximating \(u'(0)\)
  • Problem B.12: Investigate the accuracy of a simplified scheme
• Appendix: Software engineering; wave equation model
  • A 1D wave equation simulator
    • Mathematical model
    • Numerical discretization
    • A solver function
  • Saving large arrays in files
    • Using savez to store arrays in files
    • Using joblib to store arrays in files
    • Using a hash to create a file or directory name
  • Software for the 1D wave equation
    • Making hash strings from input data
    • Avoiding rerunning previously run cases
    • Verification
  • Programming the solver with classes
    • Class Problem
    • Class Mesh
    • Class Function
    • Class Solver
  • Migrating loops to Cython
    • Declaring variables and annotating the code
    • Visual inspection of the C translation
    • Building the extension module
    • Calling the Cython function from Python
  • Migrating loops to Fortran
    • The Fortran subroutine
    • Building the Fortran module with f2py
    • How to avoid array copying
  • Migrating loops to C via Cython
    • Translating index pairs to single indices
    • The complete C code
    • The Cython interface file
    • Building the extension module
  • Migrating loops to C via f2py
    • Migrating loops to C++ via f2py
  • Exercises
    • Exercise C.1: Explore computational efficiency of numpy.sum versus built-in sum
    • Exercise C.2: Make an improved numpy.savez function
    • Exercise C.3: Visualize the impact of the Courant number
    • Exercise C.4: Visualize the impact of the resolution
• References

Index

• Index