John Travis
Mississippi College
When solving $Ax=b$ for a general system of equations, one often might want to use Gaussian Elimination with scaled partial pivoting. However, one might want to solve such a system several times with a constant matrix $A$ but with a variety of "forcing term" vectors for b. Using the standard approach would replicate a significant number of calculations in the reduction of A to triangular form. There is however a way to save the intermediate steps of Gaussian Elimination so that only a final forward and backward solve would be necessary and not all of the elimination steps.
In general, one might be able to factor $A = LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. Then to solve $Ax=b$ simply replace and solve
$LUx = b$
using a twostep approach:
Solve $Ly=b$ for $y$ using forward solving.
Solve $Ux=y$ for $x$ using backward solving.
Since each of these require only $O(n^2)$ calculations, then repeated solves with various values for $b$ offer significant computational savings.
A = 2 1 3 0 1 1 2 3 0 0 4 1 5 0 0 4 The Factorization: L = 2 0 0 0 1 0.5 0 0 0 0 4 0 5 2.5 5 20.25 U = 1 0.5 1.5 0 0 1 1 6 0 0 1 0.25 0 0 0 1 and multiplying these together gives the original A: A = 2 1 3 0 1 1 2 3 0 0 4 1 5 0 0 4 A = 2 1 3 0 1 1 2 3 0 0 4 1 5 0 0 4 The Factorization: L = 2 0 0 0 1 0.5 0 0 0 0 4 0 5 2.5 5 20.25 U = 1 0.5 1.5 0 0 1 1 6 0 0 1 0.25 0 0 0 1 and multiplying these together gives the original A: A = 2 1 3 0 1 1 2 3 0 0 4 1 5 0 0 4 
b = 8 2 7 4 Note the temporary values in the vector y: y = 4 4 1.75 0.851852 With the final answer for Ax=b as x equals: x = 1.48148 0.851852 1.96296 0.851852 and to check, look at bAx: ans = 0 0 0 8.88178e16 b = 8 2 7 4 Note the temporary values in the vector y: y = 4 4 1.75 0.851852 With the final answer for Ax=b as x equals: x = 1.48148 0.851852 1.96296 0.851852 and to check, look at bAx: ans = 0 0 0 8.88178e16 

Solving using b = 8 2 7 4 y = 4 4 1.75 0.851852 yields x = 1.48148 0.851852 1.96296 0.851852 with bAx = ans = 0 0 0 8.88178e16 Solving using b = 9 2 7 5 y = 4.5 13 1.75 2.53086 yields x = 1.02469 0.197531 2.38272 2.53086 with bAx = ans = 0 0 0 7.10543e15 Solving using b = 8 2 7 4 y = 4 4 1.75 0.851852 yields x = 1.48148 0.851852 1.96296 0.851852 with bAx = ans = 0 0 0 8.88178e16 Solving using b = 9 2 7 5 y = 4.5 13 1.75 2.53086 yields x = 1.02469 0.197531 2.38272 2.53086 with bAx = ans = 0 0 0 7.10543e15 
invA = 0.197531 0.197531 0.0493827 0.160494 0.419753 0.580247 0.604938 0.283951 0.0617284 0.0617284 0.234568 0.0123457 0.246914 0.246914 0.0617284 0.0493827 Check our result by comparing to the builtin function for inverses: ans = 1.1102e16 8.3267e17 0.0000e+00 0.0000e+00 1.1102e16 1.1102e16 0.0000e+00 5.5511e17 6.9389e18 6.9389e18 2.7756e17 0.0000e+00 2.7756e17 2.7756e17 6.9389e18 0.0000e+00 invA = 0.197531 0.197531 0.0493827 0.160494 0.419753 0.580247 0.604938 0.283951 0.0617284 0.0617284 0.234568 0.0123457 0.246914 0.246914 0.0617284 0.0493827 Check our result by comparing to the builtin function for inverses: ans = 1.1102e16 8.3267e17 0.0000e+00 0.0000e+00 1.1102e16 1.1102e16 0.0000e+00 5.5511e17 6.9389e18 6.9389e18 2.7756e17 0.0000e+00 2.7756e17 2.7756e17 6.9389e18 0.0000e+00 
