Solving Linear Equations [A].[x] = [b]

[ Solving [A].[x] = [b] ] [ [L][U] Decomposition ]
[ Solving [A].[x] = [b] for multiple [b]'s ] [ Inverse of Matrix [A] ]

In this page we demonstrate use of the matrix functions for the solution of linear systems of equations. We will assume readers are familiar with:

A review of the theoretical concepts needed to solve linear matrix equations may be found in Technical Report T.R. 95-74.


SOLVING LINEAR EQUATIONS [A].[x] = [b]

Let "A" be a (nxn) matrix, and "b" be a (nx1) matrix. The most straightforward way of solving the matrix equations

    [ A ].[ x ] = [ b ]
is with
    FUNCTION           PURPOSE
    ==============================================================

    Solve( A, b )      Solve matrix equations [A].[x] = [b].

    ==============================================================

Example 1 : The script of code

    /* [a] : Setup matrices [A] and [b] */ 

    A = [ 1 , 2; 3 , 4 ];
    B = [ 5 ; 11 ];

    PrintMatrix( A, B );

    /* [b] : Solve matrices and print solution vector */ 

    X = Solve (A, B);
    PrintMatrix( X );
generates the output
    MATRIX : "A"

    row/col                  1            2   
            units                             
       1            1.00000e+00  2.00000e+00
       2            3.00000e+00  4.00000e+00

    MATRIX : "B"

    row/col                  1   
            units                
       1            5.00000e+00
       2            1.10000e+01

    MATRIX : "X"

    row/col                  1   
           units                
       1            1.00000e+00
       2            2.00000e+00

Example 2 : Now let's solve a system of equations having units. The script of code (this example only works for ALADDIN 2.0, scheduled for release in July 1997):

    /* [a] : Structure Calculation */

    E = 200 MPa;
    I = 0.4 m^4;
    A = 0.1 m^2;
    L = 5.0 m;

    /* [b] : Setup and print 3x3 stiffness matrix */

    stiff = Matrix([3,3]);
    stiff = ColumnUnits( stiff, [N/m, N/m, N/rad]);
    stiff = RowUnits( stiff,     [m], [3] );

    stiff[1][1] =  E*A/L;
    stiff[2][2] = 12*E*I/L^3;
    stiff[2][3] = -6*E*I/L/L;
    stiff[3][2] = -6*E*I/L/L;
    stiff[3][3] =  4*E*I/L;

    PrintMatrix(stiff);

    /* [c] : Setup and print 3x1 force vector  */

    force = [ 0.0 N ; -5 N; 0.0 N*m ];
    PrintMatrix(force);

    /* [d] : Compute and print displacement matrix */

    displ = Solve(stiff, force ); 
    PrintMatrix(displ);
    PrintMatrixCast(displ, [ cm; deg ] );

generates the output:

    MATRIX : "stiff"

    row/col                  1            2            3   
            units          N/m          N/m        N/rad   
       1            4.00000e+06  0.00000e+00  0.00000e+00
       2            0.00000e+00  7.68000e+06 -1.92000e+07
       3        m   0.00000e+00 -1.92000e+07  6.40000e+07

    MATRIX : "force"

    row/col                  1   
            units                
       1        N   0.00000e+00
       2        N  -5.00000e+00
       3      N.m   0.00000e+00
    
    MATRIX : "displ"

    row/col                  1   
            units                
       1        m   0.00000e+00
       2        m  -2.60417e-06
       3      rad  -7.81250e-07

    MATRIX : "displ"

    row/col                  1   
        units                
       1       cm   0.00000e+00
       2       cm  -2.60417e-04
       3      deg  -4.47623e-05

In the final block of matrix output, radians are scaled to degrees, and translational displacements are expressed in "cm."


[L][U] DECOMPOSITION

In the method of LU decomposition, matrix "A" is decomposed into a product of lower and upper triangular matrices.

    Step 1 : Decompose [ A ]  ====>  [ L ].[ U ]

The solution vector "x" is computed via forward substitution and then backward susbstitution

    Step 2 : Solve [ L ].[ z ] = [ b ] : Forward  Substitution
    Step 3 : Solve [ U ].[ x ] = [ x ] : Backward Substitution

Matrix functions are provided for the LU decomposition and forward/backward substitution.

    FUNCTION             PURPOSE
    ==============================================================

    Decompose( A )       Decompose "A" into a product of lower and
                         upper triangular matrices.

    Substitution ( LU )  Compute forward/backward substitution on
                         decomposed LU matrix product.

    ==============================================================

Example 3 : The script of ALADDIN commands

    /* [a] : Setup matrices [A] and [b] */ 

    A = [ 1 , 2; 3 , 4 ];
    B = [ 5 ; 11 ];

    PrintMatrix( A, B );

    /* [b] : Solve matrices and print solution vector */ 

    LU = Decompose (A);
    X  = Substitution ( LU, B);
    PrintMatrix( X );

produces output identical to Example 1.


SOLVING [A].[x] = [b] for multiple [b]'s

In the method of LU decomposition, solutions to multiple right-hand vectors can be computed once the matrix [A] has been decomposed.

Example 4 : The script

    /* [a] : Setup matrices [A], [B1], and [B2] */

    A = [ 1 , 2; 3 , 4 ];
    B1 = [ 5 ; 11 ];
    B2 = [ 1 ;  2 ];

    PrintMatrix( A );

    /* [b] : Decompose "A" into the "LU" matrix product */

    LU = Decompose (A);

    /* [c] : Solve matrix equations via forward/backward substitution */

    X1  = Substitution ( LU, B1);
    PrintMatrix( B1, X1 );

    X2  = Substitution ( LU, B2);
    PrintMatrix( B2, X2 );
solves two sets of equations A.x = B, with A as defined in Examples 1 and 2. The abbreviated output is:
    SOLVING [A].[X1] = [B1]              SOLUTION VECTOR
    ================================================================

    MATRIX : "B1"                        MATRIX : "X1"

    row/col                  1           row/col                  1   
            units                                units                
       1            5.00000e+00             1            1.00000e+00
       2            1.10000e+01             2            2.00000e+00

    SOLVING [A].[X2] = [B2]              SOLUTION VECTOR
    ================================================================

    MATRIX : "B2"                        MATRIX : "X2"

    row/col                  1           row/col                  1   
            units                                units                
       1            1.00000e+00             1            0.00000e+00
       2            2.00000e+00             2            5.00000e-01


INVERSE OF MATRIX [A]

The inverse of a square matrix "A" is computed with

    FUNCTION           PURPOSE
    ==============================================================

    Inverse( A )       Compute inverse of square matrix [A]

    ==============================================================

Example 5 : The script of code

    /* [a] : Setup matrix [A] and compute [A]^-1 */ 

    A = [ 1,  2,  3;
          7, 10,  3;
          3,  7,  8 ];

    Ainv = Inverse (A);
    PrintMatrix( A , Ainv);

    /* [b] : Check results by computing [A].[Ainv] */ 

    PrintMatrix( A*Ainv );

defines a (3x3) matrix [A], and computes its inverse. We check that [A].[Ainv] equals [I], a (3x3) identity matrix. The output is:

    MATRIX : "A"

    row/col                  1            2            3   
            units                                          
       1            1.00000e+00  2.00000e+00  3.00000e+00
       2            7.00000e+00  1.00000e+01  3.00000e+00
       3            3.00000e+00  7.00000e+00  8.00000e+00

    MATRIX : "Ainv"

    row/col                  1            2            3   
            units                                          
       1            2.68182e+00  2.27273e-01 -1.09091e+00
       2           -2.13636e+00 -4.54545e-02  8.18182e-01
       3            8.63636e-01 -4.54545e-02 -1.81818e-01

    MATRIX : "UNTITLED"

    row/col                  1            2            3   
            units                                          
       1            1.00000e+00  0.00000e+00  0.00000e+00
       2            8.88178e-16  1.00000e+00  4.44089e-16
       3           -8.88178e-16  0.00000e+00  1.00000e+00

Example 6 : Now let's add units to matrix X in Example 5.

If [A] is defined as

    /* [a] : Setup matrix [A] and compute [A]^-1 */ 

    A = ColumnUnits ([ 1,  2,  3;
                       7, 10,  3;
                       3,  7,  8 ], [ N/m, N/m, m ]);

then the inverse of A is

    MATRIX : "Ainv"

    row/col                  1            2            3   
            units                                          
       1      m/N   2.68182e+00  2.27273e-01 -1.09091e+00
       2      m/N  -2.13636e+00 -4.54545e-02  8.18182e-01
       3      1/m   8.63636e-01 -4.54545e-02 -1.81818e-01


Developed in April 1996 by Mark Austin
Last Modified February 14, 1997
Copyright © 1996-1997, Mark Austin, Department of Civil Engineering, University of Maryland