$$ \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\acc}{\boldsymbol{a}} \newcommand{\rpos}{\boldsymbol{r}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\q}{\boldsymbol{q}} \newcommand{\g}{\boldsymbol{g}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\stepzero}{*} \newcommand{\stephalf}{***} \newcommand{\stepone}{**} \newcommand{\baspsi}{\psi} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integer}{\mathbb{Z}} $$

 

 

 

Table of contents

Preface
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
            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
            General mechanical vibrating system
            A sliding mass attached to a spring
            A jumping washing machine
            Motion of a pendulum
            Dynamic free body diagram during pendulum motion
            Motion of an elastic pendulum
            Vehicle on a bumpy road
            Bouncing ball
            Two-body gravitational problem
            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
            Exercise 1.28: Implement the PEFRL algorithm
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
      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
            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
            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
            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
      Implementation
            Scalar computations
            Vectorized computations
            Verification
            Visualization
      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
            Elastic waves in a rod
            Waves on a membrane
            The acoustic model for seismic waves
            Sound waves in liquids and gases
            Spherical waves
            The linear shallow water equations
            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
      Linearization at the differential equation level
            Explicit time integration
            Backward Euler scheme and Picard iteration
            Backward Euler scheme and Newton's method
            Crank-Nicolson discretization
      1D stationary nonlinear differential equations
            Finite difference discretization
            Solution of algebraic equations
      Multi-dimensional nonlinear PDE problems
            Finite difference discretization
            Continuation methods
      Operator splitting methods
            Ordinary operator splitting for ODEs
            Strang splitting for ODEs
            Example: Logistic growth
            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
      Finite differences of \( t^n \)
            Software
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
            Model with damping and nonlinearity
            Extension to quadratic damping
            The general model formulated as first-order ODEs
      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
            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