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:

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[ref]Mathworks and youtube both have excellent resources for a number of techniques and implementations respectively.[/ref].

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.

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.

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).

!========== OCP5930 Assignment 2_FTCS Source ===========!
! Author	Thomas Bryce Kelly
! Email		tbk14@my.fsu.edu
! See this and more code over at <http://misc.tkelly.org/>
! Checkout https://tkelly.org/ for other information
!
! Copyright	GNU GENERAL PUBLIC LICENSE V3.0
!			see <http://www.gnu.org/licenses/gpl.html> for details.
! ==================================================!

program main

real*8, dimension(25,2) :: v
real*8 c, dx, dt, pi, r		        !Declaration of variables
logical::exist				!Does the output file exist
integer length

!===================== Set Values=============================!
length = 25
pi = 3.1415926535897			!define
c = 0.1				!Give a value
dx = 1.0/25.0				!Let = 1/25 -- dimensionless
dt = dx/(20*c)				!Define so that CFL criteria satisfied -- dimensionless
!==============================================================!

r =c*dt/dx
write (*,*) "r = ", r, "dx = ", dx, "dt = ", dt

!~~~~~~~~~~~~~~~~~~~~~~~~File IO Setup~~~~~~~~~~~~~~~~~~~~~~~~~!
!Time to open the output file 
inquire(file="assignment_2_FTCS.csv", exist=exist)		!Test if file exists
!Dont want IO error...
if(exist) then
  ! Erase the old file before creating new one.
  open(12, file="assignment_2_FTCS.csv", status="old")
  close(12, status="delete")
endif
! Create new file 
open(12, file="assignment_2_FTCS.csv", status="new", action="write")
!Write the header
write (12,*) "X,Y,t"
!===========================Data Time===========================!
!Initialize array
do i=1,length
  v(i,1) = sin(2*pi*(i*dx))+1
  write(12,*) 2*i*dx-1,",", v(i,1),",", 0
end do

!Forward in time:
do k=1,10000	!Loop for each time: t=k*dt
  !new first point = old first point + c*dt*(far right - 2* old first point + old point two)/dx^2
  v(1,2) = v(1,1)+r*( -1*v(length,1)+v(2,1) )			!Left end
  write (12,*) 2*dx-1,",", v(1,2),",", k*dt
  
  !write (*,*) "The slope on the left end point is ", v(length,k)-2*v(1,k)+v(2,k)
  
  do i = 2,length-1	!Loop for each point
    ! [new point] = [oldpoint]+ -c*dt([old point - 1]+[oldpoint+1]-2[old point])/dx^2
    v(i, 2) = v(i,1)+r*(-1*v(i-1,1)+v(i+1,1))
    !Write the value
    !write (*,*) i*dx,",", v(i,2),",", k*dt
    write (12,*) 2*i*dx-1,",", v(i,2),",", k*dt
  end do

  !Right end point must be done separately (special case)
  v(length,2) = v(length,1)+r*( -1*v(length-1,1)+v(1,1) )			!Right 
  write (12,*) 2*length*dx-1, ",", v(length,2),",", k*dt
  
  do j=1,length
  	v(j,1)=v(j,2)
  end do
end do

close(12)	!Close the output file 
stop
end program

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 [ref]Check out my post on the value of R (here) and a site that helped me to make my R script (here).[/ref] , 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:

#<-<-<-<-<-<-<-<-<-<- OCP5930 Assignment 2 R Script <-<-<-<-<-<-<-<-<-<-<-!
# Author Thomas Bryce Kelly
# Email tbk14@my.fsu.edu
# See this and more code over at <https://github.com/tbrycekelly>
# Checkout https://tkelly.org/ for other information
#
# Copywrite GNU GENERAL PUBLIC LICENSE V3.0
# see <http://www.gnu.org/licenses/gpl.html> for details.
# <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-!

print("Loading FTCS Source.")
FTCS <- read.csv("~/Dropbox/FSU/OCP5930/Fortran/assignment_2_FTCS.csv", header<-TRUE,sep<-",")
print("Loading CTCS Source.")
CTCS <- read.csv("~/Dropbox/FSU/OCP5930/Fortran/assignment_2_CTCS.csv", header<-TRUE,sep<-",")
print("Loading analytical Source.")
ANA <- read.csv("~/Dropbox/FSU/OCP5930/Fortran/assignment_1.csv", header<-TRUE,sep<-",")

delta <- 0.02

frames <- as.numeric(commandArgs(TRUE)[1])
step <- as.numeric(commandArgs(TRUE)[2])

err_SINE <- 4.60979
err_STEP<-3.0

CTCS_drift <- vector()
FTCS_drift <- vector()
CTCS_mu <- vector()
FTCS_mu <- vector()
CTCS_max <- vector()
FTCS_max <- vector()

setwd('Some Directory')#Where I want the stuff

print("Analyzing Data")

for(i in 0:frames){
 if(i %% step <-<- 0){
 # creating a name for each plot file with leading zeros
 if (i < 10) {name <- paste('0000',i,'plot.png',sep<-'')}
 if (i < 100 && i ><- 10) {name <- paste('000',i,'plot.png', sep<-'')}
 if (i<1000 && i ><- 100) {name <- paste('00', i,'plot.png', sep<-'')}
 if (i<10000 && i ><- 1000) {name <- paste('0', i,'plot.png', sep<-'')}
 if (i<100000 && i ><- 10000) {name <- paste('', i,'plot.png', sep<-'')}
 
 #saves the plot as a .png file in the working directory
 png(name)
 
 l_FTCS <- subset(FTCS, FTCS$t < delta*(i+0.5) & FTCS$t > delta*(i-0.5))
 l_CTCS <- subset(CTCS, CTCS$t < delta*(i+0.5) & CTCS$t > delta*(i-0.5)) 
 l_ANA <- subset(ANA, ANA$t < delta*(i+0.5) & ANA$t > delta*(i-0.5)) 
 
 
 mu_CTCS <- 0.0
 mu_FTCS <- 0.0
 
 CTCS_MX <- -Inf
 CTCS_MY <- -Inf
 FTCS_MA <- -Inf
 FTCS_MY <- -Inf
 ANA_MX <- -Inf
 ANA_MY <- -Inf
 
 for(j in 2:25) {
 #Update Mu values
 if (j >1) {
 mu_CTCS <- mu_CTCS + sqrt((l_CTCS[j,1]-l_CTCS[j-1,1])**2 + (l_CTCS[j,2]-l_CTCS[j-1,2])**2)
 mu_FTCS <- mu_FTCS + sqrt((l_FTCS[j,1]-l_FTCS[j-1,1])**2 + (l_FTCS[j,2]-l_FTCS[j-1,2])**2)
 }
 
 #Check for new maximums
 if (l_FTCS[j,2] > FTCS_MY) {
 FTCS_MY <- l_FTCS[j,2]
 FTCS_MX <- j
 }
 if (l_CTCS[j,2] > CTCS_MY) {
 CTCS_MY <- l_CTCS[j,2]
 CTCS_MX <- j
 }
 if (l_ANA[j,2] > ANA_MY) {
 ANA_MY <- l_ANA[j,2]
 ANA_MX <- j
 }
 
 }
 #Update Global Data
 CTCS_drift <- c(CTCS_drift, (CTCS_MX-ANA_MX))
 FTCS_drift <- c(FTCS_drift, FTCS_MX-ANA_MX)
 CTCS_mu <- c(CTCS_mu, mu_CTCS)
 FTCS_mu <- c(FTCS_mu, mu_FTCS)
 CTCS_max <- c(CTCS_max, CTCS_MY)
 FTCS_max <- c(FTCS_max, FTCS_MY)
 
 plot(l_ANA$X, l_ANA$Y, type<-'l', xlim <- c(-1,1), ylim <- c(0,2), ylab <-'Y Value', xlab<-'X Value', main <- paste('Signal at timestep ', i), col <- 'black')
 lines(l_CTCS$X, l_CTCS$Y, type<-'l', col <- 'blue')
 lines(l_FTCS$X, l_FTCS$Y, type<-'l', col <- 'red')
 
 text(-0.75, 2, paste("CTCS mu: ",format(mu_CTCS/err_SINE, digits<-3)))
 text(-0.75, 1.8, paste("FTCS mu: ",format(mu_FTCS/err_SINE, digits<-3)))
 dev.off()
 }
}

png("Phase_Shift.png")
plot(FTCS_drift, main<-"Phase shift", xlab<-"Time-step", ylab<-"Shift (x)", col<-"red", type<-"l")
lines(CTCS_drift, col<-"blue")
dev.off

png("MU_Shift.png")
plot(FTCS_mu, main<-"Mu shift", xlab<-"Time-step", ylab<-"Shift (mu)", col<-"red", type<-"l")
lines(CTCS_mu, col<-"blue")
dev.off

png("Amplitude_Shift.png")
plot(FTCS_max, main<-"Amplitude shift", xlab<-"Time-step", ylab<-"Shift (x)", col<-"red", type<-"l")
lines(CTCS_max, col<-"blue")
dev.off

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

convert ./IMG/*plot.png ./out.mp4

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.

!========== OCP5930 Assignment 2_CTCS Source ===========!
! Author Thomas Bryce Kelly
! Email tbk14@my.fsu.edu
! See this and more code over at <http://misc.tkelly.org/>
! Checkout https://tkelly.org/ for other information
!
! Copywrite GNU GENERAL PUBLIC LICENSE V3.0
! see <http://www.gnu.org/licenses/gpl.html> for details.
! ==================================================!

program main

real*4, dimension(1:25,1:3) :: v
real*4 c, dx, dt, pi, r, sam !Declaration of variables
logical::exist !Does the output file exist
integer length, k

!===================== Set Values=============================!
length = 25 !Number of Points
pi = 3.1415926535897 !Define <pi>
c = 0.1 !Give <c> a value
dx = 1.0/25.0 !Let <dx> = 1/25 -- dimensionless
dt = dx/(20*c) !Define <dt> so that CFL criteria satisfied -- dimensionless
r =c*dt/dx !Helper Coefficient
!==============================================================!

!Print reference info:
Write (*,*) "Constants used: r = ", r," dx = ", dx, " dt = ", dt

!===========================File IO Setup======================!
!Time to open the output file <assignment_2_CTCS.csv>
inquire(file="assignment_2_CTCS.csv", exist=exist) !Test if file exists

!Dont want IO error...
if(exist) then
 !Erase the old file before creating new one.
 open(12, file="assignment_2_CTCS.csv", status="old")
 close(12, status="delete")
endif

! Create new file <assignment_2_CTCS.csv>
open(12, file="assignment_2_CTCS.csv", status="new", action="write")

!Write the header
write (12,*) "X,Y,t"

!===========================Data Time==========================!

!Initialize array

!This is Sine initialization values:
do i=1,length
 v(i,1) = sin(2*pi*(i*dx))+1
 write(12,*) 2*i*dx-1,",", v(i,1),",", 0
end do

!This is Square wave initialization values
!do i=1,length
! if (i > 0.25 * length) then
! v(i,1) = 0.0
! else
! v(i,1) = 1.0
! end if
! write(12,*) 2*i*dx-1,",", v(i,1),",", 0
!end do

!Compute the first time step (Forward in Time)
 !new first point = old first point + c*dt*(far right - 2* old first point + old point two)/dx^2
 v(1,2) = v(1,1)+r*( -1*v(length,1)+v(2,1) ) !Left end
 write (12,*) 2*dx-1,",", v(1,2),",", k*dt
 
 !Middle section calculations
 do i = 2,length-1 !Loop for each point
 ! [new point] = [old point] + -c*dt([old point - 1]+[old point + 1] - 2[old point])/dx^2
 v(i, 2) = v(i,1)+r*(-1*v(i-1,1)+v(i+1,1))
 !Write the value
 write (12,*) 2*i*dx-1,",", v(i,2),",", k*dt
 end do

 !Right end point calculations
 v(length,2) = v(length,1)+r*( -1*v(length-1,1)+v(1,1) )
 write (12,*) 2*length*dx-1, ",", v(length,2),",", k*dt


!Centered in time (main runtime):
do k=2,1000000 !Loop for each time: t=k*dt
 !new first point = old old first point + c*dt*(far right - 2* old first point + old point two)/dx^2
 v(1,3) = (0.02*v(1,2)+0.98*v(1,1))+r*( -1*v(length,2)+v(2,2) ) !Left end
 write (12,*) 2*dx-1,",", v(1,3),",", k*dt

 do i = 2,length-1 !Loop for each point
 ! [new point] = [old old point] + -c*dt([old point - 1]+[old point + 1] - 2[old point])/dx^2
 v(i, 3) = (0.02*v(i,2)+0.98*v(i,1))+r*(-1*v(i-1,2)+v(i+1,2))
 !Write the value
 write (12,*) 2*i*dx-1,",", v(i,3),",", k*dt
 end do

 !Right end point must be done seperately (special case)
 v(length,3) = (0.02*v(length,2)+0.98*v(length,1))+r*( -1*v(length-1,1)+v(1,1) ) !Right endss
 write (12,*) 2*length*dx-1, ",", v(length,3),",", k*dt
 do i=1,length
 v(i,1) = v(i,2)
 v(i,2) = v(i,3)
 end do
 
end do

close(12) !Close the output file <assignment_2_CTCS.csv>
stop
end program

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.

#/bin/bash

#========== OCP5930 Assignment 2 Bash Script ===========!
# Author Thomas Bryce Kelly
# Email tbk14@my.fsu.edu
# See this and more code over at <https://github.com/tbrycekelly>
# Checkout https://tkelly.org/ for other information
#
# Copywrite GNU GENERAL PUBLIC LICENSE V3.0
# see <http://www.gnu.org/licenses/gpl.html> for details.
# ==================================================!

COUNT=2000
STEP=2

echo "Compile CTCS"
gfortran ./assignment_2_CTCS.f90
./a.out

echo "Compiling FTCS"
gfortran ./assignment_2_FTCS.f90
./a.out

echo "Compiling Analytical"
gfortran ./Assignment_1_Source.f90
./a.out

echo "Clearing output directory"
rm ./IMG/*.png

echo "Making images..."

Rscript Assignment_2.R $COUNT $STEP

echo "Compiling movie"
convert ./IMG/*plot.png ./out.mp4

echo "\\\\\\\\\\\\\\\\\\\\"
echo "!------DONE!--------!"
echo "////////////////////"

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 [latex]\frac{\partial f}{\partial t} = \frac{\partial f}{\partial x}[/latex], 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 [latex]x\in(-1,1)[/latex] 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. [latex]x=n\cdot dx[/latex] where n is an integer and dx is the distance between points on our grid—think of graph paper).

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

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

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

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

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…

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

Substituting back:

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

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.