Linear Algebra & Systems (scipy.linalg)

Almost every physical simulation—from finite element stiffness equations to heat conduction and circuit analysis—culminates in a linear system A x = b.

scipy.linalg wraps standard LAPACK/BLAS routines. It is more complete, robust, and often faster than numpy.linalg.

1. Golden Rule: Never Invert a Matrix Explicitly

In textbooks, you often see x = A^(-1) b. In production code, almost never call linalg.inv(A).

• Computing linalg.inv(A) requires about $2n^3$ operations, whereas LU factorization requires $2n^3 / 3$ (3x faster).
• Calculating an explicit inverse and multiplying by b degrades numerical accuracy and inflates residuals ||b - A x||.

Solving Systems and Comparing Residuals

2. Factor Once, Solve Many Times (lu_factor, lu_solve)

When the system matrix A remains fixed across multiple load cases or time steps, factorizing once with lu_factor reduces every subsequent solve to $O(n^2)$ back-substitution:

3. Symmetric Positive Definite Systems (assume_a="pos", Cholesky)

Stiffness and covariance matrices are Symmetric Positive Definite (SPD). Cholesky factorization (A = L L^T):

• Halves the floating-point operations compared to standard LU.
• Requires no row pivoting.
• Automatically raises LinAlgError if the matrix is not positive definite, giving a free physical validity check.

x = linalg.solve(A, b, assume_a="pos")

4. Eigenvalues, SVD & Matrix Exponentials

lam, Phi = linalg.eigh(K, M)  # Real, sorted eigenvalues

U, s, Vt = linalg.svd(A)
rank = np.linalg.matrix_rank(A)

from scipy.linalg import expm
x_t = expm(A * t) @ x0

Test eigenvalue decomposition on an oscillator system:

5. Try It Yourself

Problem: Solve a Banded Tridiagonal System

Solve the tridiagonal system A x = b where A has diagonal [2, 2, 2] and off-diagonals [-1, -1], with b = [1, 1, 1].

import numpy as np
from scipy import linalg
 
A = np.array([[2.0, -1.0, 0.0],
              [-1.0, 2.0, -1.0],
              [0.0, -1.0, 2.0]])
b = np.array([1.0, 1.0, 1.0])
 
x = linalg.solve(A, b, assume_a="pos")
print("x:", np.round(x, 4))