Solving Linear Equations [A].[x] = [b]
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.
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."
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.
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
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