FTCS scheme

The FTCS scheme is possibly the easiest scheme to understand, hence we discuss it first. We will discuss this scheme twice in basically the same way, firstly in a very explicit form, where it is obvious in the computational loop how the differential equation is being processed, and secondly in a matrix form which we will find is more powerful and enables us to easily understand the changes and extensions to the FTCS scheme introduced later.

FTCS stands for Forward Time Centred Space. What does this mean? Well, it means that we calculate the time derivative in a ``forward'' manner (using effectively Equation (3.2)) and calculate the spatial derivative in a ``centred'' manner (using effectively Equation (3.4) and/or Equation (3.5)). It is best to understand this scheme in the context of an example. Let us take for this example, the diffusion equation:

$$\frac{\partial }{\partial t} A(x,t) = \kappa \frac{\partial ^2}{\partial x^2} A(x,t).$$ (3.6)

In this equation, the scalar field^3.4 $A(x,t)$ is a function of both space and time. The variable $\kappa$ is just a constant describing the rate of diffusion of the field $A$. Note that we have a first derivative in time and a second derivative in space, one may think that it would be simpler to discuss as a first example an equation with just a first derivative in both time and space, but it turns out that such an equation (the form of which is known as the advection equation) is actually harder to calculate numerically than the diffusion equation. By the end of these notes you may either know the diffusion equation really well, or be completely sick of it, as we will use it in most of the sections and examples that follow.

Okay, to use the FTCS scheme on the diffusion equation, we need to calculate the first derivative in time using the forward first derivative equation (Equation (3.2)) and the second derivative in space using the centred second derivative equation (Equation (3.5)). Substituting these into Equation (3.6) we get:

$$\frac{A(x,t+\Delta t) - A(x,t)}{\Delta t} = \kappa \frac{A(x+h,t) + A(x-h,t) - 2 A(x,t)}{h^2}$$ (3.7)

which is the discretised form of the diffusion equation.

Let's simplify things slightly by introducing the notation that spatial indices are represented by a subscript, and time indices are represented by a superscript, such that the value of the function $A$ at the point $x_i$ and at time $t_n$ is represented as $A_i^n$, or in another form

$$A_i^n = A(x_i, t_n).$$ (3.8)

The discretised first derivative in time in this notation becomes

$$\frac{\partial }{\partial t} A(x,t) = \frac{A(x_i, t_{n+1}) - A(x_i, t_n)}{\Delta t} = \frac{A_i^{n+1} - A_i^n}{\Delta t},$$ (3.9)

the discretised second derivative in space becomes

$$\frac{\partial ^2}{\partial x^2} A(x,t) = \frac{A(x_{i+1},t_n) + A(x_{i-1},t_n) - 2 A(x_i,t_n)}{h^2} = \frac{A_{i+1}^n + A_{i-1}^n - 2 A_i^n}{h^2},$$ (3.10)

and the diffusion equation is now rewritten as:

$$\frac{A_i^{n+1} - A_i^n}{\Delta t} = \kappa \frac{A_{i+1}^n + A_{i-1}^n - 2 A_i^n}{h^2}.$$ (3.11)

Since we know the value of the function at the current point in time ($t_n$ which is $A(x_i,t_n) \forall i$), what we are really interested in is the value of the function at the next point in time ($t_{n+1}$, which is $A(x_i,t_{n+1}) \forall i$). We therefore solve Equation (3.11) for $A_i^{n+1}$ to obtain

$$A_i^{n+1} = A_i^n + \frac{\kappa \Delta t}{h^2} (A_{i+1}^n + A_{i-1}^n - 2 A_i^n),$$ (3.12)

which is the equation that we are now going to convert into Matlab code.

I have already written Matlab code to solve this differential equation (it is adapted from Professor Drummond's code, which was originally adapted from the code in Garcia), and I now give the code listing and will give a discussion of the various parts of the program line-by-line. The code listing for diffuse1.m is^3.5:

   1 % diffuse1.m - Program to solve the diffusion equation
   2 %              using the FTCS method.
   3 
   4 clear;             % clear memory
   5 help diffuse1.m    % print help information
   6 
   7 %%%%%%%% INITIALISATION %%%%%%%%%%
   8 
   9 % set up input variables
  10 tmax = 1;                          % window size in t
  11 nt = 100;                          % number of grid points in t
  12 delta_t = tmax/nt;                 % step-size in t
  13 nplots = 21;                       % number of plot points in t
  14 nplotstep = nt/nplots;             % plot step-size in t
  15 
  16 xmax = 2;                          % window size in x
  17 nx = 11;                           % number of grid points in x
  18 delta_x = xmax/(nx + 1);           % step-size in x
  19 
  20 kappa = 1;                         % coupling constant
  21 c = kappa*delta_t/delta_x^2;       % constant of propagation
  22 
  23 % initial conditions and boundary conditions
  24 A = zeros(nx,1);                   % allocates and zeros a column vector
  25 A(floor((nx+1)/2)) = 1/delta_x;    % sets up an initial delta 'spike' in A
  26 iplot = 1;                         % initialise the plotting index
  27 
  28 % set up plotting matrices
  29 xplot = (1:nx)*delta_x - xmax/2;   % row vector of x-coordinates
  30 tplot(iplot) = 0;                  % assign first time-coordinate
  31 Aplot(:,iplot) = A(:);             % assign first value of A
  32 
  33 %%%%%%%% MAIN LOOP %%%%%%%%%%
  34 
  35 % main loop
  36 for i = 1:nt
  37   
  38   % evolve our solution using the FTCS method
  39   A(2:(end-1)) = A(2:(end-1)) + c*( A(1:(end-2)) + A(3:end) - 2*A(2:(end-1)) );
  40   
  41   % save current data for plotting
  42   if (rem(i,nplotstep) < 1)
  43     
  44     iplot = iplot + 1;             % increment plotting index
  45     Aplot(:,iplot) = A(:);         % record current value of A for plotting
  46     tplot(iplot) = i*delta_t;      % record current time value
  47     
  48     % print some progress information
  49     fprintf('%g out of %g steps completed\n',i,nt);
  50     
  51   end
  52   
  53 end
  54 
  55 %%%%%%%% PLOTTING/OUTPUT %%%%%%%%%%
  56 
  57 % plot stuff
  58 figure(1)                          % opens a new figure window
  59 subplot(121)                       % divides window into a 1x2 matrix 
  60                                    % and assigns next plot to the 
  61                                    % left-hand frame (frame 1)
  62 mesh(tplot,xplot,Aplot)            % wire mesh plot of the solution
  63 view([-70 30])                     % change the viewing direction
  64 xlabel('t')                        % label the x-axis
  65 ylabel('x')                        % label the y-axis
  66 title('Diffusion of a delta spike')% give the thing a title
  67 
  68 subplot(122)                       % assigns next plot to right-hand 
  69                                    % frame (frame 2)
  70 
  71 % find contours of 3D plot and assign to the vector cs
  72 cs = contour(xplot,tplot,flipud(rot90(Aplot)),0:0.1:2.5);
  73 xlabel('x')                        % label the x-axis
  74 ylabel('t')                        % label the y-axis
  75 clabel(cs,0:5)                     % label the contours
  76 title('Contour of a delta spike')  % give the thing a title
  77 
  78 subplot(111)                       % return plotting mode to one 
  79                                    % figure per window

Notice that the code is separated into basically three parts: an initialisation section, a main computational loop, and a plotting section. This is a very simple structure that is quite useful in the programs you will be writing in this section of the course.

Going through the code line-by-line we have

lines 1-2:

give the name of the program and what it does.

line 4:

clear the memory (a good idea to do at the start of each program).

line 5:

prints help information. What this does is print any lines of code that appear at the top of the file in comments before any commands are entered. So what it does in this case is merely echo the top two lines of the file out to the terminal window.

lines 9-21:

intialise input variables. This code is hopefully fairly self explanatory. A couple of points to note though: we are not using a variable called $h$ for the step size in $x$, but $\Delta x$; and we have wrapped up all of the constant terms $\kappa \Delta t/\Delta x^2$ into a variable called c. This variable can give us an indication of whether the solution will be stable or not (more on this point later).

line 24:

this line merely allocates some memory for the vector A. We allocate its memory by calling the zeros function with the size of each dimension of the array as its arguments.

line 25:

we now set up an initial delta spike in A (basically a discrete version of a Dirac delta function). The delta spike is dependent upon the step size in $x$ only. We put the delta spike at the mid-point of the array by using the value offloor((nx+1)/2). What does the floor function do? Well, it looks at the number in its argument (in this case, the value of(nx+1)/2) and if the number isn't an integer, it rounds it down to the nearest integer. Why do we want to do this? Well, Matlab wants integer indices to its arrays, and it's possible that(nx+1)/2 could be a non-integer fraction. Just for your information, there is also a function called ceil that takes the ceiling of the number, i.e. it always rounds a number up to the nearest integer.

line 29:

this sets up the vector of $x$-coordinate positions. Note that we are using the Matlab colon operator (:) to generate the vector. The notation 1:nx means to create a vector from 1 to nx in integer steps. In the case shown here, we are setting up our $x$-coordinates to be from -xmax/2 toxmax/2.

line 31:

here we store the initial value of A into our plotting storage vector Aplot.

line 39:

now we are in the main loop. This is the line that does all of the hard work and implements the evolution of the vectorA according to the diffusion equation. Compare the code we are using here with Equation (3.12). You may have noticed that we are not assigning to the outside points of A. Why is this? We don't assign to these points to maintain the boundary conditions. In this example we are implicitly using the Dirichlet boundary conditions that the function value is zero at the boundary. To put this into a physical context, say we are investigating the heating of a metal bar. The diffusion equation describes the way in which the temperature changes with time and position along the bar. The ends of the bar we keep at zero temperature, hence the boundary conditions of our PDE. What is cool^3.6 about this line of code is that all of the points in space are being acted upon simultaneously, because we are using Matlab's array slicing notation (i.e. in the bit where we say A(2:end-1), we are slicing out only the values of A from the second element to the second to last element, and then using only these elements in the calculation). The first term in the equation on line 39 corresponds directly to the first term in Equation (3.12). As mentioned before, the variable c rolls all of the various constants together. The next term corresponds to the $A_{i+1}^n + A_{i-1}^n - 2 A_i^n$ part of Equation (3.12). So, the code A(1:(end-2)) means to say ``all of the left-hand $x$ points'' i.e. $A_{i-1}^n$. The code A(3:end) means ``all of the right-hand $x$ points'' i.e. $A_{i+1}^n$. And the code 2*A(2:end-1) means ``all of the centre $x$ points'' i.e. $A_i^n$. A pictorial explanation may be helpful at this point.

Figure: The use of element by element addition in calculation of the centred second derivative in space in the FTCS scheme. The vector $\mbf{A}$ is split into three parts: $\mbf{A}_\mrm{left}$ representing the left-hand elements of $\mbf{A}$; $\mbf{A}_\mrm{centre}$ representing the centre elements of $\mbf{A}$, each multiplied by 2; and $\mbf{A}_\mrm{right}$ representing the right-hand elements of $\mrm{A}$. Each element is added to its corresponding partner shown by the arrows in the figure, to then give the elements of the new value of $\mbf{A}$ which is denoted here $\mbf{A}_\mrm{new}$.

Figure 3.2 shows the process of element by element addition of parts the vector $\mbf{A}$ in order to calculate the centred second derivative in space in the FTCS scheme. The vector $\mbf{A}$ is split into three parts: $\mbf{A}_\mrm{left}$ representing the left-hand elements of $\mbf{A}$; $\mbf{A}_\mrm{centre}$ representing the centre elements of $\mbf{A}$, each multiplied by 2; and $\mbf{A}_\mrm{right}$ representing the right-hand elements of $\mrm{A}$. All elements of these vectors are added to corresponding elements in the other vectors as shown by the arrows in the figure. These element by element processes then produce the elements of $\mbf{A}_\mrm{new}$ which is the new value of the vector $\mbf{A}$ at the next time step. Notice that the end elements of $\mbf{A}$ and $\mbf{A}_\mrm{new}$ are set to zero. This is the implementation of the Dirichlet boundary conditions, and is why the intermediate vectors have a length with is two elements shorter than $\mbf{A}$.

line 42:

we now do a check to see if we should take a point for plotting^3.7. We defined the variablenplotstep in the initialisation section so that we can get 21 samples of the function evolution through time. The function we use here is the rem() function, and it calculates the remainder ofi/nplotstep. What we are doing here is using the if statement to check to see if this remainder is less than one, if so, to then take a plot point.

lines 44-46:

here we just take the value of A and assign it to part of the array Aplot and add to the time recording vector (tplot) with the current time.

line 49:

this just prints out the number of steps out of the total number of steps we have completed. It's sometimes helpful to have progress information like this to make sure your program is doing something when it is running.

line 58:

this opens a new figure window, and will call itFigure 1 (hence the ``1'').

line 59:

this splits the window into two parts, so that we can plot something first on the left half of the figure, and then on the right half. The notation defines a matrix, namely: subplot(121) means make a subplot of the current figure window in a $1 \times 2$ matrix (the first 1 and 2), and put the next plot into the first element of that matrix (hence the last 1, and the first element of a $1 \times 2$ matrix is the left-hand one, and so the plot goes into the left hand half of the figure window).

line 62:

make a mesh plot of the array Aplot using the vectors tplot and xplot along the $x$ and $y$ axes respectively.

line 63:

change the angle that we look at the graph from.

lines 64-66:

put $x$ and $y$ labels on the graph, and give it a title.

line 68:

same idea as line 59, just this time we are plotting into the second element (hence the 2 at the end of the subplot command).

line 72:

generate a variable used for plotting contours of the array Aplot. Don't worry about this too much, if you're interested type help contour for Matlab's information on thecontour function.

line 78:

subplot(111) resets any subsequent plotting output to take up the whole window: a $1 \times 1$ matrix, first element.

The stability of this scheme for the diffusion equation can be shown to be dependent upon the value of the variable c. It is possible to show (and for the analysis on this point see Garcia) that for $\mathtt{}{c} \ge 0.5$ the solution will be unstable. In fact, it is possible to show for the advection equation that the FTCS scheme is unstable for all values of the time step!^3.8 One can think of the evolution of the function according to the discretised diffusion equation as being due to perturbations on the currently known step in time. If you look at Equation (3.12), you can see that the new point in time is basically the old point in time plus some perturbation. It isn't desirable in perturbation theory to have large perturbations, so if c is too big, then the perturbation is large and the solution is likely to diverge as we iterate it further. So, when you are running the sample code and you find for certain parameters that the solution ``blows up'', check the value of c and see what it is; it is quite possible that it is greater than 0.5. To fix this one needs to either increase the $x$ step width, and/or decrease the time step or coupling constant. Why is this do you reckon? Have a think on it, and look at the diffusion equation again.

Hopefully now you have a reasonable grasp of the FTCS scheme. I'll now introduce a matrix notation that will help us understand subsequent schemes.

Footnotes

... field^3.4

In general, one would work with vector fields, such that $x$ would actually be a vector of various (say, spatial) dimensions, it is easier to learn this stuff if we are content just to consider scalar fields.

... is^3.5

It is also given as the sample code at the back, but here I want to discuss it line-by-line, hence I include it here with line numbers.

... cool^3.6

Believe it or not, pun not intended.

... plotting^3.7

Have a think about why we may not want to take all of the data points.

... step!^3.8

Hence, we use different techniques to treat the advection equation.