Download Numerical Linear Algebra: Solving Linear Systems and Matrix Decomposition and more Study notes Astronomy in PDF only on Docsity!
Numerical Linear AlgebraNumerical Linear Algebra
Probably the simplest kind of problem.
Occurs in many contexts, often as part of larger
problem.
Symbolic manipulation packages can do linear
algebra "analytically" (e.g. Mathematica, Maple).
Numerical methods needed when:
(^) Number of equations very large (^) Coefficients all numerical
Linear SystemsLinear Systems
⋮ ⋮
Write linear system as:
(^) This system has n unknowns and m equations. (^) If n = m , system is closed. (^) If any equation is a linear combination of any others, equations are degenerate and system is singular.* *see Singular Value Decomposition (SVD), NRiC 2.6.
Matrix FormMatrix Form
Write system in matrix form:
where: =
⋯ ⋯ ⋮ ⋮ ⋮ ⋯ Columns Rows
Matrix Data RepresentationMatrix Data Representation
Recall, C stores data in row-major form:
a 11 , a 12 , ..., a 1 n ; a 21 , a 22 , ..., a 2 n ; ...; a m 1 , a m 2 , ..., a mn
If using "pointer to array of pointers to rows"
scheme in C, can reference entire rows by first
index, e.g. 3
rd
row = a[2].
(^) Recall in C array indices start at zero!!
FORTRAN stores data in column-major form:
a 11 , a 21 , ..., a m 1 ; a 12 , a 22 , ..., a m 2 ; ...; a 1 n , a 2 n , ..., a mn
Tasks of Linear AlgebraTasks of Linear Algebra
We will consider the following tasks:
- Solve A x = b , given A and b.
- Solve A x i = b i for multiple b i 's.
- Calculate A
- , where A - A = I , the identity matrix.
- Calculate determinant of A , det( A ).
Large packages of routines available for these
tasks, e.g. LINPACK, LAPACK (public domain);
IMSL, NAG libraries (commercial).
We will look at methods assuming n = m.
The Augmented MatrixThe Augmented Matrix
The equation A x = b can be generalized to a form
better suited to efficient manipulation:
The system can be solved by performing
operations on the augmented matrix.
The x
i
's are placeholders that can be omitted until
the end of the computation.
∣ = ⋯ ∣ ⋯ ∣ ⋮ ⋮ ⋮ ∣ ⋮ ⋯ ∣
The Generalized Matrix EquationThe Generalized Matrix Equation
Consider the generalized linear matrix equation:
Its solution simultaneously solves the linear sets:
A x 1 = b 1 , A x 2 = b 2 , A x 3 = b 3 , and AY = I ,
where the x
i
's and b
i
's are column vectors.
∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ = ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣
Gauss-Jordan EliminationGauss-Jordan Elimination
GJE uses one or more elementary row operations
to reduce matrix A to the identity matrix.
The RHS of the generalized equation becomes the
solution set and Y becomes A
Disadvantages:
- Requires all b i 's to be stored and manipulated at same time ⇒ memory hog.
- Don't always need A
Other methods more efficient, but good backup.
GJE Procedure, Cont'dGJE Procedure, Cont'd
Problem occurs if leading diagonal element ever
becomes zero.
Also, procedure is numerically unstable!
Solution: use "pivoting" - rearrange remaining
rows (partial pivoting) or rows & columns (full
pivoting - requires permutation!) so largest
coefficient is in diagonal position.
Best to "normalize" equations (implicit pivoting).
Gaussian Elimination withGaussian Elimination with
Backsubstitution Backsubstitution
If, during GJE, only subtract rows below pivot,
will be left with a triangular matrix:
(^) Solution for x 3 is then trivial: x 3 = b 3 '/ a 33
(^) Substitute into 2 nd row to get x 2
(^) Substitute x 3 & x 2 into 1 st row to get x 1
Faster than GJE, but still memory hog.
^ = "Gaussian Elimination"
LULU Decomposition, Cont'dDecomposition, Cont'd
Why is this better?
(^) Solving triangular matrices is easy: just use forward substitution for (1), backsubstitution for (2).
Problem is, how to decompose A into L and U?
(^) Expand matrix multiplication LU to get n 2 equations for n 2
- n unknowns (elements of L and U plus n extras because diagonal elements counted twice). (^) Get an extra n equations by choosing L ii = 1 ( i = 1, n ). (^) Then use Crout's algorithm for finding solution to these n 2
- n equations "trivially" (NRiC 2.3).
LULU Decomposition in NRiCDecomposition in NRiC
The routines ludcmp() and lubksb() perform LU
decomposition and backsubstitution respectively.
Can easily compute A
(solve for the identity
matrix column by column) and det( A ) (find the
product of the diagonal elements of the LU
decomposed matrix) - see NRiC 2.3.
WARNING: for large matrices, computing det( A )
can overflow or underflow the computer's
floating-point dynamic range.
Tridiagonal MatricesTridiagonal Matrices
Many systems can be written as (or reduced to):
a i x i -
- b i x i
- c i x i + = d i i = 1, n
i.e. a tridiagonal matrix:
Here a
1
and c
n
are associated with "boundary
conditions" (i.e. x
0
and x
n +
[
⋱ ⋱ ⋱ − − −
][
⋮ (^) − (^)
]
=
[
(^) (^) (^) ⋮ (^) − (^)
]
Sparse MatricesSparse Matrices
LU decomposition and backsubstitution is very
efficient for tri-di systems: O ( n ) operations as
opposed to O ( n
3
) in general case.
Operations on sparse systems can be optimized.
e.g. Tridiagonal Band diagonal with bandwidth M Block diagonal Banded
See NRiC 2.7 for various systems & techniques.
