Using the LU Matrix for Multiple Right-Hand Sides

• Many physical situations are modelled with a large set of linear equations.
• The equations will depend on the geometry and certain external factors that will determine the right-hand sides.
• For example, in electrical circuit problems, the resistors at the circuit ($A$ matrix) are unchanged with the varying applied voltages ($b$ vector). (e.g., Kirchhoff's Rule)
• If we want the solution for many different values of these right-hand sides,
  • it is inefficient to solve the system from the start with each one of the right-hand-side values.
  • Using the $LU$ equivalent of the coefficient matrix is preferred.
• Suppose we have solved the system $Ax=b$ by Gaussian elimination.
• We now know the $LU$ equivalent of $A$: $A =L *U$

• We can write
  $$Ax=b$$
  $$LUx=b$$
  $$Ly=b$$

• e.g., Solve $Ax=b$, where we already have its $L$ and $U$ matrices:

  

  

• Suppose that the $b$-vector is $[6~ -7~ -2~ 0]^T$.
• We first get $y(=Ux)$ from $Ly = b$ by forward substitution:
  $$y = [6~ -11~-9~ -3.12]^T$$

• and use it to compute $x$ from $Ux = y$:
  $$x= [ -0.5~ 1~ 0.3333~ -2]^T.$$

• Exercise: $b = [1~ 4~ -3~ 1]^T$ $\Longrightarrow$ $x=[0.0128~ -0.5897~ -2.0684~ 2.1795]^T$

Phyton Code:

    1 import numpy as np
    2 from scipy.linalg import lu
    3 A = np.array([[0.0, 2.0, 0.0, 1.0], [2.0, 2.0, 3.0, 2.0], [4.0, -3.0, 0.0, 1.0], [6.0, 1.0, -6.0, -5.0]])
    4 P, L, U = lu(A)
    5 print("SciPy LU-decomposition: P - Permutation Matrix \n", P)
    6 print("SciPy LU-decomposition: L - Lower Triangular with unit diagonal elements \n", L)
    7 print("SciPy LU-decomposition: U - Upper Triangular \n", U)
    8 def forward(L, b):
    9     y=np.zeros(np.shape(b),dtype=float)
   10     for i in range(len(b)):
   11         y[i]=np.copy(b[i])
   12         for j in range(i):
   13             y[i]=y[i]-(L[i, j]*y[j])
   14         y[i] = y[i]/L[i, i]
   15     return y
   16 b = np.array([[6.0], [-7.0], [-2.0], [0.0]])
   17 # b = np.array([[1.0], [4.0], [-3.0], [1.0]])
   18 y=forward(L,b)
   19 print("y vector from Ly=b by forward substitution :", np.transpose(y))
   20 def backward(U, y):
   21     x=np.zeros(np.shape(y),dtype=float)
   22     ylen=len(y)-1
   23     x[ylen] =y[ylen]/U[ylen, ylen] # Print the last stage x value
   24     for i in range(ylen-1,-1,-1):
   25         x[i]=np.copy(y[i])
   26         for j in range(ylen,i,-1):
   27             x[i]=x[i]-(U[i, j]*x[j])
   28         x[i] = x[i]/U[i, i]
   29     return x
   30 x=backward(U,y)
   31 print("x vector from Ux=y by backward substitution :", np.transpose(x))