Introduction to Modern Scientific Programming and Numerical Methods

Posted on September 6th, 2020
Previous Article :: Next Article

We all know that the ability to use computers to solve mathematical relationships is a fundamental skill for anyone planning a career in science or engineering. For this reason, numerical analysis is part of the core curriculum in just about every undergraduate engineering department. However, my observation has been that under most of these curricula, students (at least outside the computer science department) are introduced only to numerical methods and programming under the MATLAB® environment. While MATLAB® is a great tool for algorithm prototyping, being an interpreted language, it is simply too slow for large-scale computation. This is especially true for problems that cannot be easily cast into a matrix form, such as those encountered in stochastic particle simulations. There are few opportunities for students to gain exposure to other commonly used languages, such as C/C++ (and its derivatives Java, Javascript and CUDA), Python, or even FORTRAN, which is still used by many legacy codes. By focusing solely on MATLAB®, students also fail to learn about basic programming paradigms such as variable types and dynamic memory allocation. These courses also do not cover hardware technologies, including computer clusters, graphic cards, microcontrollers, or field programmable gate arrays. For most non-CS major engineering students or professionals doing computational research, practical programming is a self-taught process.

It was essentially these skills that I tried to expose students to in my ASTE-499 Applied Scientific Computing course at USC. That course was taught mainly off slides, but in the future, students taking this course – and similar courses at different universities – will have a textbook at their disposal. We recently got a contract from CRC Press (the publisher of my Plasma Simulations by Example text) for Introduction to Modern Scientific Programming and Numerical Methods, due in July 2021! This book is co-authored with Prof. Joseph Wang at USC, and Dr. Robert Martin at the Air Force Research Laboratory (AFRL). I put this page together to solicit feedback. Below you will find the proposed table of contents. Some material will likely end up getting re-ordered but it summarizes the topics we plan to include. Please leave a comment below (or just email me) with any suggestions. Also let us know if you would like to review draft copies.

Table of Contents

1. Introduction to Scientific Computing

This chapter introduces basic concepts from scientific computing by demonstrating how to simulate a ball dropped from some height and falling under the action of gravity (and drag) and bouncing on surface impact (example differential equation).
1.1 Problem Description
We start by describing the problem. Here we introduce limits, Taylor series, differentiation, and numerical round off errors.
1.2 Pseudo code
Next, a pseudo code is developed. Basic programming concepts, such as variables, data types, for loops, if statements, functions, and input/output are introduced. We also discuss binary representation and finite-precision computer arithmetic.
1.3 Programming Languages
Various programming languages, including Assembly, Basic, Pascal, Fortran, Matlab, Python, C/C++, Java, Javascript, and R are discussed. We demonstrate how to implement a basic version of the simulation code in each language and discuss the pros and cons. We learn about interpreted vs. compiled languaged.
1.4 Integration Methods
Numerical results are compared to theoretical model. We also discuss visualization options, including use of Python/Matlab plotting and the use of Tecplot and Paraview. We discuss convergence, consistency, and alternate integration methods including forward and backward Euler, Leapfrog, Crank-Nicolson, multistep, Runge Kutta. These examples are developed primarily in Matlab and Python.
1.5 Additional Physics
Additional physics is included to increase complexity. We mainly focus on adding aerodynamic drag and surface reflection. Root finding is introduced.
1.6 Modularity and Data Reduction
The simulation is expanded to include multiple bouncing balls. This section allows us to introduce arrays, functions, and function arguments. We also discuss pass-by-reference and pass-by-value. Additional visualization options such as scatter plots and animations are discussed.

2. Numerical Analysis and Linear Algebra

In this chapter, we cover many elementary topics from Numerical Analysis, including equation solving and linear algebra.
2.1 Root Finding
We start by discussing root finding (solving equations of one variable). We cover Bisection Method, Fixed point iteration, Newton’s method. Error analysis and rate of convergence is also discussed.
2.2 Differentiation and Integration
We discuss Taylor Series in more detail and discuss “big O” notation. We then cover numerical integration (quadrature) including Trapezoidal and Simpson rule, and Gaussian Quadrature. Double (or triple) integrals are also covered.
2.3 Matrix Algebra
We switch to Linear Algebra, and introduce the concept of a matrix and a vector. We learn about matrix multiplication, and see how to implement it in Matlab and also languages where it is not natively supported (such as C++). Einstein (index) notation is also discussed. We also learn about different matrix types, such as banded, sparse, singular, symmetric, diagonally dominant, and positive definite. We demonstrate how to use rotation and translation matrix for coordinate transformation.
2.4 Solving Systems of Equations
We introduce matrix inverse for solving systems of equations. We also cover topics such as determinants and cofactors and see how to compute inverse of small matrices.
2.5 Direct Methods
We cover Gaussian elimination, a direct solution method for larger systems for which computing a matrix inverse is not practical. We discuss the Thomas algorithm for tridiagonal systems. These examples are worked out in Matlab, Python and C++.
2.6 Iterative Methods
We discuss iterative methods for solving matrix problems, including Jacobi, and Gauss-Seidel, and Successive Over Relaxation. We learn about residue norm and convergence.
2.7 Newton-Raphson Method
Non-linear systems of equations are introduced and we discuss the popular Newton-Raphson method
2.8 Conjugate Gradient
The preconditioned conjugate gradient method is discussed and convergence is compared to Jacobi
2.9 Multigrid
The multigrid method is covered next. We cover high and low frequency errors.
2.10 LU and Cholesky Decomposition
LU and Cholesky matrix decomposition is discussed and illustrated with Matlab and C++ example
2.11 Eigenvalues
Eigenvalues and eigenvectors are covered. Some popular algorithms such as power method and QR decomposition are illustrated.
2.12 Numerical Libraries
We review popular numerical libraries such as NumPy/SciPy and BLAS/LAPACK and see how to use them in our codes

3. Interpolation and Signal Processing

This chapter focus on methods for fitting curves to large amounts of data and reducing data trends
3.1 Polynomial Fits
We cover Lagrange polynomials, Hermite interpolation, Bezier splines, Richardson extrapolation
3.2 Least Squares Fit
The Least Squares method is covered next
3.3 Fourier Series
Fourier Series and the FFT method is covered. We demonstrate the method using built-in functions from Matlab and Python and with a custom C++ code.
3.4 Filtering
Methods for filtering data are covered.
3.5 Probability Distributions
Probability distribution functions are covered. We learn about cumulative distribution function and stochastic (random) sampling. Earth Mover’s method is discussed for comparing distribution functions.

4. Object Oriented Programming

This chapter we introduce advanced C++ syntax and learn about object oriented programming
4.1 Data Encapsulation
Here we discuss benefits of object oriented programming, including data encapsulation, access control, storage, and operator overloading without getting into the actual implementation details
4.2 Structures and Classes
Structs and classes are introduced here. They form the foundation of object oriented programming. We learn about access control, constructors, and destructors.
4.3 Pointers, References, and Dynamic Memory Allocation
We then start with a crash course in C++. We start by covering pointers, references, and dynamic memory allocation. These topics, while not difficult, are are a source of confusion to many new students of C++.
4.4 Operator Overloading
We learn about implementing operators for manipulating our custom classes. These operators allow us to write the code in a more concise, mathematical way.
4.5 Templates
We see how templates allow us to generalize the functionality of a class to operate on different data types. We implement a “vec3” data object for storing 3D vectors.
4.6 Storage Containers
Next we discuss various storage containers and see how they are implemented inthe C++11 standard library. These include vector, linked list, map,and set. We also learn about octrees.
4.7 Lambda Functions and Algorithms
We discuss anonymous lambda functions and see how they can be used with standard library algorithms to perform manipulations such as data sorting.
4.8 Putting it Together
We close the chapter by utilizing the new skills to write a multiple-bouncy ball simulation from Chap. 1 as well as a matrix solver. We compare the performance to Python/MATLAB version.

5. Interactive Web Applications

In this chapter we discuss using Javascript to develop simulations that run in a web browser. Not only are such simulations interactive by the basic nature of web sites, every computer has a web browser installed, allowing us to develop code even on systems that may not have a more standard compiler available, such as lab computers used for data acquisition.
5.1 HTML
We start by describing the basic syntax of HTML files. We learn about basic element types such as <head>, <body>, <p>, <h?>, <a>, <ul>, <li>, <input>, and <div>.
5.2 CSS
We next learn how to change the visualization representation using styles. We learn about the rule selectors in CSS files.
5.3 Javascript
Next we learn the basic syntax of Javascript and see how to use it to dynamically modify the webpage and how to process user input
5.4 Canvas
We then introduce the <canvas> element that can be used for 2D and 3D rendering. We mainly focus on the 2D “context”. We learn how to draw shaded circles, lines, and rectangles. Discussion of webGL (for parallel 3D rendering) is left for Chapter 10.
5.5 Javascript Examples
We use the above material to develop several interactive simulations covering methods discussed previously. We start with a simulation that adds a new bouncing ball at mouse click location. We then write code for filtering and interpolating .csv data that is drag-and-dropped from the file system.

6. Documentation and Code Testing

This chapter covers additional programming topics which are not commonly discussed but knowledge of which is critical for developing large-scale codes that can be utilized and expanded by other researchers.
6.1 Testing and Visualization
Here we start by introducing code validation and verification, and sensitivity and uncertainty analysis. We also discuss various methods for visualizing data, including isosurfaces, scatter plots, and virtual probes
6.2 Debugging
Techniques for code debugging are discussed. We also introduce the use of automated tools such as valgrind for finding memory leaks.
6.3 Coding Standards
We discuss the need for coding standards for code maintainability. Several popular standards are introduced.
6.4 Test Suites
GTest and JUnit is discussed here. We learn about benefits of automated unit testing
6.5 Version Control
We introduce version control using CVS,SVN, Mercurial, but mainly focusing on Git (current favorite).
6.6 Documentation Systems
In-line documentation using Doxygen is discussed. We include a brief description of LaTeX, focusing on mathematical formulas. We see how to include LaTeX syntax in Doxygen.

7. Partial Differential Equations

This chapter introduces basic solution approaches for selected partial differential equations. We illustrate the basic methods in Matlab and Python, but the more advanced methods are illustrated mainly using C++.
7.1 Classification
We learn about elliptic, parabolic, and hyperbollic partial differential equations. We discuss their derivation, and learn about initial and boundary conditions
7.2 Finite Difference Method
The finite different method is discussed. We learn about discretization, mesh nodes and cells, and matrix representation. First and second order differencing is discussed.
7.3 Heat Equation
We then start covering common solution methods to model PDEs. We start by discussing steady (elliptic) and un-stead (parabolic) heat (or diffusion) equation solved using Finite Difference. We write an example code in Matlab, C++, and Javascript. We cover schemes such as forward-time centered-space, RK, and Crank-Nicolson schemes. Stability and Von Neumann analysis are covered.
7.4 Finite Volume Method
The finite volume method is covered. We learn about the Divergence and Stokes Theorem, and also see how it naturally retrieves coefficients in cylindrical coordinates. We see that for uniform Cartesian meshes, the method reduces to the Finite Difference. A FVM example is developed in C++.
7.5 Finite Element Method
The FEM is discussed briefly. We use the method to develop a simple FEM solver of the heat equation in C++.
7.6 Wave Equation
We then discuss the parabolic wave equation and cover several popular solution methods.
7.7 Beam Equation
The beam equation (example of a higher-order PDE) is covered and several popular methods are illustrated.
7.8 Advection-Diffusion Equation and the SIMPLE Method
Another model equation from aerodynamics is discussed. We discuss some approaches for the Advective Term starting with the Upwind method. We then cover the SIMPLE method used for incompressible fluids.
7.9 Vlasov Equation
Vlasov Equation covering the evolution of a velocity distribution function is introduced. We cover a simple first order integration scheme.

8. Particle Methods

This chapter introduces numerical methods based on Lagrangian Mechanics and stochastic (particle) methods.
8.1 Random Numbers and Monte Carlo
We learn about using random numbers to sample values from Maxwellian velocity distribution function
8.2 Gas Kinetics
We modify the “bouncy-balls” example from Ch. 4 to simulate free molecular flow gas dynamics. We learn about using a computational mesh to sample velocity moments for computing density, stream velocity, and temperature.
8.3 Direct Simulation Monte Carlo
The Direct Simulation Monte Carlo (DSMC) method is introduced next to add intermolecular collisions
8.4 Particle in Cell
We cover Poisson’s equation (which has the same form as steady heat equation covered in Ch. 7) and electrostatics. We see how to use the PIC method to simulate ionized gas (plasma) dynamics.
8.5 Molecular Dynamics
We discuss the molecular dynamics (MD) method that is used in computational chemistry to model molecular structures
8.6 Smoothed Particle Hydrodynamics
In this section we describe the SPH method that blends Eulerian and Lagrangian description of fluid motion

9. Parallel Processing

This chapter introduces the reader to code optimization and high performance computing. These examples are worked out in C++ and CUDA
9.1 Profiling and Memory Access
We learn how to use timers and profilers to find parts of the code that run slowly. We also demonstrate how re-arranging memory access can lead to drastic performance gain, and learn about CPU cache.
9.2 Multithreading
We discuss vector operations and see how to split them up using threads. We cover race condition, atomics, and locks (and how to avoid needing them). We see how some library functions are not thread safe and lead to serialization and slow performance (i.e. Random). We also demonstrate that there is finite penalty for thread creation that may make multithreading not be efficient for small computations. Parallel efficiency is compared.
9.3 Message Passing Interface
We learn about MPI and discuss several popular matrix solution methods such as red-black ordering and domain decomposition. We also discuss approaches for particle simulations. We briefly describe setting up a cluster using Amazon Web Services and provide a short list of basic Linux commands.
9.4 CUDA
Graphics card processing using CUDA is discussed next. We cover memory transfer, pinned memory, and streams.
9.5 webGL
Returning to interactive Javascript, we see how to use webGL to add GPU processing to web applications

10. Optimization and Machine Learning

This chapter introduces optimization methods for large problems and machine learning
10.1 Optimization Basics
We discuss the need for optimization – often we do not have the full parameter space for the simulation and need to perform parameter space search. We cover the brute force and shooting methods. We see this approach is inefficient for a large parameter space.
10.2 Gradient Descent
The gradient descent method is described next. This follows up on the Conjugate Gradient discussion for matrix solving.
10.3 Genetic Algorithms
Genetic algorithms are described next. We develop simulation to find the shortest path between points.
10.4 Neural Networks
Deep neural networks, which form the foundation of machine learning, are covered. We cover basic activation functions such as ReLU and logistic curve. Cost function is also covered.
10.5 Back-propagation
We next see how back-propagation is used to update activation function weights

11. Embedded Systems

This chapter discusses programming microcontrollers and FPGAs, as they allow us to develop data acquisition sensors for code validation.
11.1 Introduction to Arduino
We learn about the basics of the popular Arduino microcontroller platform, which is programmed using a simplified form of C++.
11.2 Sensors and Electronics
Here we cover basic electronic components including resistors, capacitors, transistors, diodes, and LEDs. We also discuss breadboards and soldering.
11.3 Basic DAQ system
Putting it all together, we learn how to create a simple data acquisition system that collects pressure and humidity data from multiple sensors and stores the data onto an SD card
11.4 Edge Computing
We learn how to use TensorFlow Lite to add machine learning to microcontrollers to enable “edge computing”
11.5 FPGAs and Verilog
Here we cover what FPGAs are and introduce the syntax of Verilog. We learn how to write our own hardware instructions for performing rapid math operations. We write code for rapidly “mux”-ing data passed by an Arduino.

Appendix: Programming Syntax

This appendix summarizes common programming syntax such as variable and function declaration, flow control, memory allocation, random number sampling, and input output in Matlab, Python, Fortran, C++, R, Java, and Javascript

Example Codes

Various code examples are posted here:

14 comments to “Introduction to Modern Scientific Programming and Numerical Methods”

  1. Philip Mainwaring
    September 6, 2020 at 8:38 pm

    This looks great. As a former CS guy in aero, this is pretty much everything I would hope to see. A few random suggestions: validation cases for testing consistency, versioning (i.e. bake in to your results what build produced the result), C++11 style pointers, build systems, file formats like json & xml for configuration (please do not invent new text file formats for your solver, it makes me cry…). I would say also that some discussion of UI options wouldn’t be amiss. PyQT, etc. Maybe interprocess comm too. There’s a lot of merit to doing solvers in one language and your problem setup and viz in another (C++ and python being a common example). I also wouldn’t mind seeing some EM-specific examples, and maybe some inside baseball stuff like using ghost cells, specific treatment of BCs, etc to help improve results. Though this may be in your other book!

    Anyhow these are all relatively minor points. Looking forward to this one!

    • September 7, 2020 at 10:54 am

      Thanks Philip for many great suggestions! I plan to have a description of various common file formats in the text, somewhere. However, utilizing either is not trivial in C++, at least not without using external libraries. At least Java has an XML parser (as used by my 2D code Starfish), but I ended up writing my own JSON-like parser for CTSP.

  2. KEN DANIEL
    September 15, 2020 at 8:20 pm

    Hi Lubos,
    Here are my wild thoughts – no need to respond. Good to hear you are doing well.
    Looks like you didn’t leave anything out. How many thousands of pages are you planning on? I think you need two books. One for programming principles for engineers/scientists that includes all the computer stuff and basic numerical analysis techniques that are the most commonly used in engineering. The other is numerical analysis book with all the other numerical techniques you mention using a single programming language. The second book is similar to a few already on the market.
    I’m not sure where the chapter on embedded systems fits.
    You may want to spend some time on “make” and tool chains. Coming from matlab linking will be a foreign concept.
    I’m not sure about the logic of putting optimization in the same chapter as machine learning. I think you are talking about regression machine learning which fits more with the section on interpolation. Also, I’m not sure how interpolation fits with signal processing.
    You may want to add an example of using a large off-the-shelf code like Sparta / Paraview?
    Is the audience for the book people who need to understand what is going on inside of a code like Sparta or people who can build and augment large complex codes for a particular domain?
    Cheers,
    Ken

    • September 15, 2020 at 9:00 pm

      Thanks Ken. So we actually went back and forth on whether to do a single book or a multi-volume, and we settled on a single volume as otherwise, there is not much difference from having multiple books on related subjects. The total number of pages is going to be only about 450 (so about 40 pages per chapter). I realize this may sound too little, but the point of this text is to introduce the reader to these various topics (“I didn’t even know what I didn’t know”) and refer them to other texts for the actual details. I do plan on having a brief section on data visualization with Paraview. For your two other questions, I guess I don’t see the topics to be that incompatible. Interpolation (or perhaps what I meant is smoothing/fitting) involves looking for trends in noisy data which is essentially signal processing, no? And with optimization and machine learning, the idea is that we are looking for some parameter – or a set of relationships – that will get us a desired answer given some set of inputs, with us letting the computer find it on its own. But perhaps there is a deeper philosophical distinction between the two I am not aware of.

  3. Varun Saxena
    September 15, 2020 at 10:23 pm

    When will the book be out for sale? Or it is and I missed reading that somewhere

    • September 16, 2020 at 8:10 am

      Not until late next year. It is due to the publisher in July 2021, so likely targeting November/December release.

  4. Francisco
    October 2, 2020 at 1:41 pm

    Hi, Chemical Engineer self taught programmer here. My opinion is similar to the one given by Ken Daniel. My biggest concern is that you try to cover too many topics, that could leave you in some cases with a single page to “cover” a topic.
    But everyone will clear out once the printed version is ready, I hope to see how you manage to solve this problem.
    Best regards.

    • October 3, 2020 at 11:53 am

      Thank you Francisco. Yes, this is definitely a concern. This is why I am not trying the book to be fully inclusive of all methods/etc, but to serve as the proper introductory material, with references to more detailed texts. For example, there are 500+ page books covering CFD, and here we will try to just cover basics of Eulerian methods in about 40 pages.

  5. Greg
    October 2, 2020 at 3:35 pm

    I feel like Newton-Raphson is a bit out of order… If I were to structure chapter 2, I might have something like:

    a) Solving Linear Systems of equations
    1) Direct methods
    2) Stationary methods
    3) Steepest descent and conjugate gradient methods
    4) Krylov methods
    5) Geometric MultiGrid
    6) Algebraic Multigrid
    b) Solving Nonlinear Systems of Equations
    1) Newton-Raphson
    2) Newton-Krylov
    3) Nonlinear conjugate gradient (I know you have this in optimization, but could fit here as well).
    4) (others as relevant)

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    Also, some other things to consider might be automatic differentiation (aka algorithmic differentiation), complex-step, symbolic calculations, scientific file formats like HDF5 & NetCDF, Verification & Validation (at least what they are, and how they’re different).

    Perhaps also some discussion on how to present results / data, e.g. when to use log scales, when to use a diverging colormap vs when to use a “single-sided” colormap (e.g. white-to-red), when to use a contour plot vs a surface plot, etc.

    Would be good to present Amdahl’s Law in the parallel section.

    • October 3, 2020 at 11:51 am

      Thanks Greg, these are all good suggestions. But I think that Krylov methods are much more complicated to implement than N-R, no? For instance, I have implemented PCG many times, but every time it was just following a recipe in a numerical book. I don’t actually know how to derive the solver (yet?), but N-R is fairly straightforward.

  6. Greg
    October 4, 2020 at 9:20 pm

    Lubos,

    I think for the most part, the majority of users of scientific computing will be calling a library that already exists for the more complicated methods. I think what’s more important for these more involved solvers is to give an introduction to *how* they work, *when* to use them, and *why* they fail. To your earlier point, it seems like this book is a bit of an introductory book – with references to the vast 500+ page treatises on the subject. Thus I think it’s more important to give the reader a full overview of linear solvers THEN give a full overview of nonlinear solvers, rather than jump between linear -> nonlinear -> linear -> nonlinear (optimization). Especially since Newton methods are effectively just wrappers around linear solvers.

    Might also be a good idea to talk preconditioning…

    • October 5, 2020 at 8:41 am

      Sounds good. I will send you the draft copy of this chapter once ready, if interested, to get additional feedback.

  7. Greg
    October 8, 2020 at 12:01 pm

    Lubos,

    I’d be happy to review a draft of the chapter. Feel free to direct email me (I presume you have my email) if you’d like to.

  8. Vignesh Ramakrishnan
    February 8, 2021 at 11:06 am

    Please mention semi-discretisation schemes for PDEs(Heat PDE, Wave PDE) along with the state Laplacian matrix generation for Cartesian, semi-Cartesian and polar coordinates evaluated with examples involving various initial and boundary conditions.

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre lang="" line="" escaped="" cssfile=""> In addition, you can use \( ...\) to include equations.