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.

1-D Linear Convection

Welcome to Space and Time — Introduction of Finite-difference solutions of PDEs, the second module of "Practical Numerical Methods with Python!". In the first module, we looked into numerical integration methods for the solution of ordinary differential equations, using the phugoid model of glider flight as a motivation. In this module, we will study the numerical solution of partial differential equations (PDEs), where the unknown is a multi-variate function, using the finite-difference method. The coming IPython notebooks will cover the linear and non-linear convection equations, the diffusion equation, and Burgers' equation. We hope you will enjoy them!

In this IPython notebook, we introduce the 1-D Linear Convection equation, which is the most basic model that can be used to learn something about hyperbolic PDEs. It is surprising that this little equation can teach us so much! Here it is:

\begin{equation}\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0\end{equation}

With given initial conditions (understood as a wave), the equation represents the propagation of that initial wave with speed $c$, without change of shape. Let the initial condition be $u(x,0)=u_0(x)$. Then the exact solution of the equation is $u(x,t)=u_0(x-ct)$.

We discretize this equation in both space and time, using the Forward Difference scheme for the time derivative and the Backward Difference scheme for the space derivative. Consider discretizing the spatial coordinate $x$ into points that we index from $i=0$ to $N$, and stepping in discrete time intervals of size $\Delta t$.

From the definition of a derivative (and simply removing the limit), we know that:

\begin{equation}\frac{\partial u}{\partial x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x}\end{equation}

Our discrete equation, then, is:

\begin{equation}\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n - u_{i-1}^n}{\Delta x} = 0, \end{equation}

where $n$ and $n+1$ are two consecutive steps in time, while $i-1$ and $i$ are two neighboring points of the discretized $x$ coordinate. If there are given initial conditions, then the only unknown in this discretization is $u_i^{n+1}$. We can solve for our unknown to get an equation that allows us to advance in time, as follows:

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

Now let's try implementing this in Python.

We'll start by importing a few libraries to help us out.

Now let's define a few variables; we want to define an evenly spaced grid of points within a spatial domain that is 2 units of length wide, i.e., $x_i\in(0,2)$. We'll define a variable nx, which will be the number of grid points we want and dx will be the distance between any pair of adjacent grid points.

We also need to set up our initial conditions. The initial velocity $u_0$ is given as $u = 2$ in the interval $0.5 \leq x \leq 1$ and $u = 1$ everywhere else in $(0,2)$ (i.e., a hat function).

Here, we use the function ones() defining a numpy array which is nx elements long with every value equal to 1.

Now let's take a look at those initial conditions using a Matplotlib plot. We've imported the matplotlib plotting library as plt and the plotting function is called plot, so we'll call plt.plot. To learn about the myriad possibilities of Matplotlib, explore the Gallery of example plots.

Here, we use the syntax for a simple 2D plot: plot(x,y), where the x values are evenly distributed grid points:

Why doesn't the hat function have perfectly straight sides? Think for a bit.

Now it's time to implement the discretization of the convection equation using a finite-difference scheme.

For every element of our array u, we need to perform the operation $u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$

We'll store the result in a new (temporary) array un, which will be the solution $u$ for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

We first initialize our placeholder array un to hold the values we calculate for the $n+1$ timestep, using once again the NumPy function ones().

Then, we may think we have two iterative operations: one in space and one in time (we'll learn differently later), so we'll start by nesting one loop inside the other. Note the use of the nifty range() function. When we write: for i in range(1,nx) we will iterate through the u array, but we'll be skipping the first element (the zero-th element). Why?

Note—We will learn later that the code as written above is quite inefficient, and there are better ways to write this, Python-style. But let's carry on.

Now let's try plotting our u array after advancing in time.

OK! So our hat function has definitely moved to the right, but it's no longer a hat. What's going on?

The solution differs from the expected hat because the discretized equation is an approximation of the continuous differential equation that we want to solve. This dissipative effect reduces the gradients on the solution, and we call it numerical diffusion or artificial viscosity.

Any numerical scheme is subject to numerical diffusion, but we can reduce its effect by using smaller mesh spacings. Try is out!