Introduction to Modern Scientific Programming and Numerical Methods

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:

• Chapter 7 Javascript/WebGL heat equation solver: heat-webgl.html