Although the numerical analyst or the engineer is able to write his own code, there’s no doubt that it’s safer and easier to use a code already done. And if this code is part of a well conceptualized library, there’s no reason to reinvent the wheel.

To exemplify this idea, let’s solve the tridiagonal system shown in the post TDMA at glance using the SciPy library.

See the Python code:

import numpy as np
from scipy.linalg import solve_banded

def SolveBanded(A, D):
  # Find the diagonals
  ud = np.insert(np.diag(A,1), 0, 0) # upper diagonal
  d = np.diag(A) # main diagonal
  ld = np.insert(np.diag(A,-1), len(d)-1, 0) # lower diagonal

  # simplified matrix
  ab = np.matrix([
    ud,
    d,
    ld,
  ])

  return solve_banded((1, 1), ab, D )

if __name__ == "__main__":
  A = np.array([
     [20, -5, 0, 0],
     [-5, 15, -5, 0],
     [ 0, -5, 15, -5],
     [ 0, 0, -5, 10]
])
D = np.array([1100, 100, 100, 100])
print "SciPy - Solve Banded"
print "A:", A
print "D:", D
print "Result:", SolveBanded(A, D)

As you can see, this code is very simple. This simplicity allow us to focus our attention on what realy matters when we are solving numerical problems, since we are (almost) sure that an eventual problem is not in the library.

References

  • SymPy Development Team (2016). SymPy: Python library for symbolic mathematics
    http://www.sympy.org