Direct solvers of system of algebraic equations.
Jump to
Gauss elimination
\begin{equation}
\begin{split}
a_{11}x_{1}+a_{12}x_{2}+.......................+ a_{1n}x_{n} = b1\\
a_{21}x_{1}+a_{22}x_{2}+.......................+ a_{2n}x_{n} = b2\\
.....................................................................\\
.....................................................................\\
a_{n-1,1}x_{1}+a_{n-1,2}x_{2}+....+ a_{n-1,n}x_{n} = bn-1\\
a_{n,1}x_{1}+a_{n,2}x_{2}+....................+ a_{n,n}x_{n} = bn
\end{split}
\end{equation}Step 1:Making the first column zero
R_{2} = R_{2} - \frac{a_{21}}{a_{11}} R_{1}\\
R_{3} = R_{3} - \frac{a_{31}}{a_{11}} R_{1}\\
...........................\\
R_{i} = R_{i} - \frac{a_{i1}}{a_{11}} R_{1}\\Step 2:Making the second column elements (entries below diagonal) zero
a_{ij} = a_{ij} - \frac{a_{i2}}{a_{22}} a_{2j}\\
In general, kth column :
a_{ij} = a_{ij} - \frac{a_{ik}}{a_{kk}} a_{kj}\\
Forward elimination:
\begin{split}
\textcolor{#d0d0d0}{the \,column \, that \,is \,making \,zero}\text{for} \quad \text{ } k = 1:n-1\\
\textcolor{#d0d0d0}{march \,in \, k^{th} \,row\,below \,diagonal}\text{ for } \text{ } i = k+1:n\\
\textcolor{#d0d0d0}{get \,the \,ratio}\text{ } R = \frac{a_{ik}}{a_{kk}}\\
\textcolor{#d0d0d0}{march \,in \,i^{th}\,row\,right \,to \,diagonal}\text{ }\text{ }for j = k+1:n+1\\
\text{ }\text{ }\text{ }\text{ } a_{ij} = a_{ij} - R.a_{kk}\\
end\\
end\\
end
\end{split}There are three loops in the forward elimination. Hence the order of computation is n^{3}
\begin{equation}
\begin{split}
a_{11}x_{1}+a_{12}x_{2}+.......................+ a_{1n}x_{n} = b_{1}\\
a_{22}x_{2}+.......................+ a_{2n}x_{n} = b_{2}\\
....................................................\\
.................................................\\
a_{n-2,n-2}x_{n-2}+ a_{n-2,n-1}x_{n-1}+ a_{n-2,n}x_{n}= b_{n-2}\\
a_{n-1,n-1}x_{n-1}+ a_{n-1,n}x_{n} = b_{n-1}\\
a_{n,n}x_{n} = b_{n}
\end{split}
\end{equation}Backward substitution
\begin{split}
x_{n} = \frac{b_{n}}{a_{n,n}}\\
x_{n-1} = \frac{b_{n-1} -a_{n-1,n} x_{n} }{a_{n-1,n-1}}\\
x_{n-1} = \frac{b_{n-2} - [a_{n-2,n} x_{n} +a_{n-2,n-1}x_{n-1}]}{a_{n-2,n-2}}\\
...........................................................\\
x_{n-i} = \frac{b_{n-i} - \sum_{j=0}^{i-1} a_{n-i,n-j} x_{n-j} }{a_{n-i,n-i}}
\end{split}
Algorithm for backward substitution:
\begin{split}
x_{n} = \frac{b_{n}}{a_{n,n}}\\
\text{ for } \text{ } i = 1:n-1\\
sum =0\\
\text{ for } \text{ } j = 0:i-1\\
sum = sum + a_{n-i,n-j} x_{n-j} {a_{n-i,n-i}}\\
end\\
x_{n-i} = \frac{b_{n-i} - sum}{a_{n-i,n-j} }\\
end
\end{split}Here, the order of computation is n^{2}. So, when the matrix is triangular, the order reduces from n^{3} to n^{2}. This motivate to use triangular matrix instead of the full matrix.
If A is not diagonally dominant, pivotization (interchanging rows) shall be employed to make so. Else, the Gauss elimination would fail as there is a division by diagonal element which shall blow up.
LU decomposition
Taking a [4 x 4] matrix as example, LU decomposition is given by
\begin{bmatrix}
a_{11} & a_{12} & a_{13} &a_{14}\\
a_{21} & a_{22} & a_{23} &a_{24}\\
a_{31} & a_{32} & a_{33} &a_{24}\\
a_{41} & a_{42} & a_{43} &a_{24}
\end{bmatrix}
=
\begin{bmatrix}
l_{11} & 0 & 0 &0\\
l_{21} & l_{22} & 0 &0\\
l_{31} & l_{32} & l_{33} &0\\
l_{41} & l_{42} & l_{43} &l_{24}
\end{bmatrix}
\begin{bmatrix}
1 & u_{12} & u_{13} &u_{14}\\
0 & 1 & u_{23} &u_{24}\\
0 & 0& 1 &u_{24}\\
0 & 0 & 0 &1
\end{bmatrix}By comparing,
l_{i1} = a_{i1},
u_{1i} = \frac{a_{1i}}{l_{11}}\\..........................\begin{array}{ccccc}
AX = b & \to & \text{LU Decomposition}& \to &LUX = b &\text{[o(n}^{3}\text{)]}\\
L\underbrace{UX}_{Z}= b &\to & \text{forward substitution} &\to &Z = L^{-1}b &\text{[o(n}^{2}\text{)]}\\
UX= Z & \to & \text{back substitution} & \to & X = U^{-1}Z &\text{[o(n}^{2}\text{)]}\\
\end{array}LU decomposition also has the order of n^{3} hence it is equally as intensive as Gauss elimination. But, if A is symmetric (A^{T}=A) and positive definite (A^{T}vA > 0) computation of LU factorisation reduces by a factor of \frac{1}{2}.
The advantage of LU decomposition is when A is constant, but b changes, the factorisation need not be carried out again, making the overall order n^{2}
TDMA
\begin{bmatrix}
a_{1} & b_{1} & 0 &0&0&0\\
c_{1} & a_{2} & b_{2} &0&0&0\\
0 & c_{2} & a_{3} &b_{3}&0&0\\
0 & 0 & c_{3} &a_{4}&.&0\\
0 & 0 & 0 &.&.&b_{n-1}\\
0 & 0 &0 &0& c_{n-1}&a_{n}\\
\end{bmatrix}
\begin{bmatrix}
\phi_{1} \\
\phi_{2} \\
\phi_{3} \\
\phi_{4} \\
. \\
\phi_{n}\\
\end{bmatrix}
=
\begin{bmatrix}
d_{1} \\
d_{2} \\
d_{3} \\
d_{4} \\
. \\
d_{n}\\
\end{bmatrix}\begin{equation}
a_{i} \phi_{i} + b_{i} \phi_{i+1}+c_{i} \phi_{i-1} = d_{i}
\end{equation}a_{1} \phi_{1} = -b_{1} \phi_{2}+d_{1}\\
\phi_{1} = \frac{-b_{1} }{a_{1} }\phi_{2}+\frac{d_{1} }{a_{1} }\\
\text{ }\\
a_{2} \phi_{2} = -b_{2} \phi_{3}-c_{2} \phi_{1}+d_{2}\\
\text{ }\\
a_{2} \phi_{2} = -b_{2} \phi_{3}-c_{2} (\frac{-b_{1} }{a_{1} }\phi_{2}+\frac{d_{1} }{a_{1} })+d_{2}\\
........\\
\begin{equation}
\phi_{i} = P_{i}\phi_{i+1}+Q_{i}
\end{equation}for i = i-1 :
\begin{equation}
\phi_{i-1} = P_{i-1}\phi_{i}+Q_{i-1}
\end{equation}substituting (5) in (3),
a_{i} \phi_{i} + b_{i} \phi_{i+1} +c_{i} \phi_{i-1}=d_{i}\\
a_{i} \phi_{i} + b_{i} \phi_{i+1} +c_{i} (P_{i-1}\phi_{i}+Q_{i-1})=d_{i}\\
\text{ }\\
\phi_{i} = \frac{-b_{i} }{a_{i}+c_{i}P_{i-1}}\phi_{i+1}+ \frac{d_{i}+c_{i}Q_{i-1} }{a_{i}+c_{i}P_{i-1}}\\
\text{ }\\
\text{ but we have}\\
\phi_{i} = P_{i}\phi_{i+1}+Q_{i}\\
\text{ hence, comparing both equations :}\\
\begin{equation}
P_{i} = \frac{-b_{i} }{a_{i}+c_{i}P_{i-1}} \text{ and } Q_{i} = \frac{d_{i}-c_{i}Q_{i-1} }{a_{i}+c_{i}P_{i-1}}
\end{equation}atex]c_{1}=0[/katex]:
\begin{equation}
P_{1} = -\frac{b_{1} }{a_{1}} \text{ and } Q_{1} = \frac{d_{1} }{a_{1}}
\end{equation}P_{2} , Q_{2} to P_{n} , Q_{n} can be found using (6). P_{n} will be zero since b_{n} = 0 hence , \phi_{n} = Q_{n} . (\phi_{i} = P_{i}\phi_{i+1}+Q_{i}) Then using (4), all \phi_{i} can be calculated(back substitution). The algorithm is given below, which is the order of n^{3}
Algorithm
- Compute P_{1} and Q_{1}using Eq. (7)
- For i = 2, 3,…, N use forward recursion to compute the values of P_{i} and Q_{i} using (6)
- Set \phi{N} = Q_{N}
- For i = N-1, N-2,…3,2,1 use backward recursion (4) to get \phi{i}
Python Code
def tdma(a,b,c,d,T):
N=len(a)
P=np.ones(N)
Q=np.ones(N)
P[0]=-b[0]/a[0]
Q[0]= d[0]/a[0]
for i in range(1,N):
denom=1/(a[i]+c[i]*P[i-1])
P[i] = -b[i]*denom
Q[i] = (d[i]-c[i]*Q[i-1])*denom
T[N-1]=Q[N-1]
for i in range(N-2,-1,-1):
T[i]=P[i]*T[i+1]+Q[i]