[ 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