Introduction to Finite Differences
Finite Elements in the Geosciences
Over the last several decades, computing has played an increasingly important role in the earth sciences. This trend has developed along with the rapid increase in computational power available to scientists, as well as the increasing availability of numerical codes that can be used to solve problems. However, to use these numerical codes effectively, scientists must have a basic understanding of how these codes work. This class will introduce students to basic numerical approaches to solving problems in the earth sciences. We will focus primarily on finite difference and finite element methods. Example problems will be insprired by the solid earth geopysics, but will be general enough that the understanding gained from this class can be used in broader fields of the earth sciences.
Homework: To catch up to speed on programming in Python, we will read the first several chapters of S20 and/or LL20. Please tell us which book you find more useful for this introduction.
Summary: We decided to proceed using LL20, starting with Chapter 6. The earlier chapters of LL20 and the entirety of S20 are useful references for implementing Python code, and it would be useful to read them and implement some of the examples. We do not plan to go through these chapters in detail, unless we need to learn something specific from them.
Homework: For next time, read LL20 Chapter 6: Computing Integrals and Testing Code. Carefully and write your own Python codes to implement the examples in this chapter. For a specific exercise, compute the integral of the sine function from 0 to pi (as in exercise 6.4a on page 169) using your own code to implement the trapezoidal and midpoint methods. Make your own version of the table on page 145, showing how the solution changes as you increase n. For a bonus application, try implementing the monte carlo method to solve this integral. Be prepared to share your code, and to demonstrate how it works, with the group (e.g., using share screen on Zoom).
Read Ahead: Next time, we will continue with LL20 Chapters 7 and 8. It would be a good idea to familiarize yourself with these chapters in advance of our next meeting.
Summary: We talked about methods for numerical integration (LL20 Chapter 6), using trapezoid rule and the midpoint method. We discussed the monte carlo method, and generally talked about programming techniques.
Homework: For next time, read LL20 Sections 8.1-8.3. These sections give an introduction to writing code to solve ordinary differential equations, with applications to exponential problems. Please write code to implement the problems discussed in these sections. Exercises 8.1-8.3 on pages 273-274 will help with this.
Summary: We went through LL20, sections 8.1 to 8.3 and discussed some examples from the exercises assigned.
Homework: For next time, read sections 8.4 and 8.5 of LL20 (osciallating problems), and do exercises 8.9 and 8.11.
Summary: We went through LL20, sections 8.4 to 8.5. We reviewed the differences between the Forward Euler method and more advanced finite difference methods such as Euler-Cromer (first order) or Runge-Kutta (second or even fourth order) methods. We tested our codes to see how the solution is affect by changing resolution.
Homework: For next time, read chapter 9 of LL20 (partial differential equations) and prepare a coding solution for exercise 9.2b.
Summary: We went through LL20, chapter 9. We reviewed the diffusion equation, and how it is implemented in finite-differences. We payed special attention to the boundary conditions, including how they are implemented in the finite difference solution. We discussed the temporal resolution (time step size) that is necessary to obtain a stable solution, for a given spatial resolution.
Homework: For next time, we begin the LL17 book, which goes into more detail about finite difference methods. We start with Chapter 1 (vibrational problems), sections 1.1 to 1.7. The book goes into more detail than we need about some of the methods. Thus, you can pass through the following subsections relatively quickly: 1.3.2-1.3.5, 1.4.6, 1.5.6, 1.74.
We will not assign an exercise this week, so we can focus on the methods. However, it would be good to familiarize yourself with some of the codes associated with this chapter, and try to get some of them to work. Here is the website that provides some code for this book: [LL17 codes]
Note that the relevant code for this chapter can be found using the "vib" link.
Summary: We went through LL17, chapter 1. We discussed the stability or instability of finite difference solutions to the vibrational ODE, how we can determine if a solution is stable, and the conditions that make a solution unstable. We discussed the different solution schemes (e.g., Centered Finite Difference, Forward Euler, Backward Euler, Crank-Nicolson, Runge-Kutta, Forward-Backward Euler-Cromer) and their implementation.
Homework: For next time, we continue with the LL17 book, moving into solutions to the wave equation in Chapter 2. In particular, we will discuss sections 2.1 to 2.4 in detail, but we will omit sections 2.2.3, 2.2.4, 2.3.4, and 2.4.4. To practice implementation, let's do exercises 2.1 and 2.4 (in section 2.5). Also, let's read sections 2.6.1 to 2.6.3 (about reflecting boundary conditions), section 2.7 (about variable wave velocity), and section 2.11 (about FD methods in 2D and 3D).
Summary: We went through LL17, chapter 2. We discussed using a stencil to implement a finite difference solution to the wave equation in 1 dimension (and also in 2 dimensions). The discussed the importance of the Courant condition, which specifies the relationship between spatial and temporal resolutions, and leads to an unstable solution if it is violated. We implented standing waves using a finite difference scheme, and discussed implementation of boundary conditions.
Homework: For next time, we will finish our reading of the LL17 book, looking at the diffusion equation (Chapter 3) and advection (Chapter 4). For Chapter 3, let's read sections 3.1 and 3.2, as well as 3.5.1 to 3.3.5 (heterogeneous media) and 3.6.1-3.6.2 (two dimensions). For Chapter 4, we will read sections 4.1.1 to 4.1.4.
Summary: We went through LL17, sections 3.1-3.2, 3.5, 3.6, and 4.1. In doing so we discussed the diffusion equation and in particular the difference between the explicity solution (which we used before for the wave equation) and the implicity solution. The implicity solution is useful for the diffusion equation because it allows for longer timesteps for a given spatial resolution (larger F parameter). However, the implicit solution requires inversion of a generally sparse or banded matrix, and there are efficient ways to do this. The advection equation (section 4.1) is most readily solved with a simple upwind differences scheme.
Homework: For next time, we start with finite elements using the LM19 book. Let's read chapter 1, and chapter 2 sections 2.1-2.4, as well as a quick look at section 2.5. We have a three-part exercise: (a) Write your own numerical integration program (as we did early in the semester) to reproduce the example in section 2.2.3 and plot the result (you should get something similar to figure 2.2). (b) Now apply this same numerical integration method to approximate the function y=exp(x^2) using a linear function in the span [0,2]. (c) Finally, find the analytical expressions for the 3 Lagrange polynomials of order 2 for an element with points at [0.0, 0.5, 1.0], and plot them.
Summary: We went through chapters 1 and 2 of LM19, and sections 2.1-2.4 in particular. We discussed how to approximate vectors using basis vectors (easier for orthogonal vectors) and we expanded on this idea to include approximation of functions using (preferably orthogonal) basis functions. We examined the example of the least squares method, which is related to the Taylor series expansion of polynomial functions. Better sets of basis functions (because they have a higher degree or orthogonality) are the sin and cos functions (e.g., Fourier series), as well as the Lagrange polynomials (good for interpolation problems) and Bernstein polynomials. We set up the formalism for expressing any given function as a sum of a set of basis functions, where each basis function is multiplied by a given coefficient. We discussed the method for computing these coefficients, and used the exercise in section 2.2.3 as an example.
Homework: For next time, we will continue with the LM19 book, moving into finite element basis functions in chapter 3. Let's read sections 3.1, 3.2, and 3.5, and look quickly at section 3.6. Let's do exercises 3.1 and 3.3 (which can be done on paper), and exercise 3.5 (it is sufficient to just map out what should be programmed) in section 3.8.
Summary: This week we spend most of our time discussing section 3.1 of LM19. This discussion involved defining the finite element basis functions, which are simply linear functions for first-order functions (P1), and quadratics for 2nd order (P2). These functions are defined only locally for each element (each element has 2 or more nodes associated with it). We how we can take integrals of combinations of these P1 or P2 functions to obtain "local" versions of the A matrix and b vector. These local versions can then be assembled into the "global" A matrix and b vector. We looked at implementation (section 3.2), and extensions into 2 and 3 dimensions (3.6)
Homework: For next time, we will continue with LM19, reading sections 4.1 (approximating differential equations), 5.1 (computing with finite elements) and 6.1 (time-dependence). We will not assign an exercise this week - our goal is to understand the reading.
Summary: We fininshed our reading of LM19, and discussed sections 4.1, 5.1, and 6.1, extending our analysis of basis functions to the solutions of differential equations. We first discussed how we can use basis functions generally to approximate solutions to differential equations (sections 4.1), including integration by parts to convert the differential equation into a "weak form" that allows approximation of the solution using basis functions are not second order differentiable. Next, we applied this understanding to the finite element basis function (section 5.1) and discussed how to construct a matrix system for P1 elements. Finally, we discussed how to set up finite differences in time to solve time-dependent problems with finite element solutions in space (section 6.1). This involves constructing both the mass matrix M and the stiffness matrix K. Both of these usually only need to be constructed once because they are dependent on the finite element mesh and physical properties of the material (e.g., viscosity or a diffusivity constant), which do not change with time for most problems.
Homework: We will meet next week to discuss the projects, so it would be helpful to get started on them this week so we can answer questions. We also discussed how it would be useful to do an exercise about constructing the mass and stiffness matrices. Here is an exercise:
Consider the linear system described in equation (6.16), which solves for time-dependent diffusion. Let's assume a diffusion coefficient of 1 mm^2/s and a domain width of 5 mm. First, break the domain into 5 equal-sized elements (with nodes at 0.0, 1.0, 2.0, 3.0, 4.0, and 5.0) and construct M and K. Next, let's assume that we need higher resolution on the left side, so we shift the nodes slightly (now nodes are at 0.0, 0.5, 1.2, 2.2, 3.5, and 5.0). Construct the M and K matrices for this non-uniform mesh.
Project writeups are due by 13.00 on Friday, May 28.
Introduction to Finite Differences