I want to share with you a step by step guide to how I structured, simulated, and compiled my results for comparing two basic finite difference approximation schema. The system which we will be modeling is given by this differential equation:

Screen Shot 2015-01-27 at 09.24.10

You are not alone if you are feeling uneasy upon looking at this equation; but don’t worry, we’ll get through it together. This equation is classified as a ‘hyperbolic, linear, first order differential equation’ which simply means (A) hyperbolic means that the solution is given by the initial values given to it—we’ll talk about it—; (B) linear means that there are no exponents; (C) and first order means the derivative is just the first derivative.

I promise we won’t dwell on the maths here, but we need a foundation of reference before we walk through it conceptually. The outline of this post will be (1)How to generate the data—FTCS, (2) How we can analyze that data,  (3) How to generate a better set of data—CTCS, and (3) How to make our lives easier. There’s a technical Addendum to follow.

Before we get too far, here are some of my results:

Generating the Data

To generate the data points that I will later use to form plots, animations, and analyses of the schema, I used Fortran. Fortran may seem like an odd choice of programing language, but it is still regularly used in the majority of fluid dynamic models to do the heavy numerical lifting. For this analysis, I’ll be looking at the output of two different Finite Difference Approximation (FDA) schema: Forward in Time, Centered in Space (FTCS) and Centered in Time, Centered in Space (CTCS). Don’t let the names and acronyms be worrisome, they are quite simple to use once you see an example.

While there is a tremendous amount of mathematical underpinnings for FDAs, I am neither willing nor qualified to discuss it here. Instead here are several references(1)Mathworks and youtube both have excellent resources for a number of techniques and implementations respectively..

FTCS Scheme

First we will look into the FTCS scheme since you can imagine this to be a basic version of the CTCS scheme that we’ll see later. As an Initial Value Problem (IVP), at each time step the updated solution is dictated by the previous solution. In other words, if we the first solution, say f(0), then we can find f(1), and then with f(1), we can get f(2), and so on for as long as we like. Since this is an approximation, we will also be interested in understanding how accurate our solution is as a function of time (i.e. Does the difference grow ad infinitum, does it stay fixed, or decrease?)

Now onto the algorithm. To keep from getting technical, I’ll simply show the formula and then tell you what it does.

Screen Shot 2015-01-27 at 08.49.55

The notation here is simply f(x=i, t=n+1), so we have that the point i in the next time step, n+1, is the sum of the old value of point i, f(i,n), and the scary looking bit. The right hand side is simply the change in f over the change in x times our time step. This is simply a result provided by the differential equation we are trying to model. Here is the same equation generalized.

Screen Shot 2015-01-27 at 08.55.24

If this isn’t making sense and you would like to know more about this, check out the addendum at the bottom.

Armed with the first algorithm, we can now take a look at the Fortran code that will initialize our first time step, f(0), and then progress forward in time one step at a time using this algorithm. (Feel free to skip past this block if you don’t care about the code).

This program will compute and generate an output csv file of all the points for as many time steps as you want. Here it is set for 10,000 steps, but I have run it for up to 10,000,000 without issue.

Data Analysis

Now that we have our raw data, we’ll need to analyze it and create graphical representations so that other people can understand our results. For this I am using an R script I wrote (2)Check out my post on the value of R (here) and a site that helped me to make my R script (here). , and what it does it read in the data and then for each time step it generates a plot and exports it. That way we can then run a quick compile call to stack all the plots on top of each other and make a animation. Without further ado, here is the R script code:

And with one more step, I can have my animation.

And here we go:

In this animation the black line is the analytical (exact) solution, the red line is the FTCS approximation that we have been working with, and the blue line is the CTCS approximation (see below). Perhaps you can spot some trends.

CTCS Scheme

The CTCS scheme can be seen as a simple extension of the FTCS scheme where you included information from the last two time steps instead of just the last one. This has the benefit of being stable (whereas the FTCS is unstable and eventually blows up). In the animation above, you can see that the CTCS signal propagates in a stable fashion whereas the FTCS signal starts expanding. Checkout the Fortran code for the CTCS scheme below.

 

Automating

One last thing we might want to consider if we’re planning on running this set of programs any number of times under different conditions or values would be to automate this workflow a bit. To that end, I wrote a simple and straightforward shell script.

All of the code here, and much more, can be found over in my misc directory (learn more about it here).


 

Addendum

Let’s start from the beginning. We are asked to model \frac{\partial f}{\partial t} = \frac{\partial f}{\partial x}, and in order to do so we need two pieces of information: what’s the domain we’re going to use and what are the initial conditions.

For the domain we will use the cyclic domain of x\in(-1,1) wherein the domain is rolled up onto itself so that the point to the left of -1 is actual next to +1. Sorry if I’m not explaining it very well, but essentially it is like playing some of those old school games where if you try to move off the screen you reappear on the opposite side. For the second piece of information, we will initialize with a simple sine wave; but just as with our choice of domain, there are interesting dynamics that occur with other decisions too.

Let’s walk through the formula’s we’ll use of the FTCS (Forward in Time, Centered in Space) scheme. Since this is a simulation we’ll be using discretized values (i.e. x=n\cdot dx where n is an integer and dx is the distance between points on our grid—think of graph paper).

Let’s say that f(n\Delta x,m\Delta t) \mapsto f(i, j) then in general we can write our finite difference equation as:

f(i,j+1)=f(i,j) + \Delta t (\frac{\partial f}{\partial t})

The only source of information available to us in regards to the value of \frac{\partial f}{\partial t} comes from our differential equation. Substituting in we then have:

f(i,j+1)=f(i,j) + \Delta t (\frac{\partial f}{\partial x})

This change may not look like anything at first, but take a look at how it can be rewritten in terms of f. Since the derivative of f with respect to x is simply the slope, let’s do just that: rise over run…

\frac{\partial u}{\partial x} \approx \frac{f(i+1,j)-f(i-1,j)}{2\cdot\Delta x}

Substituting back:

f(i,j+1)=f(i,j) + \Delta t (\frac{f(i+1,j)-f(i-1,j)}{2\cdot\Delta x})

And that’s about it. Looking at this equation we can talk about it conceptually as:

To find the value of point i in the next time step, let’s take the old value of i and add some quantity that will estimate it’s change in time. In other words, try to figure out where that value was going (increase/decrease). To do that we are using the relationship between its change in time and change in space and the two closest points to calculate it.

My narrative might not be perfect, but I hope you can see that we’re using information contained within three points at time-step k to calculate the value of a point in time-step k+1.

 

Notes   [ + ]

1. Mathworks and youtube both have excellent resources for a number of techniques and implementations respectively.
2. Check out my post on the value of R (here) and a site that helped me to make my R script (here).