Matrix form of FTCS scheme

To find the matrix form of the FTCS scheme we first recall the discretised form of the FTCS scheme:

$$\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.13)

If we stare at this equation for long enough we might be able to see exactly why we can turn this into a matrix form. Let's pay attention to just the right-hand side of Equation (3.13). We can see that we have three terms of the solution $A(x,t)$ at three points in space, but at the same point in time. This suggests that we could write $A(x,t)$ as a vector with each element representing the function value at each position in space and for the current point in time. Also, to ``pick out'' the relevant points from the vector $\mbf{A}$ we could multiply it by a square matrix^3.9 with relevant numbers down its diagonals. To make this more clear, let's consider the following matrix:

$$\begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix},$$ (3.14)

if we then multiply this matrix by the column vector:

$$\begin{pmatrix}v_1\\ v_2\\ v_3 \end{pmatrix}$$ (3.15)

we get out

$$\begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix} ... ...1\\ v_2\\ v_3 \end{pmatrix} = \begin{pmatrix}v_1\\ v_2\\ v_3 \end{pmatrix}.$$ (3.16)

At which point you go ``duh, Paul, that's because it was the identity''. Well, yes, it was. However, if we wanted to ``pick out'' the bottom two elements of the vector we would use the following matrix:

$$\begin{pmatrix}0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0 \end{pmatrix}$$ (3.17)

and we get

$$\begin{pmatrix}0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0 \end{pmatrix} ... ...v_1\\ v_2\\ v_3 \end{pmatrix} = \begin{pmatrix}v_2\\ v_3\\ 0 \end{pmatrix}.$$ (3.18)

So what? you may ask. Well, look at the vector output. It has what were the bottom two elements now in the top two positions of the new vector, and has therefore picked out the elements that we wanted. It has also shifted the elements of the vector up one (in some sense). As we will hopefully see soon, this is a handy property that we will use. Now, imagine that each element in the vector corresponds to a point in space. Also imagine that we want to calculate some function such that the new ``current point in space'' is equal to the difference between the points in space on either side of it. Using the vector we have at the moment, we therefore want the middle element ($v_2$) to become $v_1 - v_3$. We don't know what the other elements should become at this stage since there is no point ``outside'' them. So we want a vector output that looks something like this:

$$\begin{pmatrix}?\\ v_1 - v_3\\ ? \end{pmatrix}.$$ (3.19)

We could do this by using the following matrix:

$$\begin{pmatrix}0 & -1 & 0\\ 1 & 0 & -1\\ 0 & 1 & 0 \end{pmatrix}.$$ (3.20)

Let's see if it works...

$$\begin{pmatrix}0 & -1 & 0\\ 1 & 0 & -1\\ 0 & 1 & 0 \end{pmatrix... ...2\\ v_3 \end{pmatrix} = \begin{pmatrix}-v_2\\ v_1 - v_3\\ v_2 \end{pmatrix}.$$ (3.21)

Well, it did basically what we wanted. The central element is now equal to $v_1 - v_3$, but the top element is $-v_2$ and the bottom element is $v_2$--why is this? Well, imagine that there are zeroes above and below the vector. So, to calculate what the top element would be, we would have ``something like'' $0-v_2$ for the top element and ``something like'' $v_2-0$ for the bottom element. Hence we get the output we see in Equation (3.21).

Note that the process here is much like what we want to do with the diffusion equation, namely: add the left- and right-hand points of a vector and subtract twice the central point. Using the (oversimplified) $3 \times 3$ matrices above we would do something like this:

$$\begin{pmatrix}-2 & 1 & 0\\ 1 & -2 & 1\\ 0 & 1 & -2 \end{pmatri... ...= \begin{pmatrix}v_2 - 2 v_1\\ v_1 - 2 v_2 + v_3\\ v_2 - 2 v_3 \end{pmatrix}.$$ (3.22)

From the central element of the output vector we can see that this has implemented what we wanted to achieve: $v'_2 = v_1 + v_3 - 2 v_2$--look like the diffusion equation? Now we hopefully have a feel for the structure of the $\mbf{D}$ matrix for the diffusion equation.

Another way of looking at this is if we take the right-hand side of the discretised diffusion equation:

$$\kappa \frac{A_{i+1}^n + A_{i-1}^n - 2 A_i^n}{\Delta x^2},$$ (3.23)

then we can rewrite it as

$$\sum_{j=1}^N D_{ij} A_j^n,$$ (3.24)

where the $D_{ij}$ are the elements of some matrix $\mbf{D}$ whose elements can be written as

$$D_{i,j} = \frac{\kappa}{\Delta x^2} \bigl( \delta_{i+1,j}^n + \delta_{i-1,j}^n - 2 \delta_{i,j}^n \bigr),$$ (3.25)

where the $\delta_{i,j}$ are Kronecker delta functions and are equal to 1 if and only if $i = j$; zero otherwise. If we now think of the $A_i^n$ in the discretised diffusion equation as elements in a column vector, we can rewrite the diffusion equation in the following form:

$$\frac{\mbf{A}^{n+1} - \mbf{A}^n}{\Delta t} = \frac{\kappa}{\Delta x^2} \mbf{D}\mbf{A}^n.$$ (3.26)

We wish to find the value of the vector $\mbf{A}$ at the next point in time, which is at $t_{n+1}$, hence we solve this equation for $\mbf{A}^{n+1}$:

$$\mbf{A}^{n+1} = \mbf{A}^n + \frac{\kappa \Delta t}{\Delta x^2} \mbf{D}\mbf{A}^n,$$ (3.27)

$$= \Bigl( \mbf{I}+ \frac{\kappa \Delta t}{\Delta x^2} \mbf{D}\Bigr) \mbf{A}^n,$$ (3.28)

$$= \mbf{M}\mbf{A}^n,$$ (3.29)

where $\mbf{I}$ is the identity matrix and $\mbf{M}$ is just a matrix to make the equations look nice and simple. So, our differential equation problem has reduced to finding the $\mbf{D}$ matrix and then constructing the $\mbf{M}$ matrix from that and then repeatedly applying $\mbf{M}$ to evolve the vector $\mbf{A}$. This in its essence is all one has to code to use this representation.

Now, given that we have the boundary conditions that the function is constant at the left- and right-hand edges (Dirichlet boundary conditions) and are that they equal to zero, we construct the following $\mbf{D}$ matrix:

$$\mbf{D}= \begin{pmatrix}0 & 0 & 0 & \cdots & 0\\ 1 & -2 & 1 & \c... ...ts & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 0 \end{pmatrix}.$$ (3.30)

For an $\mbf{A}$ vector (at time instance $n$) such as

$$\mbf{A}^n = \begin{pmatrix}A_1^n\\ A_2^n\\ A_3^n\\ \vdots\\ A_N^n \end{pmatrix},$$ (3.31)

verify for yourself that one application of the matrix $\mbf{D}$ onto $\mbf{A}^n$ gives

$$\mbf{D}\mbf{A}^n = \begin{pmatrix}0\\ A_1^n - 2 A_2^n + A_3^n\\ ... ... A_3^n + A_4^n\\ \vdots\\ A_{N-2}^n - 2 A_{N-1}^n + A_N^n\\ 0 \end{pmatrix}.$$ (3.32)

Notice how the top and bottom rows of zeroes in $\mbf{D}$ ensure that the ``edge'' elements of $\mbf{A}$ stay zero. Thus the boundary conditions can be made implicit in the $\mbf{D}$ matrix. Handy eh?

What if we wanted periodic boundary conditions? Such a physical example would be if we wished to model the diffusion of a wave through time and space. So, if we had periodic boundary conditions, we could let the wave travel across the ``window'' in space we have defined (to the left, say) and it would appear from the right traveling to the left. This would allow us to let the wave ``see'' a longer ``distance'' without having to have really large windows in space. To do this with a variant of the $\mbf{D}$ matrix we would have to think carefully about the problem, and then realise how the $\mbf{D}$ matrix acts upon the $\mbf{A}$ vector in an element-by-element fashion. You may wish to verify that the $\mbf{D}$ matrix that would give diffusion equation type behaviour with periodic boundary conditions is

$$\mbf{D}= \begin{pmatrix}-2 & 1 & 0 & \cdots & 0 & 0 & 1\\ 1 & -2... ...& 0 & 0 & \cdots & 1 & -2 & 1\\ 1 & 0 & 0 & \cdots & 0 & 1 & -2 \end{pmatrix}.$$ (3.33)

Note the 1's in the top-right and bottom-left corners. They are like where the $\mbf{D}$ matrix has ``wrapped'' around from the other side of the matrix. It is this ``wrapped around'' kind of structure in the $\mbf{D}$ matrix that gives the periodic nature of the boundary conditions.

What are the advantages of using this matrix scheme? Well, your code is a lot easier to read, and a lot easier to change. The main bit of mental work on the programmer's part is in calculation of what the $\mbf{D}$ matrix should be to apply to the $\mbf{A}$ vector. Then, inside the main loop, all one has to do is have a line saying something like:

A = M*A;

to evolve the solution. It is literally that simple.

Ok. I hope that you have a bit of a grasp of the matrix representation of differential equations now. Using this form for the equations allows me to be somewhat less verbose in my explanations, and makes some of the mathematics a bit more obvious. We go on to investigate the implicit FTCS scheme in the next section.

Footnotes

... matrix^3.9

Why square? Well, we will get back a vector with the same dimensions as before after applying the matrix to it.