Re-wrote Thomas Algorithm chapter to include more details

Author: jiegillet

contents/thomas_algorithm/thomas_algorithm.md

@@ -1,6 +1,6 @@
 # Thomas Algorithm
 
-As alluded to in the [Gaussian Elimination chapter](../gaussian_elimination/gaussian_elimination.md), the Thomas Algorithm (or TDMA -- Tri-Diagonal Matrix Algorithm) allows for programmers to **massively** cut the computational cost of their code from $$\sim O(n^3) \rightarrow \sim O(n)$$! This is done by exploiting a particular case of Gaussian Elimination, particularly the case where our matrix looks like:
+As alluded to in the [Gaussian Elimination chapter](../gaussian_elimination/gaussian_elimination.md), the Thomas Algorithm (or TDMA, Tri-Diagonal Matrix Algorithm) allows for programmers to **massively** cut the computational cost of their code from $$ O(n^3)$$ to $$O(n)$$ in certain cases! This is done by exploiting a particular case of Gaussian Elimination, in the case where our matrix looks like:
 
 $$
 \left[
@@ -14,27 +14,83 @@ $$
 \right]
 $$
 
-By this, I mean that our matrix is *Tri-Diagonal* (excluding the right-hand side of our system of equations, of course!). Now, at first, it might not be obvious how this helps; however, we may divide this array into separate vectors corresponding to $$a$$, $$b$$, $$c$$, and $$d$$ and then solve for $$x$$ with back-substitution, like before.
+This matrix shape is called *Tri-Diagonal* (excluding the right-hand side of our system of equations, of course!).
+Now, at first, it might not be obvious how this helps. Well, first of all, the system easier to encode: we may divide it into four separate vectors corresponding to $$a$$, $$b$$, $$c$$, and $$d$$ (in some implementation, you will see the missing $$a_0$$ and $$c_n$$ set to zero to get four vectors of the same size).
+Secondly, equations this short are easy to solve analytically, let's see how.
 
-In particular, we need to find an optimal scale factor for each row and use that. What is the scale factor? Well, it is the diagonal $$-$$ the multiplicative sum of the off-diagonal elements.
-In the end, we will update $$c$$ and $$d$$ to be $$c'$$ and $$d'$$ like so:
+We'll apply the mechanisms introduced in the Gaussian Elimination chapter.
+Our first goal is to eliminate the $$a_i$$ terms and set the diagonal values $$b_i$$ to $$1$$. The $$c_i$$ and $$d_i$$ terms will be transformed into $$c'_i$$ and $$d'_i$$.
+The first row is particularly easy to transform since there is no $$a_0$$, we simply need to divide the row by $$b_0$$:
 
 $$
+\left\{
 \begin{align}
-c'_i &= \frac{c_i}{b_i - a_i \times c'_{i-1}} \\
-d'_i &= \frac{d_i - a_i \times d'_{i-1}}{b_i - a_i \times c'_{i-1}}
+c'_0 &= \frac{c_0}{b_0} \\
+d'_0 &= \frac{d_0}{b_0}
 \end{align}
+\right.
+$$
+
+Let's assume that we found a way to transform the first $$i-1$$ rows. How would we transform the next one? We have
+
+$$
+\begin{array}{ccccccc|c}
+    &  & \ddots & & & & &  \\
+(i-1) & & 0 & 1 & c'_{i-1} & & & d'_{i-1} \\
+(i)   & &   & a_i & b_i & c_i & & d_i \\
+      & &   &   &   &  \ddots &  &
+\end{array}
 $$
 
-Of course, the initial elements will need to be specifically defined as
+Let's transform row $$(i)$$ in two steps.
+**Step one**: eliminate $$a_i$$ with the transformation $$(i)^* = (i) - a_i \times (i-1)$$:
 
 $$
+\left\{
 \begin{align}
-c'_0 = \frac{c_0}{b_0}
-d'_0 = \frac{d_0}{b_0}
+a^*_i &= 0 \\
+b^*_i &= b_i - a_i \times c'_{i-1} \\
+c^*_i &= c_i \\
+d^*_i &= d_i - a_i \times d'_{i-1}
 \end{align}
+\right.
 $$
 
+**Step two**: get $$b'_i=1$$ with the transformation $$(i)' = (i)^* / b^*_i $$:
+
+$$
+\left\{
+\begin{align}
+a'_i &= 0 \\
+b'_i &= 1 \\
+c'_i &= \frac{c_i}{b_i - a_i \times c'_{i-1}} \\
+d'_i &= \frac{d_i - a_i \times d'_{i-1}}{b_i - a_i \times c'_{i-1}}
+\end{align}
+\right.
+$$
+
+Brilliant! With the last two formula, we can calculate all the $$c'_i$$ and $$d'_i$$ in a single pass, starting from row $$1$$, since we already know the values of $$c'_0$$ and $$d'_0$$.
+
+Of course, what we really need are the solutions $$x_i$$. It's back substitution time!
+
+If we express our system in terms of equations instead of a matrix, we get
+
+$$
+x_i + c'_i \times x_{i+1} = d'_i
+$$
+
+plus the last row that is even simpler: $$x_n = d'_n$$. One solution for free!
+Maybe we can backtrack from the last solution? Let's (barely) transform the above equation:
+
+$$
+x_i = d'_i - c'_i \times x_{i+1}
+$$
+
+and that's all there is to it. We can calculate all the $$x_i$$ in a single pass starting from the end.
+
+Overall, we only need two passes, and that's why our algorithm is $$O(n)$$!
+The transformations are quite easy too, isn't that neat?
+
 ## Example Code
 
 {% method %}
@@ -48,27 +104,6 @@ $$
 [import, lang:"haskell"](code/haskell/thomas.hs)
 {% endmethod %}
 
-This is a much simpler implementation than Gaussian Elimination and only has one for loop before back-substitution, which is why it has a better complexity case.
-
 <script>
 MathJax.Hub.Queue(["Typeset",MathJax.Hub]);
 </script>
-$$
-\newcommand{\d}{\mathrm{d}}
-\newcommand{\bff}{\boldsymbol{f}}
-\newcommand{\bfg}{\boldsymbol{g}}
-\newcommand{\bfp}{\boldsymbol{p}}
-\newcommand{\bfq}{\boldsymbol{q}}
-\newcommand{\bfx}{\boldsymbol{x}}
-\newcommand{\bfu}{\boldsymbol{u}}
-\newcommand{\bfv}{\boldsymbol{v}}
-\newcommand{\bfA}{\boldsymbol{A}}
-\newcommand{\bfB}{\boldsymbol{B}}
-\newcommand{\bfC}{\boldsymbol{C}}
-\newcommand{\bfM}{\boldsymbol{M}}
-\newcommand{\bfJ}{\boldsymbol{J}}
-\newcommand{\bfR}{\boldsymbol{R}}
-\newcommand{\bfT}{\boldsymbol{T}}
-\newcommand{\bfomega}{\boldsymbol{\omega}}
-\newcommand{\bftau}{\boldsymbol{\tau}}
-$$