org.simulator.math Class MatrixOperations

java.lang.Object org.simulator.math.MatrixOperations

public class MatrixOperations
extends Object

Class used to perform matrix operations, focusing on finding vector solutions to the vector equation F(x) = 0.

TODO: Add calculation of eigenvectors. This isn't hard to implement as we've already calculated the eigenvalues and this means all we need to do is solve for x in the matrix equation (A-lambda*I)*x = 0. However, if we have complex eigenvalues then we have a problem. Java doesn't have complex number support natively.

Class method summary:

Numerical Methods. Only call these directly if you know what you're doing.
• mnewt - Uses a straightforward Newton-Rhapson method to compute the equilibrium solution iteratively starting at a reasonably close guess.
• ludcmp - Performs LU decomposition of a matrix (A = L*U) to prepare it for solving the linear equation A*x=B for x.
• lubksb - Performs back substitution on an LU decomposed matrix to solve the linear equation A*x=B for x.
• fdjac - Numerically approximates the jacobian for a system of equations using the method of forward differences.
• fmin - Not currently used. Computes 1/2 the dot product of a vector function with itself at a specified x vector.
• elmhes - Reduces the given matrix to upper Hessenberg form.

Utility Functions. Call these methods!

Notes: The majority of the code in this class comes from Numerical Recipes in C, 2nd Edition. The code was translated from C to Java by Eric Harley. Mostly this amounted to figuring out how to deal with function pointers to user defined functions and re-indexing things starting at 0 instead of 1.

References: _Numerical Recipies in C_ (2nd Edition) by Press, Teutolsky, Vetterling, and Flannery Cambridge University Press, 1992

Since:
0.9
Version:
\$Rev: 168 \$
Author:
Chris Moore, Roland Keller

Nested Class Summary
static class MatrixOperations.MatrixException
This exception is thrown when errors in the computation of matrix-related solutions, their eigenvalues or eigenvectors.

Field Summary
static double ALF
Ensures sufficient decrease in function value
static double EPS
Approximate square root of the JVM precision
static double STPMX
Scaled maximum step length allowed in line searches
static double TINY
Extremely small value.
static double TOLF
Convergence criterion on function values
static double TOLMIN
Criterion deciding whether spurious convergence to a minimum of min
static double TOLX
Convergence criterion on delta X

Constructor Summary
MatrixOperations()

Method Summary
static void balance(double[][] a)
Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix with identical eigenvalues.
static void elmhes(double[][] a)
Reduces a[1..n][1..n] to upper Hessenberg form
static int hqr(double[][] a, double[] wr, double[] wi)
Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n].
static void lubksb(double[][] a, int[] indx, double[] b)
Solves the set of n linear equations AX = B.
static double ludcmp(double[][] a, int[] indx)
Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output,l arrand as in equation (2.3.14) above; indx[1..n] is an output vector that records the row permutation effected by the partial pivoting; d (return value) is output +/- 1 depending on whether the number of row interchanges was even or odd respectively.
static double sign(double a, double b)
Returns the value of a or |a| with the same sign as b

Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

ALF

public static final double ALF
Ensures sufficient decrease in function value

Constant Field Values

EPS

public static final double EPS
Approximate square root of the JVM precision

Constant Field Values

STPMX

public static final double STPMX
Scaled maximum step length allowed in line searches

Constant Field Values

TINY

public static final double TINY
Extremely small value. Nigh zero.

Constant Field Values

TOLF

public static final double TOLF
Convergence criterion on function values

Constant Field Values

TOLMIN

public static final double TOLMIN
Criterion deciding whether spurious convergence to a minimum of min

Constant Field Values

TOLX

public static final double TOLX
Convergence criterion on delta X

Constant Field Values
Constructor Detail

MatrixOperations

public MatrixOperations()
Method Detail

ludcmp

public static double ludcmp(double[][] a,
int[] indx)
throws MatrixOperations.MatrixException
Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output,l arrand as in equation (2.3.14) above; indx[1..n] is an output vector that records the row permutation effected by the partial pivoting; d (return value) is output +/- 1 depending on whether the number of row interchanges was even or odd respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix.

Parameters:
a - the matrix to be decomposed
indx - the array to put the return index into
Returns:
+1 or -1 signifying whether the number of row interchanges is odd or even
Throws:
MatrixOperations.MatrixException

lubksb

public static void lubksb(double[][] a,
int[] indx,
double[] b)
Solves the set of n linear equations AX = B. Here a[1..n][1..n] is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx[1..n] is input as the permutation vector returned by ludcmp. b[1..n] is input as the right hand side vector B, and returns with the solution vector X. a, and indx are not modified by this routine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that b will begin with many zero elements, so it is efficient for use in matrix inversion.

Parameters:
a - the matrix to be solved as described
indx - the array returned by ludcmp
b - the vector to be solbed as described

balance

public static void balance(double[][] a)
Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix with identical eigenvalues. A symmetric matrix is already balanced and is unaffected by this procedure. The parameter RADIX should be the machine's floating-point radix.

Parameters:
a - the matrix to be balanced

elmhes

public static void elmhes(double[][] a)
Reduces a[1..n][1..n] to upper Hessenberg form

Parameters:
a - the matrix to be reduced

sign

public static double sign(double a,
double b)
Returns the value of a or |a| with the same sign as b

Parameters:
a - the input as specified above
b - the input as specified above
Returns:
the value of a or |a| with the same sign as b

hqr

public static int hqr(double[][] a,
double[] wr,
double[] wi)
throws MatrixOperations.MatrixException
Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n]. On input a can be exactly as output from elmhes (� 11.5); on output it is destroyed. The real and imaginary parts of the eigenvalues are returned in wr[1..n] and wi[1..n], respectively.

Parameters:
a - the input matrix
wr - the array specified in the function description
wi - the array specified in the function description
Returns:
0 iff success
Throws:
MatrixOperations.MatrixException