/*
* @(#)Equations.java
*
* $Date: 2012-07-03 01:10:05 -0500 (Tue, 03 Jul 2012) $
*
* Copyright (c) 2011 by Jeremy Wood.
* All rights reserved.
*
* The copyright of this software is owned by Jeremy Wood.
* You may not use, copy or modify this software, except in
* accordance with the license agreement you entered into with
* Jeremy Wood. For details see accompanying license terms.
*
* This software is probably, but not necessarily, discussed here:
* http://javagraphics.java.net/
*
* That site should also contain the most recent official version
* of this software. (See the SVN repository for more details.)
*/
package com.bric.math;
import java.util.Arrays;
import java.util.Comparator;
/** This contains a static method to solve a system of simple linear equations.
*
*/
public class Equations {
/** Generating the text in certain exceptions is not a trivial task if you
* expect those exceptions to be called millions of times.
*
*/
public static boolean VERBOSE_EXCEPTIONS = true;
private static Comparator<double[]> coefficientComparator = new Comparator<double[]>() {
public int compare(double[] d1, double[] d2) {
int v1 = 0;
int v2 = 0;
int a;
for( a = 0; a<d1.length; a++) {
if(d1[a]==1) {
v1 = a;
a = d1.length;
}
}
for( a = 0; a<d2.length; a++) {
if(d2[a]==1) {
v2 = a;
a = d2.length;
}
}
return v1-v2;
}
};
public static String toString(double[][] d) {
String s = "";
for(int a = 0; a<d.length; a++) {
s = s+toString(d[a])+"\n";
}
return s.trim();
}
public static String toString(double[] d) {
String s = "[";
for(int a = 0; a<d.length; a++) {
if(a==0) {
s = s+" "+d[a];
} else {
s = s+", "+d[a];
}
}
return s+" ]";
}
/** Given a matrix of variable coefficients for a linear system of equations,
* this will solve for each variable.
* <P> For example, if you have the three equations:
* <BR> 2x + y + z = 1
* <BR> 6x + 2y + z = -1
* <BR> -2x + 2y + z = 7
* <P> You can call:
* <code><BR>double[][] coeffs = { { 2, 1, 1, 1},</code>
* <code><BR> { 6, 2, 1, -1},</code>
* <code> <BR> { -2, 2, 1, 7} };</code>
* <code><BR>Equations.solve(coeffs,true);</code>
* <P> Then this method will simplify the matrix to:
* <code><BR>{ { 1, 0, 0, -1}, </code>
* <code><BR> { 0, 1, 0, 2}, </code>
* <code><BR> { 0, 0, 1, 1} }</code>
* <P> This indicates that:
* <BR> 1*x+0*y+0*z = -1
* <BR> 0*x+1*y+0*z = 2
* <BR> 0*x+0*y+1*z = 1
* <P> This method uses Gaussian Elimination and back substitution.
* <P> Note that due to the limits of double-precision, some safeguards were
* built-in that may round values less than .00000001 to zero. If you
* need better precision, you may have to copy and paste the algorithm
* in this code and use a more precise number format.
* @param coefficients the matrix of coefficients. It must be N x (N+1) units in size.
* @throws IllegalArgumentException if it is impossible to solve for each variable. In this case an
* exception is thrown, and no guarantee is made regarding the final state of the
* <code>coefficients</code> matrix.
* <P> Also an exception may be thrown if the matrix is not correctly sized.
*/
public static void solve(double[][] coefficients) {
solve(coefficients,true);
}
/** Given a matrix of variable coefficients for a linear system of equations,
* this will solve for each variable.
* <P> For example, if you have the three equations:
* <BR> 2x + y + z = 1
* <BR> 6x + 2y + z = -1
* <BR> -2x + 2y + z = 7
* <P> You can call:
* <code><BR>double[][] coeffs = { { 2, 1, 1, 1},</code>
* <code><BR> { 6, 2, 1, -1},</code>
* <code> <BR> { -2, 2, 1, 7} };</code>
* <code><BR>Equations.solve(coeffs,true);</code>
* <P> Then this method will simplify the matrix to:
* <code><BR>{ { 1, 0, 0, -1}, </code>
* <code><BR> { 0, 1, 0, 2}, </code>
* <code><BR> { 0, 0, 1, 1} }</code>
* <P> This indicates that:
* <BR> 1*x+0*y+0*z = -1
* <BR> 0*x+1*y+0*z = 2
* <BR> 0*x+0*y+1*z = 1
* <P> This method uses Gaussian Elimination and back substitution.
* <P> Note that due to the limits of double-precision, some safeguards were
* built-in that may round values less than .00000001 to zero. If you
* need better precision, you may have to copy and paste the algorithm
* in this code and use a more precise number format.
* @param coefficients the matrix of coefficients. It must be N x (N+1) units in size.
* @param sort If this is <code>true</code>, then you're guaranteed for your final solution
* to be sorted such that the 1's form a diagonal line from [0,0] to [N,N]. If this
* is <code>false</code>, then the solution <i>may</i> not form a diagonal line of 1's.
* <P> (If there are no zeroes in the matrix, the diagonal line will be there by default,
* but if zeroes are present then some rows may have to be skipped and the normal order
* is disrupted.)
* <P> The sort is performed by calling <code>Arrays.sort()</code>.
* @throws IllegalArgumentException if it is impossible to solve for each variable. In this case an
* exception is thrown, and no guarantee is made regarding the final state of the
* <code>coefficients</code> matrix.
* <P> Also an exception may be thrown if the matrix is not correctly sized.
*/
public static void solve(double[][] coefficients,boolean sort) {
if(coefficients==null) throw new NullPointerException("The coefficients matrix is null.");
int size = coefficients.length;
//this will keep track of which rows have been processed
//in an ideal world, we would like to go sequentially down the list,
//but if a coefficient somewhere is already zero, then we may have to skip around.
boolean[] b = new boolean[coefficients.length];
//this keeps track of which row we solved, from first to last
//we use this to back substitute
int[] order = new int[b.length];
int ctr = 0;
int row = 0;
int a, i;
double t;
int errorCounter = 0;
//println(coefficients);
while(ctr<b.length) {
if(coefficients[row].length!=size+1)
throw new IllegalArgumentException("The matrix must be N x (N+1) units long. The matrix provided is "+size+" x "+coefficients[row].length+" units.");
if(b[row]==false && Math.abs(coefficients[row][ctr])>.0000000001) {
errorCounter = 0;
//we haven't considered this row yet
//first we make the value at [row][ctr] a 1.0, by multiplying the entire row by a constant:
t = 1/coefficients[row][ctr];
for(a = 0; a<coefficients[row].length; a++) {
coefficients[row][a] *= t;
}
coefficients[row][ctr] = 1.0; //to make sure, despite possible arithmetic rounding, that we get 1.0
b[row] = true;
for(a = 0; a<coefficients.length; a++) {
if(b[a]==false) {
//this is a row we need to adjust
t = coefficients[a][ctr];
for(i = 0; i<coefficients[a].length; i++) {
coefficients[a][i] -= coefficients[row][i]*t;
}
}
}
//we've now considered this row
order[ctr++] = row;
}
errorCounter++;
row++;
row = row%(coefficients.length);
if(errorCounter>coefficients.length) {
if(VERBOSE_EXCEPTIONS) {
throw new IllegalArgumentException("The coefficient matrix cannot be solved. Either it has infinitely many solutions, or zero solutions:\n"+toString(coefficients));
} else {
throw new IllegalArgumentException("The coefficient matrix cannot be solved. Either it has infinitely many solutions, or zero solutions.");
}
}
//System.out.println("\trow = "+row);
//println(coefficients);
}
//System.out.println("\t\tmoving on...");
//coefficients now contains an array that looks like this:
// 1 a b c ... x
// 0 1 d e ... x
// 0 0 1 f ... x
// 0 0 0 1 ... x
// .......... 1 x
//so now we have to clear out all the other values by backtracking.
//we stored the order we calculated everything in in the "order" array
//the last element in the order array points to the row
//in our matrix that reads "0 ... 0 1 z". In other words,
// it is the only solved row.
//So we can skip it, and move to the row that reads: "0 ... 1 a b".
ctr = 0;
for(a = order.length-2; a>=0; a--) {
row = order[a];
for(i = coefficients[row].length-2; i>a; i--) {
t = coefficients[row][i]*coefficients[order[i]][coefficients[row].length-1];
coefficients[row][coefficients[row].length-1] -= t;
coefficients[row][i] = 0;
}
//System.out.println("\trow = "+row);
//println(coefficients);
}
if(sort) {
Arrays.sort(coefficients,coefficientComparator);
}
}
}