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

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:

• Concepts from "linear algebra" that determine when a system of equations will have a unique solution, multiple solutions, and no solutions.
• The method of LU decomposition for solving systems of linear equations.

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
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

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