From 24d26989c98d1f349bee26edf4fa308be3727b5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9mie=20Gillet?= Date: Sat, 4 Aug 2018 11:26:51 +0900 Subject: [PATCH 1/2] Re-wrote Thomas Algorithm chapter to include more details --- contents/thomas_algorithm/thomas_algorithm.md | 95 +++++++++++++------ 1 file changed, 65 insertions(+), 30 deletions(-) diff --git a/contents/thomas_algorithm/thomas_algorithm.md b/contents/thomas_algorithm/thomas_algorithm.md index 2be5b674e..3e17e3759 100644 --- a/contents/thomas_algorithm/thomas_algorithm.md +++ b/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. - -$$ -\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}} -$$ From 08d9a20ec297968d57c05bbfda6d3658d9d7e443 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9mie=20Gillet?= Date: Sun, 5 Aug 2018 22:08:47 +0900 Subject: [PATCH 2/2] Included Leios' suggestions --- contents/thomas_algorithm/thomas_algorithm.md | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/contents/thomas_algorithm/thomas_algorithm.md b/contents/thomas_algorithm/thomas_algorithm.md index 3e17e3759..0ec058de6 100644 --- a/contents/thomas_algorithm/thomas_algorithm.md +++ b/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 $$ 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: +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[ @@ -15,10 +16,10 @@ $$ $$ 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. +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. -We'll apply the mechanisms introduced in the Gaussian Elimination chapter. +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$$: @@ -43,6 +44,7 @@ $$ $$ Let's transform row $$(i)$$ in two steps. + **Step one**: eliminate $$a_i$$ with the transformation $$(i)^* = (i) - a_i \times (i-1)$$: $$