Simulating shockwaves and partial differential equations

Preamble: Run the cells below to import the necessary Python packages

Open this notebook in Google Colab:

ODEs vs PDEs¶

While there are many strategies available for integrating an ODE, typically when we first encounter a new ODE to integrate, we can use built-in integration tools. The only exception is in the case of special types of differential equations (delay, integro-differential, etc.) PDEs, however, usually are not so straightforward. Each new PDE we encounter will have its own phenomenology, strategies, and limitations, and so we often have to analyze the structure of a PDE before picking a numerical strategy for solving it.

However, common PDE (Navier Stokes, Laplace, etc.) often have significant research behind them, and there may be external libraries that work in certain regimes, or which can help with meshing, etc. But even in very restricted regimes, the design of these tools is expensive and narrow. Using finite element tools like Ansys, COMSOL, etc. is a whole specialization in itself.

However, when figuring out how to tackle a PDE, it can be helpful to identify its broad properties. A particularly important designation is whether the PDE is elliptical, hyperbolic, and parabolic.

where $u' = \frac{\partial u}{\partial x}$ and $\dot{u} = \frac{\partial u}{\partial t}$ in the table

We can see that the designation refers to the relationship between the order of time derivatives and space derivatives in the PDE. In general, in systems in which they are both the same order (such as first order in Burgers’ equation, or second order in the wave equation), the PDE will often exhibit waves and traveling solutions. In contrast when the derivatives are mixed, we tend to get spreading or coarsening solutions that distort as they change scales, as occurs in the diffusion equation.

In order to demonstrate the connections between simulating ODEs and simulating ODES, we will first implement ODE solver, similar to the ones used in previous lectures. We will structure this code with a base class BaseFixedStepIntegrator that will be inherited by specific solvers. This will allow us to easily switch between different solvers in the future.

Just to test that everything is working as intended, we can try integrating a simple system, a pendulum.

Determining boundary conditions¶

The boundary conditions are the conditions on the solution at of the PDE at edges of the domain, which includes both space and time for PDEs. In this sense, boundary conditions generalize the initial conditions we need to integrate ODEs. The boundary conditions are usually specified as functions of the integration variable(s) (like concentration or probability) as a function of the domain variables (e.g. starting or ending time, or position along an edge of the domain).

In a boundary value problem, we specify the boundary conditions as $u(\mathbf{x}, t) = u_0(\mathbf{x})$ for all $\mathbf{x} \in \partial \Omega$, where $\Omega$ denotes our solution domain (square, disk, etc). These are known as Dirichlet boundary conditions, where we specify the value of $u$ at the boundary. We can also have Neumann boundary conditions, where we instead specify the value of the derivative of $u$ at the boundary. For example, we could specify $\nabla u(\mathbf{x}, t) \cdot \mathbf{\hat{t}} = f_0(\mathbf{x})$ for all $\mathbf{x} \in \partial \Omega$, where $\mathbf{\hat{t}}$ is a unit vector normal to the boundary. This type of boundary condition specifies a fixed flux of the scalar field through the boundary. We can also have mixed or Robin boundary conditions, where we specify a combination of Dirichlet and Neumann boundary conditions. $u(\mathbf{x}, t) + \nabla u(\mathbf{x}, t) \cdot \mathbf{\hat{t}} = g_0(\mathbf{x})$ for all $\mathbf{x} \in \partial \Omega$.

Burgers’ equation¶

Burgers’ equation is nonlinear partial differential equation describing the time evolution of the speed $u(x, t)$ of a fluid in one dimension. This is a toy model of how steep waves form in shallow water in the ocean.

We will solve Burgers’ equation on a 1D domain of length $L$, so $x \in [0, L]$. Burger’s equation is a one-dimensional partial differential equation. The integration variable is $u(x,t)$, and the domain variables are $x, t$. Our boundary conditions and initial conditions will thus consist of the initial value of the field $u(x, 0)$ and the values of the velocity or acceleration at the edges of the domain $x$. Because this equation has first order derivatives in both space and time, we consider it a hyperbolic equation. The solution of Burgers’ equation tend to evolve towards shocks, which are sharp discontinuities in the velocity.

Semi-discretization: The Method of lines¶

In a semi-discrete method, we explicitly discretize the spatial domain using finite differences, spectral projection, or something else. We then solve the resulting set of coupled ODEs in continuous time using a built-in ODE solver. This solver might also discretize time (implicitly or explicitly); however, this is handled independently of how we choose to discretize space. Mathematically, we can think of discretizing space as splitting a continuous PDE into a large set of coupled ODEs. The initial field value becomes the initial conditions, and the boundary conditions become a condition on how the ODEs are coupled together. The number of coupled ODEs depends on how finely we discretize space.

We will solve semi-discretized PDEs is using a built-in ODE solver. We therefore separate meshing (dealing with the domain) from timestepping (evolving forward in time)

Finite difference operators¶

The finite difference operators on particular choice of spatial discretization scheme. Finite difference operators essentially split space into small pieces with size $\Delta x$. They then approximate spatial derivatives in the PDE on these discrete block.

The central first-order finite difference operators in 1D has the form

The central second-order finite difference operators in 1D has the form

We can see that, as the mesh becomes finer $\Delta x \to 0$, the quality of the approximation increases.

Semi-discretizing Burgers’ equation¶

Using the finite-different approximations to approximate the spatial derivatives, Burgers’ equation becomes

We have a hyperparameter to set, which is the number $N$ of spatial points that we use to discretize the domain. $\Delta x \equiv L / N$. After performing this discretization, the field variable becomes $\mathbf{u}(t) \in \mathbb{R}^N$. We can therefore write Burgers’ equation as a multivariate ODE, taking advantage of the fact that $u(x + \Delta x, t)$ corresponds to the next index into the vector $\mathbf{u}(t)$

The index $i \in 0, 1, ..., N$ tracks the different mesh points on the spatial site.

Boundary conditions¶

What happens when we get to the last value of the array? We can’t look up $u_{i+1}(t)$ if $i = L$. We now need to define, physically, what happens at the edges of our domain.

Periodic boundary conditions are simple to implement, but are not usually well-motived physically. The involve setting the condition that $u_{L+1}(t) = u_{0}(t)$, and that $u_{-1}(t) = u_{L}(t)$

Dirichlet boundary conditions describe many physical systems, but may be less stable numerically. The involve setting a fixed value, like $u_{L+1}(t) = 0$, and $u_{-1}(t) = 0$

Neumann boundary conditions are also very physical, and here they lead to mirror boundary conditions, where the walls of the domain reflect oncoming shocks. Since these involve conditions on the first derivative of the field with respect to space, we have to set these in terms of finite differences: $u_{L+1}(t) - u_{L-1}(t) = 0$, which is equivalent to the reflection symmetry $u_{L+1}(t) = u_{L-1}(t)$. The same goes for the other side of the domain, $u_{-1}(t) = u_{1}(t)$

We can now try numerically solving this system with the method of lines, in order to explore the behavior of Burgers’ equation.

Explore the behavior of this PDE¶

1. Try playing with the parameters $\nu$ and $\Delta x$ to see how the solution evolves. What causes the shocks to live longer?

2. Try changing the initial condition to a Gaussian. What happens?

3. Try changing the boundary conditions to reflective (Neumann) or Dirichlet. What happens?

4. What happens as we change the timestep $\Delta t$ or the number of mesh points $N_{mesh}$

The CFL condition for hyperbolic PDEs¶

The Courant–Friedrichs–Lewy (CFL) condition is a condition on mesh size $\Delta x$ and the time step $\Delta t$ that must be satisfied for the numerical solution of a hyperbolic PDE to be stable. Suppose that we have the wave equation in one dimension:

where $c$ determines the speed of travelling waves. A full discretization of this equation in space and time would be

The CFL condition determines the maximum value of $\Delta t$ that can be used for a given $\Delta x$ and $c$

where the dimensionless Courant number $C_{max} \sim 1$ is determined by our iteration scheme. The exact value of $C_{max}$ can be interpreted as the number of mesh cells that a particle of $\rho$ can traverse per integration time step (for an Euler time-stepping scheme and fixed mesh size, this is exact). Realistically, we want to stay a little bit below the CFL limit, because we don’t want to be right on the edge of instability. A typical “safety factor” would be $C = 0.7$.

Variants of the CFL condition also hold for other classes of PDEs. For example, the Neumann condition is a meshing condition for parabolic (diffusive) PDE: $D (\Delta t) / (\Delta x)^2 \sim 1$.

Applying the CFL condition to Burgers’ equation¶

While we typically use the CFL condition here for hyperbolic PDE (wave equations), in practice it’s a good rule-of-thumb for PDEs that admit travelling solutions, like solitons or Fisher waves of advance. In general, we can often find $c$ analytically for a PDE by defining the envelope $u(x)$ of a travelling solution wave, and then plugging the travelling solution $u(x - ct)$ into the PDE and solving for $c$. For Burgers’ equation, we can pick a well-defined speed $c$ associated with the shock. For our parameter choices above, $c \approx 3$.

It looks like our fronts travel about 1 spatial unit per time unit, and so we estimate that the relevant speed scale $c = 1$

If we set $N_{mesh} = 300$, then $\Delta x = \pi / 150$ and the CFL condition tells us that $\Delta t < \Delta x / c$, or that $\Delta t < \pi / 150 \approx 0.02$

The sparsity pattern of a PDE¶

When we discretize a PDE into a set of $N$ ODE, these ODE have coupling that depends on the geometry of the spatial domain and the discretization scheme. We can think of this relationship as defining a graph relating the various ODE that we are solving. For our finite-differences discretization of the first derivative, the adjacency matrix relating each of the N lattice sites to each other site looks like this:

The boundary conditions manifest as a condition on the first and last rows and columns of this graph. For example, the pattern above assumes that the field $u(x,t)$ is clamped at zero on the boundaries. In the original PDE, this is the condition $u(0, t) = u(L, t) = 0$, while in the spatially-discretized PDE, this corresponds to assuming that any index that reaches past the edges of the discretized domain will find a zero, $u_{-1} = u_{N+1} = 0$. However, if we wanted to capture periodic boundary conditions, then we would need to modify this matrix,

We can come up with similar procedures for reflecting or mixed boundary conditions. Since discretization removes explicit spatial variables from the problem, boundary conditions always appear as connectivity conditions on the resulting graph. Notice that our graph is highly sparse, with $\sim \mathcal{O}(N)$ nonzero elements. This sparsity can be exploited to avoid having to perform the $\sim \mathcal{O}(N^2)$ operations that are typically required when working with $N \times N$ matrices.

We can also write the second-order finite difference operator in a similar manner to the first-order finite difference operator: $$

$$

Notice that the second-order operator has elements along its diagonal. This is a consequence of the center value $u_i(t)$ appearing in the finite difference approximation of teh second derivative.

Burgers’ equation in the spectral basis¶

Recall that our central idea is to convert PDEs into a series of coupled ODES, and then solve this separately. Discretizing space is one way to do this. Each of the resulting ODEs describes a small, contiguous region of the domain.

However, we are free to pick other methods of partitioning the domain. An common alternative is to discretized our field by projecting it onto a series of basis functions $\phi_k(x)$ defined over the entire domain. If we think of the full spatial solution of a PDE as a superposition of many basis functions, the corresponding coupled ODEs describe the amplitude of each basis function over time.

A common example in classical physics is travelling waves. A sinusoidal travelling wave can be writted as a sum of two standing waves, whose relative amplitudes are each changing over time.

For Burgers’ equation, our original PDE was

Let’s define a set of Fourier basis functions,

where the angular wavenumber $k$ parametrizes the basis set. If our domain is $[0, L]$, then the values of $k$ become a countable set of frequencies that satisfy the periodic boundary conditions, $k = \frac{2 \pi n}{L}$. The Fourier basis functions are orthonormal, so we can write the solution as a linear combination of the basis functions

where $c_k(t)$ is the amplitude of the $k$ th Fourier mode at time $t$.

Let’s see what the profile of the Fourier basis functions, $c_k(t)$, looks like at various times. Notice that when we take the Fourier transform in Python, we truncate and take the first part of the array, since the second half is just the negative frequencies. Recall that these exist due to the phase of the complex exponential in the full definition of the Fourier transform. Since we only care about the magnitude of modes, we can just take the first half of the array returned by the Fourier transform.

We can now substituting the Fourier series expression for $u(x, t)$ into the Burgers’ equation PDE. We recall our rules for the Fourier transform of the derivatives of a function, which allow us to write the derivatives of $u(x,t)$ in terms of the amplitudes $c_k$ and wavenumbers $k$ of the Fourier modes

Substituting these into the PDE, we get

Notice that the nonlinear term is a bit more involved, since it involves a product of two sums. If we we working with a linear PDE, we wouldn’t have this complication. For example, if diffusion dominates the dynamics (or $\nu \gg 1$), then we would be able to split this equation term-by-term into a set of $N$ coupled ODEs in continuous time

For linear equations like the heat equation, the solution method is extremely simple: convert to the spectral domain, determine the time evolution of the coefficients $c_k(t)$, and then convert back to the physical domain. The modes of the heat equation are decoupled in the spectral domain, so we can solve the ODEs for each mode independently and very efficiently.

Nonlinear terms require pseudospectral methods¶

For nonlinear systems like Burgers’ equation or reaction-diffusion equations, it’s often not possible to convert the entire PDE into a concise spectral form. Instead, we need to “roundtrip” the solution between the physical and spectral domains at every evaluation. This is a bit more complicated, but it’s still possible to solve the problem efficiently. Methods that involve a roundtrip between the physical and spectral domains are called pseudo-spectral methods.

Spectral methods as preconditioning¶

Recall that the convergence time of linear dynamical systems was related to the condition number of the matrix. That’s because high-condition number matrices are closer to having degenerate column spaces, and thus act like lower-dimensional dynamical systems under repeated iteration. As a result, iterative algorithms like the power method are affected by a matrix’s condition number.

We can think of spectral methods as a way of preconditioning the system, by transforming the system into a basis where the condition number is lower. This is why spectral methods are so effective for linear partial differential equations. For nonlinear PDEs, the condition number of the system is not the only factor that determines the convergence time, but spectral transforms can sometimes allow us to solve the problem with larger timesteps.

What are the appropriate basis functions?¶

Usually, we pick basis functions that form an orthogonal set, and which are consistent with our boundary conditions. This is because we want to be able to represent the solution as a linear combination of the basis functions, and overlapping basis functions leads to degeneracy. Additionally, orthogonal basis functions usually lead to sparse coupling among the ODEs describing our discretized PDE.

So far, we’ve been using the Fourier basis functions, which describe systems with periodic boundary conditions. However, we can also use other basis functions, which describe systems with different boundary conditions.

Chebyshev basis functions¶

For Dirichlet boundary conditions (where the field value is clamped at a constant value at the boundary), we can use the Chebyshev basis functions, which are a set of functions that vanish or approach constant values at the boundary. The Chebyshev basis functions are a set of orthogonal functions that are defined on the interval $[-1, 1]$. They are defined recursively as

Legendre basis functions¶

Another option for Dirichlet boundary conditions, are the Legendre basis functions, which are a set of functions that vanish or approach constant values at the boundary. The Legendre basis functions are a set of orthogonal functions that are defined on the interval $[-1, 1]$. They are defined recursively as

Some specificities of Burgers’ equation¶

In Burgers’ equation, the shock waves are discontinuities in the derivative of the velocity, but not in the density. As with other PDEs with predominant advection, the shock waves propagate ballistically, i.e. at a constant speed of sound in the fluid. As a result, when the diffusion term is weak ($\nu \to 0$) certain initial conditions $u_0(x)$ in the shape of a shock wave will get carried with the flow, and will not change shape,

Burgers’ equation is a unique nonlinear PDE because it can be transformed into a linear PDE by introducing a new set of coordinates. The Cole-Hopf transformation consists of substituting the variable $u(x)$ with $\phi(x)$ using the following defintion

Inserting this into the viscous Burgers equation gives

After some algebraic manipulation, we get

This is a linear diffusion equation, which can easily be solved in one dimension. We can solve this simplified heat equation, and then cast back into the original variables by performing the inverse transformation

Unlike our chaotic systems, the complexity of the dynamics of the Burgers equation is not intrinsic, but rather due to our choice of dynamical variables. This is a preview of a common problem we will encounter in machine learning: is a problem intrinsically hard, or do we just need to pick the right coordinates?

Appendix¶

• This is William’s leftover code, which might get used in a future version of the tutorial