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$

Author:

Chris Moore, Roland Keller

Nested Class Summary

Nested Classes

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

Field Summary

Fields

Modifier and Type	Field and Description
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

Constructors

Constructor and Description
MatrixOperations()

Method Summary

Modifier and Type	Method and Description
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

equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

ALF

public static final double ALF

Ensures sufficient decrease in function value

See Also:

Constant Field Values

EPS

public static final double EPS

Approximate square root of the JVM precision

See Also:

Constant Field Values

STPMX

public static final double STPMX

Scaled maximum step length allowed in line searches

See Also:

Constant Field Values

TINY

public static final double TINY

Extremely small value. Nigh zero.

See Also:

Constant Field Values

TOLF

public static final double TOLF

Convergence criterion on function values

See Also:

Constant Field Values

TOLMIN

public static final double TOLMIN

Criterion deciding whether spurious convergence to a minimum of min

See Also:

Constant Field Values

TOLX

public static final double TOLX

Convergence criterion on delta X

See Also:

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