-
-
Save cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9 to your computer and use it in GitHub Desktop.
import numpy as np | |
## Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver | |
def TDMAsolver(a, b, c, d): | |
''' | |
TDMA solver, a b c d can be NumPy array type or Python list type. | |
refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm | |
and to http://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) | |
''' | |
nf = len(d) # number of equations | |
ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy arrays | |
for it in xrange(1, nf): | |
mc = ac[it-1]/bc[it-1] | |
bc[it] = bc[it] - mc*cc[it-1] | |
dc[it] = dc[it] - mc*dc[it-1] | |
xc = bc | |
xc[-1] = dc[-1]/bc[-1] | |
for il in xrange(nf-2, -1, -1): | |
xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il] | |
return xc |
Thanks for this! I've forked your gist to make it work in 3.0 and increase the speed by a factor of 2 by using numba.
License? (I also nagged ofan666)
(Thanks for posting this, by the way.)
I've added a GNU General Public License as a separate gist. Does this help? (thanks for checking by the way!)
Thanks! That helps.
How to be in case when on diagonal line during computing we get zero?
for example:
A = np.array([[6,2,0,0],[3,1,4,0],[0,1,7,5],[0,0,3,4]],dtype=float)
How to be in case when on diagonal line during computing we get zero?
for example:
A = np.array([[6,2,0,0],[3,1,4,0],[0,1,7,5],[0,0,3,4]],dtype=float)
@BatoNam, there is condition, when you could use this method. Conditional is
abs(a[i]) + abs(c[i]) <= abs(b[i]) for each i
xc is just a reference to bc
It's pointless.
Example: