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.
• 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$.
• Consider now that we want to solve $Ax=b$ with some new $b$-vector. We can write
  $$Ax=b=L*U*x=b$$

• The product of $U$ and $x$ is a vector, call it $y$. Now, we can solve for $y$ from $Ly = b$ and this is readily done because $L$ is lower-triangular and we get $y$ by forward-substitution. Call the solution $y = b'$.
• Going back to the original $LUx = b$, we see that, from $Ux=y=b'$, we can get $x$ from $Ux = b'$, which is again readily done by back-substitution ($U$ is upper-triangular).
• i.e., Solve $Ax=b$, where we already have its $L$ and $U$ matrices:
  $$\left[ \begin{array}{rrrr} 1&0&0&0\\ 0.66667 & 1 & 0 & 0 ... ... & 6.8182 &5.6364\\ 0& 0 &0 & 1.5600 \\ \end{array} \right]$$
  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.$$

• Now, if we want the solution with a different $b$-vector;
  $$bb = [1 4 -3 1]^T$$
  we just do $Ly = bb$ to get
  $$y = [1, 3.3333,-1.8182,3.4]^T$$
  and then use this $y$ in $Ux = y$ to find the new $x$:
  $$x=[0.0l28, -0.5897, -2.0684, 2.1795]^T$$