3-2 TDMA （Thomas’algorithm）[python]
What is TDMA？
Solve an equation whose coefficient matrix is tridiagonal
$$\left[\begin{array}{cccccc}b_0&c_0&&&& \\ a_0&b_0&c_1&&& \\ &\ddots& \ddots & \ddots && \\&& \ddots & \ddots & \ddots &\\ &&&a_{n-3}&b_{n-2}&c_{n-2} \\ &&&&a_{n-2} &b_{n-1} \end{array}\right] \left[\begin{array}{c} u_0 \\ u_1 \\ \vdots \\ \vdots \\u_{n-2} \\ u_{n-1} \end{array}\right] = \left[\begin{array}{c} d_0 \\ d_1 \\ \vdots \\ \vdots \\d_{n-2} \\ d_{n-1} \end{array}\right] $$
Consider a tridiagonal matrix of n × n matrix as above. We want to find the u matrix.
Note that the number of b and d are n each, and a and c are n-1 each.
① First, the formula of the first row of the matrix can be written and transformed as follows.
$$b_{0}c_{0}+c_{0}u_{1}=d_1 ⇔ u_0+\frac{c_0}{b_0}u_1=\frac{d_1}{b_0}$$
By rearranging \(e_0=\frac{c_0}{b_0} ,f_0=\frac{d_1}{b_0}\) the following formula is obtained.
$$u_0=f_0-e_{0}u_{1}$$
② Substitute the result of the previous row into the next row.
$$a_{0}u_{0}+b_{1}u_{1}+c_{1}u_{2}=d_1$$
$$⇔ a_{0}(f_{0}-e_{0}u_{1})+b_{1}u_{1}+c_{1}u_{2}=d_1$$
$$⇔ (b_{1}-a_{0}e_{0})u_{1}+c_{1}u_{2}=d_1-a_{0}f_{0}$$
$$⇔ u_{1}+\frac{c_{1}}{b_{1}-a_{0}e_{0}}u_{2}=\frac{d_{1}-a_{0}f_{0}}{b_{1}-a_{0}e_{0}}$$
By rearranging\(e_{1}=\frac{c_{1}}{b_{1}-a_{0}e_{0}} ,f_1=\frac{d_{1}-a_{0}f_{0}}{b_{1}-a_{0}e_{0}}\), the following formula is obtained.
$$u_{1}=f_{1}-e_{2}u_{2}$$
③ By repeating the operation of ② in the same way, We can obtain the following formula in the second row from the bottom of the matrix.
$$u_{n-2}=f_{n-2}-e_{n-2}u_{n-1}$$
Also, the formula in the bottom row of the matrix is as follows.
$$a_{n-2}u_{n-2}+b_{n-1}u_{n-1}=d_{n-1}$$
By combining these two equations, \(u_{n-1}\)can be obtained as follows.
$$a_{n-2}(f_{n-2}-e_{n-2}u_{n-1} )+b_{n-1}u_{n-1}=d_{n-1}$$
$$⇔ (b_{n-1}-a_{n-2}e_{n-2})u_{n-1}=d_{n-1}-a_{n-2}f_{n-2}$$
$$⇔ u_{n-1}=\frac{d_{n-1}-a_{n-2}f_{n-2}}{b_{n-1}-a_{n-2}e_{n-2}}$$
④ Finally, by repeatedly assigning n-2 to 0 in the following formula created in ②, all\(u\)can be obtained.
$$u_{i}=f_{i}-e_{i}u_{i+1} \hspace{12pt} (i=n-2,n-1,\cdots 0)$$
Summery
After all, this method performs forward elimination and back substitution like Gaussian elimination. However, the calculation amount is small because only the tridiagonal part is calculated.
We will test TDMA, using python. Please, consider the following matrix.
$$\left[\begin{array}{ccc}3&1&0 \\ 1&4&2 \\ 0&2&5 \end{array}\right] \bf u = \left[\begin{array}{c} 5 \\ 15 \\ 19 \end{array}\right] $$
If we get\(u= \left[\begin{array}{c} 1 \\ 2 \\ 3 \end{array}\right] \), it is successful.
We’ve calculated with the python code which is at the end of this article. The result is as follows.
[1.0, 2.0000000000000004, 2.9999999999999996]
Good results!
The order of calculation is from \(u_{n-1}\) to\(u_{0}\)so I wrote 『u.reverse()』on the 28th lines to reverse \(u\).
Next time, I will try to solve the steady-state advection-diffusion equation using this method.
by hide
import numpy
def TDMA(a, b, c, d):
"""
a,cはn-1個の要素
b,dはn個の要素
"""
# ①
e = /b[0]]
f = [d[0]/b[0]]
# ②
for i in range(1, len(d)-1):
e.append(c[i]/(b[i]-a[i-1]*e[i-1]))
f.append((d[i]-a[i-1]*f[i-1])/(b[i]-a[0]*e[0]))
# ③
u = [(d[len(d)-1]-a[len(d)-2]*f[len(d)-2]) /
(b[len(d)-1]-a[len(d)-2]*e[len(d)-2])]
# ④
for i in range(len(d)-2, -1, -1):
u.append(f[i]-e[i]*u[len(d)-2-i])
u.reverse()
return u
#-----------------------------
#-- main --
#-----------------------------
a = [1.0, 2.0]
b = [3.0, 4.0, 5.0]
c = [1.0, 2.0]
d = [5.0, 15.0, 19.0]
kai = TDMA(a, b, c, d)
print(kai)
Tweet
Leave a Reply