6. Class Numerical Linear Algebra, Part 1
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, etc.). Numerical methods needed when: Number of equations very large. One or more coecients numerical.
Linear Systems
Write linear system as: a11 x1 + a12 x2 + + a1n xn = b1 a21 x1 + a22 x2 + + a2n xn = b2 . .=. . . . am1 x1 + am2 x2 + + amn xn = bm This system has n unknowns and m equations. If n = m, system is closed. If m n and any equation is a linear combination of any others, equations are degenerate and system is singular.
Numerical Constraints
Numerical methods have their own problems when: 1. Equations are degenerate within round-o error. 2. Accumulated round-o errors swamp solution (magnitudes of as and xs vary wildly). For n, m < 50, single precision usually OK (but why bother?). For n, m < 200, double precision usually OK. For 200 < n, m < few thousand, solutions possible only for sparse systems (lots of as zero).
1
Matrix Form
Write system in matrix form: Ax = b, where:
A=
a11 a21 . . .
a12 a22 . . .
a1n a2n . . .
am1 am2 amn
.
Matrix Data Representation
Recall, C stores data in row-major form: a11 , a12 , . . . , a1n ; a21 , a22 , . . . , a2n ; . . . ; am1 , am2 , . . . , amn . If using pointer to array of pointers to rows scheme in C, can reference entire rows by rst index, e.g., 3rd row = a[2]. (!) Recall in C array indices start at zero! FORTRAN stores data in column-major form: a11 , a21 , . . . , am1 ; a12 , a22 , . . . , am2 ; . . . ; a1n , a2n , . . . , amn .
Note on Numerical Recipes in C
The canned routines in NRiC make use of special functions dened in nrutil.c (header nrutil.h). In particular, arrays and matrices are allocated dynamically with indices starting at 1, not 0. If you want to interface with the NRiC routines, but prefer the normal C array index convention, pass arrays by subtracting 1 from the pointer address (i.e., pass p - 1 instead of p) and pass matrices by using the functions convert matrix() and free convert matrix() in nrutil.c (see NRiC 1.2 for more information).
Tasks of Linear Algebra
We will consider the following tasks: 1. Solve Ax = b, given A and b. 2. Solve Axi = bi for multiple bi s. 3. Calculate A1 , where A1 A = 1, the identity matrix. 4. Calculate the determinant of A, det(A).
2
Large packages of routines available for these tasks, e.g., LINPACK, LAPACK, GSL (public domain), IMSL, NAG libraries (commercial). We will look at methods assuming n = m.
The Augmented Matrix
The equation Ax = b can be generalized to a form better suited to ecient manipulation: a12 a11 a1n b1 a21 a22 a2n b2 . (A|b) = . . . . . . . . . . . . an1 an2 ann bn The system can be solved by performing operations on the augmented matrix. The xi s are placeholders that can be omitted until the end of the computation. Elementary row operations The following row operations can be performed on an augmented matrix without changing the solution of the underlying system of equations: 1. Interchange two rows. 2. Multiply a row by a nonzero real number. 3. Add a multiple of one row to another row. The idea is to apply these operations in sequence until the system of equations is trivially solved. The generalized matrix equation Consider the generalized linear matrix equation:
a11 a21 a31 a41 a12 a22 a32 a42 a13 a23 a33 a43 a14 a24 a34 a44 x11 x21 x31 x41 x12 x22 x32 x42 x13 x23 x33 x43 y11 y21 y31 y41 y12 y22 y32 y42 y13 y23 y33 y43 y14 y24 y34 y44 = b11 b21 b31 b41 b12 b22 b32 b42 b13 b23 b33 b43 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 .
coefficients
solutions and inverse
RHS and identity
Its solution simultaneously solves the linear sets: Ax1 = b1 , Ax2 = b2 , Ax3 = b3 , and AY = 1, where the xi s and bi s are column vectors.
3
Gauss-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 A1. Disadvantages: 1. Requires all bi s to be stored and manipulated at same time memory hog. 2. Dont always need A1 . Other methods more ecient, but good backup. Procedure Start with simple augmented matrix as example:
a11 a12 a13 b1 a21 a22 a23 b2 a31 a32 a33 b3 Divide rst ro...