Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c) L.A. Barba, G.F. Forsyth, C. Cooper, 2014. Based on CFDPython, (c) L. A. Barba, also under CC-BY license.

Burgers' Equation

Hi there! We have reached the fifth and final IPython notebook of the series Space and Time — Introduction of Finite-difference solutions of PDEs, the second module of "Practical Numerical Methods with Python".

Throughout this module, we learned about the finite-difference solution for the linear and non-linear convection equations and the diffusion equation. It's time to combine all these into one: Burgers' equation. Before you continue, make sure you have completed the previous notebooks of this series, it will make your life easier. The wonders of code reuse!

You can read about Burgers' Equation on its wikipedia page.

Burgers' equation in one spatial dimension looks like this:

\begin{equation}\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2}\end{equation}

As you can see, it is a combination of non-linear convection and diffusion. It is surprising how much you learn from this neat little equation!

We can discretize it using the methods we've already detailed in the previous notebooks of this module. Using forward difference for time, backward difference for space and our 2nd-order method for the second derivatives yields:

\begin{equation}\frac{u_i^{n+1}-u_i^n}{\Delta t} + u_i^n \frac{u_i^n - u_{i-1}^n}{\Delta x} = \nu \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\end{equation}

As before, once we have an initial condition, the only unknown is $u_i^{n+1}$. We will step in time as follows:

\begin{equation}u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)\end{equation}

Initial and Boundary Conditions

To examine some interesting properties of Burgers' equation, it is helpful to use different initial and boundary conditions than we've been using for previous steps.

Our initial condition for this problem is going to be:

\begin{eqnarray} u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\ \phi &=& \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu} \bigg) \end{eqnarray}

This has an analytical solution, given by:

\begin{eqnarray} u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\ \phi &=& \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg) \end{eqnarray}

Our boundary condition will be:

\begin{equation}u(0) = u(2\pi)\end{equation}

This is called a periodic boundary condition. Pay attention! This will cause you a bit of headache if you don't tread carefully.

Saving Time with SymPy

The initial condition we're using for Burgers' Equation can be a bit of a pain to evaluate by hand. The derivative $\frac{\partial \phi}{\partial x}$ isn't too terribly difficult, but it would be easy to drop a sign or forget a factor of $x$ somewhere, so we're going to use SymPy to help us out.

SymPy is the symbolic math library for Python. It has a lot of the same symbolic math functionality as Mathematica with the added benefit that we can easily translate its results back into our Python calculations (it is also free and open source).

Start by loading the SymPy library, together with our favorite library, NumPy.

We're also going to tell SymPy that we want all of its output to be rendered using $\LaTeX$. This will make our Notebook beautiful!

Start by setting up symbolic variables for the three variables in our initial condition and then type out the full equation for $\phi$. We should get a nicely rendered version of our $\phi$ equation.

It's maybe a little small, but that looks right. Now to evaluate our partial derivative $\frac{\partial \phi}{\partial x}$ is a trivial task.

If you want to see the unrendered version, just use the Python print command.

Now what?

Now that we have the Pythonic version of our derivative, we can finish writing out the full initial condition equation and then translate it into a usable Python expression. For this, we'll use the lambdify function, which takes a SymPy symbolic equation and turns it into a callable function.

Lambdify

To lambdify this expression into a useable function, we tell lambdify which variables to request and the function we want to plug them in to.

Back to Burgers' Equation

Now that we have the initial conditions set up, we can proceed and finish setting up the problem. We can generate the plot of the initial condition using our lambdify-ed function.

This is definitely not the hat function we've been dealing with until now. We call it a "saw-tooth function". Let's proceed forward and see what happens.

Periodic Boundary Conditions

We will implement Burgers' equation with periodic boundary conditions. If you experiment with the linear and non-linear convection notebooks and make the simulation run longer (by increasing nt) you will notice that the wave will keep moving to the right until it no longer even shows up in the plot.

With periodic boundary conditions, when a point gets to the right-hand side of the frame, it wraps around back to the front of the frame.

Recall the discretization that we worked out at the beginning of this notebook:

\begin{equation}u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)\end{equation}

What does $u_{i+1}^n$ mean when $i$ is already at the end of the frame?

Think about this for a minute before proceeding.