$$ \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\vexd}[1]{{v_{\small\mbox{e}, #1}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \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{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \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{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$ next

Introduction to computing with finite difference methods

Introduction to computing with finite difference methods

Hans Petter Langtangen [1, 2]

[1] Center for Biomedical Computing, Simula Research Laboratory
[2] Department of Informatics, University of Oslo

Dec 14, 2013

Note: PRELIMINARY VERSION

Table of contents

Finite difference methods
      A basic model for exponential decay
      The Forward Euler scheme
            Step 1: Discretizing the domain
            Step 2: Fulfilling the equation at discrete time points
            Step 3: Replacing derivatives by finite differences
            Step 4: Formulating a recursive algorithm
      The Backward Euler scheme
      The Crank-Nicolson scheme
      The unifying \( \theta \)-rule
      Constant time step
      Compact operator notation for finite differences
Implementation
      Making a solver function
            Function for computing the numerical solution
            Integer division
            Doc strings
            Formatting of numbers
            Running the program
      Verifying the implementation
            Running a few algorithmic steps by hand
            Comparison with an exact discrete solution
      Computing the numerical error as a mesh function
      Computing the norm of the numerical error
            Scalar computing
      Plotting solutions
            Plotting multiple curves
            Experiments with computing and plotting
            Combining plot files
            Plotting with SciTools
      Creating command-line interfaces
            Reading a sequence of command-line arguments
            Working with an argument parser
      Creating a graphical web user interface
            Making a compute function
            Generating the user interface
            Running the web application
      Computing convergence rates
            Estimating \( r \)
            Implementation
            Debugging via convergence rates
      Memory-saving implementation
Software engineering
      Making a module
      Prefixing imported functions by the module name
      Doctests
      Unit testing with nose
            Basic use of nose
            Alternative assert statements
            Applying nose
            Installation of nose
            Using nose to test modules with doctests
      Classical class-based unit testing
            Basic use of unittest
            Demonstration of unittest
      Implementing simple problem and solver classes
            The problem class
            The solver class
            The visualizer class
            Combining the objects
      Improving the problem and solver classes
            A generic class for parameters
            The problem class
            The solver class
            The visualizer class
Performing scientific experiments
      Software
      Combining plot files
      Interpreting output from other programs
      Making a report
            Plain HTML
            HTML with MathJax
            LaTeX
            Sphinx
            Markdown
            Wiki formats
            Doconce
            Worked example
      Publishing a complete project
Exercises and Problems
      Exercise 1: Derive schemes for Newton's law of cooling
      Exercise 2: Implement schemes for Newton's law of cooling
      Exercise 3: Find time of murder from body temperature
      Exercise 4: Experiment with integer division
      Exercise 5: Experiment with wrong computations
      Exercise 6: Plot the error function
      Exercise 7: Compare methods for a given time mesh
      Exercise 8: Change formatting of numbers and debug
      Problem 9: Write a doctest
      Problem 10: Write a nose test
      Problem 11: Make a module
      Exercise 12: Make use of a class implementation
      Exercise 13: Generalize a class implementation
      Exercise 14: Generalize an advanced class implementation
Analysis of finite difference equations
      Experimental investigation of oscillatory solutions
      Exact numerical solution
      Stability
      Comparing amplification factors
      Series expansion of amplification factors
      The fraction of numerical and exact amplification factors
      The global error at a point
      Integrated errors
      Truncation error
      Consistency, stability, and convergence
Exercises
      Exercise 15: Visualize the accuracy of finite differences \( u=e^{-at} \)
      Exercise 16: Explore the \( \theta \)-rule for exponential growth
Model extensions
      Generalization: including a variable coefficient
      Generalization: including a source term
            Schemes
      Implementation of the generalized model problem
            Deriving the \( \theta \)-rule formula
            The Python code
            Coding of variable coefficients
      Verifying a constant solution
      Verification via manufactured solutions
      Extension to systems of ODEs
General first-order ODEs
      Generic form
      The \( \theta \)-rule
      An implicit 2-step backward scheme
      Leapfrog schemes
            The ordinary Leapfrog scheme
            The filtered Leapfrog scheme
      The 2nd-order Runge-Kutta scheme
      A 2nd-order Taylor-series method
      The 2nd- and 3rd-order Adams-Bashforth schemes
      4th-order Runge-Kutta scheme
      The Odespy software
      Example: Runge-Kutta methods
            Remark about using the \( \theta \)-rule in Odespy
      Example: Adaptive Runge-Kutta methods
Exercises
      Exercise 17: Experiment with precision in tests and the size of \( u \)
      Exercise 18: Implement the 2-step backward scheme
      Exercise 19: Implement the 2nd-order Adams-Bashforth scheme
      Exercise 20: Implement the 3rd-order Adams-Bashforth scheme
      Exercise 21: Analyze explicit 2nd-order methods
      Problem 22: Implement and investigate the Leapfrog scheme
      Problem 23: Make a unified implementation of many schemes
Applications of exponential decay models
      Scaling
      Evolution of a population
      Compound interest and inflation
      Radioactive Decay
            Deterministic model
            Stochastic model
            Relation between stochastic and deterministic models
      Newton's law of cooling
      Decay of atmospheric pressure with altitude
            Multiple atmospheric layers
            Simplification: \( L=0 \)
            Simplification: one-layer model
      Compaction of sediments
      Vertical motion of a body in a viscous fluid
            Overview of forces
            Equation of motion
            Terminal velocity
            A Crank-Nicolson scheme
            Physical data
            Verification
            Scaling
      Decay ODEs from solving a PDE by Fourier expansions
Exercises and Projects
      Exercise 24: Simulate an oscillating cooling process
      Exercise 25: Radioactive decay of Carbon-14
      Exercise 26: Simulate stochastic radioactive decay
      Exercise 27: Radioactive decay of two substances
      Exercise 28: Simulate the pressure drop in the atmosphere
      Exercise 29: Make a program for vertical motion in a fluid
      Project 30: Simulate parachuting
      Exercise 31: Formulate vertical motion in the atmosphere
      Exercise 32: Simulate vertical motion in the atmosphere
      Exercise 33: Compute \( y=|x| \) by solving an ODE
      Exercise 34: Simulate growth of a fortune with random interest rate
      Exercise 35: Simulate a population in a changing environment
      Exercise 36: Simulate logistic growth
      Exercise 37: Rederive the equation for continuous compound interest
Bibliography

Finite difference methods for partial differential equations (PDEs) employ a range of concepts and tools that can be introduced and illustrated in the context of simple ordinary differential equation (ODE) examples. This is what we do in the present document. By first working with ODEs, we keep the mathematical problems to be solved as simple as possible (but no simpler), thereby allowing full focus on understanding the key concepts and tools. The choice of topics in the forthcoming treatment of ODEs is therefore solely dominated by what carries over to numerical methods for PDEs.

Theory and practice are primarily illustrated by solving the very simple ODE \( u'=-au \), \( u(0)=I \), where \( a>0 \) is a constant, but we also address the generalized problem \( u'=-a(t)u + b(t) \) and the nonlinear problem \( u'=f(u,t) \). The following topics are introduced:

The exposition in a nutshell. Everything we cover is put into a practical, hands-on context. All mathematics is translated into working computing codes, and all the mathematical theory of finite difference methods presented here is motivated from a strong need to understand strange behavior of programs. Two fundamental questions saturate the text:

next