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.)
restart
withLinearAlgebra:
UseHardwareFloats≔false:
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:
M:=0|3|1,7|−13|−2,1|2|4
V:=5.84,−8.46,13.11
M≔0317−13−2124
V≔5.84−8.4613.11
To solve the linear system, you can use the LinearSolve command. Note that you can use separate Matrix and right-hand side Vector arguments, or simply use an augmented Matrix. To ensure that the method used is LU decomposition (the step-by-step method demonstrated in the next section), specify method='LU'.
x:=LinearSolveM,V,method='LU';x:=LinearSolveM|V,method='LU'
x≔1.6499999961.1899999982.270000002
Compute the residual to see how accurate the solution is:
M.x−V
−4.×10−9−4.×10−90.
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:
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).
P,L,U:=LUDecompositionM
P,L,U≔010100001,10001017971,7−13−2031003
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.
V2:=TransposeP.V
V2≔−8.465.8413.11
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:
V3:=ForwardSubstituteL,V2
V3≔−8.4600000005.8400000006.810000000
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:
x:=BackwardSubstituteU,V3
x≔1.6500000001.1900000002.270000000
0.0.0.
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:
To solve the linear system, you can use the LinearSolve command. To ensure that the method used is QR decomposition (the step-by-step method demonstrated in the next section), specify method='QR'.
x:=LinearSolveM,V,method='QR'
x≔1.6500000041.1900000012.269999999
2.×10−92.2×10−80.
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.
First, do a QR decomposition of M. The output components are Q (an orthogonal Matrix) and R (an upper triangular Matrix).
Q,R:=QRDecompositionM
Q,R≔05262131−91311317210−92621310−1311312106326213107131131,52−89210−20326210322621310021131131
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.
V2:=TransposeQ.V
V2≔−4.61100000020.91150381682620.3638931298131
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:
x:=BackwardSubstituteR,V2
x≔1.6500000021.1900000012.269999999
2.×10−92.×10−90.
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:
M_temp:=Matrix4,i,j→i+ij−7,shape=triangularlower
M:=M_temp.TransposeM_temp
IsMatrixShapeM,symmetric
IsDefiniteM
true
You will also need a right-hand side vector:
V:=6,1,3,−2
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.
LinearSolveM,V
x≔LinearSolveM,V,method='Cholesky'
Solving Step-By-Step
To solve the system one step at a time, start with a Cholesky decomposition:
L:=LUDecompositionM,method='Cholesky'
Use ForwardSubstitute to solve L.V2 = V:
V2:=ForwardSubstituteL,V
Use BackwardSubstitute to solve Transpose(L).x=V2:
x:=BackwardSubstituteTransposeL,V2
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.
M≔Matrix3,i,j→i+2j−8,shape=triangularlower
V≔7,8,1
x:=ForwardSubstituteM,V
x≔LinearSolveM,V
M:=Matrix3,i,j→ij+6,shape=triangularupper
V:=6,−3,14
x:=BackwardSubstituteM,V
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:
M≔ −3|1|6,0|2|7,0|0|−4
V≔1,3,5
x:=LinearSolveM,V,method=subs
Return to Index for Example Worksheets
Download Help Document
