Re-wrote Thomas Algorithm chapter to include more details (#329)

* Re-wrote Thomas Algorithm chapter to include more details

* Included Leios' suggestions

Author: jiegillet

contents/thomas_algorithm/thomas_algorithm.md

@@ -1,6 +1,7 @@
 # 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 where the matrix looks like this:
 
 $$
 \left[
@@ -14,27 +15,84 @@ $$
 \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, firstly, it makes the system easier to encode: we may divide it into four separate vectors corresponding to $$a$$, $$b$$, $$c$$, and $$d$$ (in some implementations, you will see the missing $$a_0$$ and $$c_n$$ set to zero to get four vectors of the same size).
+Secondly, and most importantly, equations this short and regular are easy to solve analytically.
 
-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 start by applying mechanisms familiar to those who have read the [Gaussian Elimination](../gaussian_elimination/gaussian_elimination.md) 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}
+$$
+
+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}
+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.
 $$
 
-Of course, the initial elements will need to be specifically defined as
+**Step two**: get $$b'_i=1$$ with the transformation $$(i)' = (i)^* / b^*_i $$:
 
 $$
+\left\{
 \begin{align}
-c'_0 = \frac{c_0}{b_0}
-d'_0 = \frac{d_0}{b_0}
+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,8 +106,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>