Solving Linear Systems
The LinearAlgebra package not only gives you tools to solve linear systems directly, but also allows you to examine each step of the process.
By default, Maple uses hardware floats, but in this worksheet, we will use software floats. (For information about hardware floats and using them effectively, see the LinearAlgebra Functions with NAG Support example worksheet.)
>


>


>



Solving with a Single Command  LU Method


The following Matrix and Vector represent a linear system M.x=V which we wish to solve for x:
>


>


 (1.1) 
To solve the linear system, you can use the LinearSolve command. Note that you can use separate Matrix and righthand side Vector arguments, or simply use an augmented Matrix. To ensure that the method used is LU decomposition (the stepbystep method demonstrated in the next section), specify method='LU'.
>


 (1.2) 
Compute the residual to see how accurate the solution is:
>


 (1.3) 


Solving the Linear System One Step at a Time  LU Method


You can solve a linear system one step at time by using the LU Method.
As in the previous section, the following Matrix and Vector represent a linear system M.x=V where we want to solve for x:
>


>


 (2.1) 
First, do an LU decomposition of M. The output components are P (the permutation Matrix), L (a unit lower triangular Matrix), and U (an upper triangular Matrix).
>


 (2.2) 
Since this is a decomposition of M, P.L.U=M. The inverse of a permutation Matrix is always its transpose. To solve the system, M.x=V, start by multiplying both sides by the inverse (transpose) of P.
>


 (2.3) 
We now have L.U.x = Transpose(P).V, or L.U.x=V2 for short.
Since L is lower triangular, use ForwardSubstitute to find a Vector V3 such that L.V3 = V2:
>


 (2.4) 
Since L.V3=V2 and L.U.x=V2, we need to find a solution vector, x such U.x=V3. Since U is upper triangular, use BackwardSubstitute to do this:
>


 (2.5) 
Compute the residual to see how accurate the solution is:
>


 (2.6) 


Solving with a Single Command  QR Method


The following Matrix and Vector represent a linear system M.x=V where we want to solve for x:
>


>


 (3.1) 
To solve the linear system, you can use the LinearSolve command. To ensure that the method used is QR decomposition (the stepbystep method demonstrated in the next section), specify method='QR'.
>


 (3.2) 
Compute the residual to see how accurate the solution is:
>


 (3.3) 


Solving the Linear System One Step at a Time  QR Method


The QR method can be used to solve the system one step at a time.
As in the previous section, the following Matrix and Vector represent a linear system M.x=V where we want to solve for x:
>


>


 (4.1) 
First, do a QR decomposition of M. The output components are Q (an orthogonal Matrix) and R (an upper triangular Matrix).
>


 (4.2) 
Since this is a decomposition of M, Q.R=M. The inverse of an orthogonal Matrix is always its transpose (its Hermitian transpose if complex values are involved). To solve the system, M.x=V, start by multiplying both sides by the inverse (transpose) of Q.
>


 (4.3) 
We now have R.x = Transpose(Q).V, or R.x=V2 for short.
Since R is upper triangular use BackwardSubstitute to solve for x:
>


 (4.4) 
Compute the residual to see how accurate the solution is:
>


 (4.5) 


Cholesky Method


The Cholesky method uses a special form of LU decomposition. It works only for symmetric (Hermitian in the case of complex entries), positive definite Matrices. When you have a Matrix of this form, the Cholesky method is faster than the LU method because it does not involve a permutation Matrix, and it uses less storage because the upper triangular Matrix is the transpose (Hermitian transpose in the case of complex entries) of the lower triangular Matrix, and thus only one of them needs to be stored.
To try this method, you will need a symmetric, positive definite Matrix. One way to create such a Matrix is to create a lower triangular Matrix, then multiply it by its transpose. Here is an example, followed by verification that it is indeed symmetric and positive definite:
>


 (5.1) 
>


 (5.2) 
>


 (5.3) 
You will also need a righthand side vector:
>


 (5.4) 

Solving with a Single Command


The system can be solved with a single call to LinearSolve. Since the Matrix has the properties required by the Cholesky method, that method will automatically be used. An explicit request for this method can also be made.
>


 (5.1.1) 
>


 (5.1.2) 
Compute the residual to see how accurate the solution is:
>


 (5.1.3) 


Solving StepByStep


To solve the system one step at a time, start with a Cholesky decomposition:
>


 (5.2.1) 
Use ForwardSubstitute to solve L.V2 = V:
>


 (5.2.2) 
Use BackwardSubstitute to solve Transpose(L).x=V2:
>


 (5.2.3) 
Compute the residual to see how accurate the solution is:
>


 (5.2.4) 



Solving Triangular Systems


Triangular systems can be solved in one step. To solve lower triangular systems, use ForwardSubstitute. To solve upper triangular systems, use BackwardSubstitute. Cannot remember which "Substitute" command to use? As long as the shape of your Matrix is triangular[upper] or triangular[lower], LinearSolve will figure out which of ForwardSubstitute and BackSubstitute to use.
>


>


 (6.1) 
>


>


 (6.2) 
>


>


 (6.3) 
>


>


 (6.4) 
If your system is triangular, but stored with a rectangular shape, specify method=subs, and LinearSolve will use ForwardSubstitute or BackwardSubstitute (as appropriate) to solve the system:
>


>


 (6.5) 
>


>


 (6.6) 
Compute the residual to see how accurate the solution is:
>


 (6.7) 
>



Return to Index for Example Worksheets
