You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by pi...@apache.org on 2003/12/27 16:22:34 UTC
cvs commit: jakarta-commons/math/xdocs/userguide analysis.xml linear.xml
pietsch 2003/12/27 07:22:34
Modified: math/src/java/org/apache/commons/math/analysis
SplineInterpolator.java
math/xdocs navigation.xml
math/xdocs/userguide analysis.xml linear.xml
Log:
Added documentation.
Revision Changes Path
1.11 +6 -8 jakarta-commons/math/src/java/org/apache/commons/math/analysis/SplineInterpolator.java
Index: SplineInterpolator.java
===================================================================
RCS file: /home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/analysis/SplineInterpolator.java,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -r1.10 -r1.11
--- SplineInterpolator.java 19 Nov 2003 03:28:23 -0000 1.10
+++ SplineInterpolator.java 27 Dec 2003 15:22:34 -0000 1.11
@@ -76,11 +76,12 @@
throw new IllegalArgumentException("Dataset arrays must have same length.");
}
+ // TODO: What's this good for? Did I really write this???
if (c == null) {
// Number of intervals. The number of data points is N+1.
int n = xval.length - 1;
// Check whether the xval vector has ascending values.
- // Separation should be checked too (not implemented: which criteria?).
+ // TODO: Separation should be checked too (not implemented: which criteria?).
for (int i = 0; i < n; i++) {
if (xval[i] >= xval[i + 1]) {
throw new IllegalArgumentException("Dataset must specify sorted, ascending x values.");
@@ -88,7 +89,7 @@
}
// Vectors for the equation system. There are n-1 equations for the unknowns s[i] (1<=i<=N-1),
// which are second order derivatives for the spline at xval[i]. At the end points, s[0]=s[N]=0.
- // Vectors are offset by -1, except the lower diagonal vector which is offset by -2. Layout:
+ // Vector indices are offset by -1, except for the lower diagonal vector where the offset is -2. Layout of the equation system:
// d[0]*s[1]+u[0]*s[2] = b[0]
// l[0]*s[1]+d[1]*s[2]+u[1]*s[3] = b[1]
// l[1]*s[2]+d[2]*s[3]+u[2]*s[4] = b[2]
@@ -102,10 +103,6 @@
// Setup right hand side and diagonal.
double dquot = (yval[1] - yval[0]) / (xval[1] - xval[0]);
for (int i = 0; i < n - 1; i++) {
- // TODO avoid recomputing the term
- // (yval[i + 2] - yval[i + 1]) / (xval[i + 2] - xval[i + 1])
- // take it from the previous loop pass. Note: the interesting part of performance
- // loss is the range check in the array access, not the computation itself.
double dquotNext =
(yval[i + 2] - yval[i + 1]) / (xval[i + 2] - xval[i + 1]);
b[i] = 6.0 * (dquotNext - dquot);
@@ -114,7 +111,8 @@
}
// u[] and l[] (for the upper and lower diagonal respectively) are not
// really needed, the computation is folded into the system solving loops.
- // Keep this for documentation purposes:
+ // Keep this for documentation purposes. The computation is folded into
+ // the loops computing the solution.
//double u[] = new double[n - 2]; // upper diagonal
//double l[] = new double[n - 2]; // lower diagonal
// Set up upper and lower diagonal. Keep the offsets in mind.
1.7 +3 -2 jakarta-commons/math/xdocs/navigation.xml
Index: navigation.xml
===================================================================
RCS file: /home/cvs/jakarta-commons/math/xdocs/navigation.xml,v
retrieving revision 1.6
retrieving revision 1.7
diff -u -r1.6 -r1.7
--- navigation.xml 15 Nov 2003 18:38:16 -0000 1.6
+++ navigation.xml 27 Dec 2003 15:22:34 -0000 1.7
@@ -18,8 +18,9 @@
<item name="Statistics" href="/userguide/stat.html"/>
<item name="Data generation" href="/userguide/random.html"/>
<item name="Linear Algebra" href="/userguide/linear.html"/>
+ <item name="Numerical Analysis" href="/userguide/analysis.html"/>
<item name="Special Functions" href="/userguide/special.html"/>
- <item name="Utilities" href="/userguide/utilities.html"/>
+ <item name="Utilities" href="/userguide/utilities.html"/>
</menu>
</body>
</project>
1.6 +164 -24 jakarta-commons/math/xdocs/userguide/analysis.xml
Index: analysis.xml
===================================================================
RCS file: /home/cvs/jakarta-commons/math/xdocs/userguide/analysis.xml,v
retrieving revision 1.5
retrieving revision 1.6
diff -u -r1.5 -r1.6
--- analysis.xml 15 Nov 2003 18:38:16 -0000 1.5
+++ analysis.xml 27 Dec 2003 15:22:34 -0000 1.6
@@ -9,15 +9,57 @@
<body>
<section name="4 Numerical Analysis">
<subsection name="4.1 Overview" href="overview">
- <p>This is yet to be written. Any contributions will be gratefully
- accepted!</p>
+ <p>
+ The numerical analysis package provides basic blocks for common tasks in numerical computation. Currently, only real valued, univariate (depending on one variable) functions are handled.
+ </p>
+ <p>
+ Available blocks:
+ <ul>
+ <li>A framework for solving non-linear equations (root finding)</li>
+ <li>Generating functions by interpolation.</li>
+ </ul>
+ </p>
+ <p>
+ Possible future additions may include numerical differentation, numerical integration and finding minimal or maximal values of a function. Functionality dealing with multivariate functions, complex valued functions or special type functions will probably go into other packages.
+ </p>
</subsection>
<subsection name="4.2 Root-finding" href="rootfinding">
<p>
- <code>org.apache.commons.math.analysis.UnivariateRealSolver</code> provides the means to
- find roots of univariate, real valued, functions. Commons-Math supports various
+ A <code>org.apache.commons.math.analysis.UnivariateRealSolver</code> provides the means to
+ find roots of univariate, real valued, functions. A root is the value where the function
+ value vanishes. Commons-Math supports various
implementations of <code>UnivariateRealSolver</code> to solve functions with differing
- characteristics.
+ characteristics. The current interface allows for computing exactly one root. If the given
+ interval contains more than one root, an indication may be given or not. The current
+ implementations all wont notify the user and return simply the first root found.
+ </p>
+ <p>
+ There are numerous non-obvious traps and pitfalls in root finding. Firstly, the usual
+ disclaimers due to the nature how real world computers calculate values apply. If the
+ computation of the function provides numerical instabilities, for example due to bit
+ cancellation, the root finding algorithms may behave badly and fail to converge or even
+ return bogus values. There wouldn't necessarily be an indication that the computed root
+ is way off the true value. Secondly, The root finding problem itself may be inherently
+ ill conditioned. There is a "domain of indeterminacy", the interval for which the function
+ has near zero absolute values around the true root, may be very large. Even worse, small
+ problems like roundoff error may cause the function value to "numerically oscillate" between
+ negative and positive values. This may again result in roots way off the true value, without
+ indication. There is not much a generic algorithm can do if such ill conditioned problems
+ are met. A way around this is to transform the problem in order to get a better conditioned
+ function.
+ </p>
+ <p>
+ The package provides several implementations off root finding algorithms, each with advantages
+ and drawbacks. The may algorithms differ in
+ <ul>
+ <li>Number of iterations for computing a specific root of a specific function.</li>
+ <li>Number of function evaluations per iteration. An algorithm needing less iterations may still
+ need multiple function evaluations, which may be more costly (for example, involve a numerical
+ integration).</li>
+ <li>Whether the interval must bracket a root or not (function values have different signs at
+ interval endpoints).</li>
+ <li>Behaviour in case of problems (indicate the error, return bogus values...).</li>
+ </ul>
</p>
<p>
In order to use the root-finding features, first a solver object must be created. It is
@@ -26,8 +68,8 @@
of the solver objects supported by Commons-Math. The typical usage of <code>UnivariateRealSolverFactory</code>
to create a solver object would be:</p>
<source>UnivariateRealFunction function = // some user defined function object
- UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
- UnivariateRealSolver solver = factory.newDefaultSolver(function);</source>
+UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
+UnivariateRealSolver solver = factory.newDefaultSolver(function);</source>
<p>
The solvers that can be instantiated via the <code>UnivariateRealSolverFactory</code> are detailed below:
<table>
@@ -42,14 +84,37 @@
methods. For a function <code>f</code>, and two domain values, <code>min</code> and
<code>max</code>, <code>solve</code> computes the value <code>c</code> such that:
<ul>
- <li><code>f(c) = 0.0</code></li>
+ <li><code>f(c) = 0.0</code> (see "function value accuracy")</li>
<li><code>min <= c <= max</code></li>
</ul>
</p>
+ <p>
+ Typical usage:
+ </p>
<source>UnivariateRealFunction function = // some user defined function object
- UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
- UnivariateRealSolver solver = factory.newBisectionSolver(function);
- double c = solver.solve(1.0, 5.0);</source>
+UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
+UnivariateRealSolver solver = factory.newBisectionSolver(function);
+double c = solver.solve(1.0, 5.0);</source>
+ <p>
+ The BrentSolver uses the Brent-Dekker algorithm which is fast and robust. This
+ algorithm is recommended for most users. If there are multiple roots in the interval, or
+ there is a large domain of indeterminacy, the algorithm will converge to a random root in
+ the interval without indication that there are problems. Interestingly, the examined
+ text book implementations all disagree in details of the convergence criteria. Also each implementation
+ had problems for one of the test cases, so the expressions had to be fudged further.
+ Don't expect to get exactly the same root values as for other implementations of this
+ algorithm.
+ </p>
+ <p>
+ The SecantSolver uses a variant of the well known secant algorithm. It may be a bit faster
+ than the Brent solver for a class of well behaved functions.
+ </p>
+ <p>
+ The BisectionSolver is included for completeness and for establishing a fall back in cases
+ of emergency. The algorithm is simple, most likely bug free and guaranteed to converge even
+ in very advert circumstances which might cause other algorithms to malfunction. The drawback
+ is of course that it is also guaranteed to be slow.
+ </p>
<p>
Along with the <code>solve</code> methods, the <code>UnivariateRealSolver</code>
interface provides many properties to control the convergence of a solver. For the most
@@ -58,28 +123,103 @@
is easily done through getter and setter methods on <code>UnivariateRealSolver</code>:
<table>
<tr><th>Property</th><th>Methods</th><th>Purpose</th></tr>
- <tr><td>Absolute accuracy</td><td>
+ <tr>
+ <td>Absolute accuracy</td>
+ <td>
<div>getAbsoluteAccuracy</div>
<div>resetAbsoluteAccuracy</div>
- <div>setAbsoluteAccuracy</div></td><td>This is yet to be written. Any contributions will be greatfully accepted!</td></tr>
- <tr><td>Function value accuracy</td><td>
+ <div>setAbsoluteAccuracy</div>
+ </td>
+ <td>
+ The Absolute Accuracy is maximal difference between the computed root and the true root
+ of the function. This is what most people think of as "accuracy" intuitively. The initial
+ value is choosen as a sane value for most real world problems, for roots in the range from
+ -100 to +100. For accurate computation of roots near
+ zero, in the range form -0.0001 to +0.0001, the value may be decreased. For computing roots
+ much larger in absolute value than 100, the absolute accuracy may never be reached because
+ the given relative accuracy is reached first.
+ </td>
+ </tr>
+ <tr>
+ <td>Relative accuracy</td>
+ <td>
+ <div>getRelativeAccuracy</div>
+ <div>resetRelativeAccuracy</div>
+ <div>setRelativeAccuracy</div>
+ </td>
+ <td>
+ The Relative Accuracy is the maximal difference between the computed root and the true
+ root, divided by the maximum of the absolute values of the numbers. This accuracy
+ measurement is more suited for numerical calculations with computers, due to the way
+ floating point numbers are represented. The default value is choosen so that algorithms
+ will get a result even for roots with large absolute values, even while it may be
+ impossible to reach the given absolute accuracy.
+ </td>
+ </tr>
+ <tr>
+ <td>Function value accuracy</td>
+ <td>
<div>getFunctionValueAccuracy</div>
<div>resetFunctionValueAccuracy</div>
- <div>setFunctionValueAccuracy</div></td><td>This is yet to be written. Any contributions will be greatfully accepted!</td></tr>
- <tr><td>Maximum iteration count</td><td>
+ <div>setFunctionValueAccuracy</div>
+ </td>
+ <td>
+ This value is used by some algorithms in order to prevent numerical instabilities. If the
+ function is evaluated to an absolute value smaller than the Function Value Accuracy, the
+ algorithms assume they hit a root and return the value immediately. The default value is
+ a "very small value". If the goal is to get a near zero function value rather than an
+ accurate root, computation may be speed up by setting this value appropriately.
+ </td>
+ </tr>
+ <tr>
+ <td>Maximum iteration count</td>
+ <td>
<div>getMaximumIterationCount</div>
<div>resetMaximumIterationCount</div>
- <div>setMaximumIterationCount</div></td><td>This is yet to be written. Any contributions will be greatfully accepted!</td></tr>
- <tr><td>Relative accuracy</td><td>
- <div>getRelativeAccuracy</div>
- <div>resetRelativeAccuracy</div>
- <div>setRelativeAccuracy</div></td><td>This is yet to be written. Any contributions will be greatfully accepted!</td></tr>
+ <div>setMaximumIterationCount</div>
+ </td>
+ <td>
+ This is the maximal number of iterations the algorith will try. If this number is exceeded,
+ non-convergence is assumed and an exception is thrown. The default value is 100, which
+ should be plenty giving that a bisection algorithm can't get any more accurate after 52
+ iterations because of the number of mantissa bits in a double precision floating point number.
+ If a number of ill conditioned problems is to be solved, this number can be decreased in order
+ to avoid wasting time.
+ </td>
+ </tr>
</table>
</p>
</subsection>
<subsection name="4.3 Interpolation" href="interpolation">
- <p>This is yet to be written. Any contributions will be gratefully
- accepted!</p>
+ <p>
+ A <code>org.apache.commons.math.analysis.UnivariateRealInterpolator</code> is used to
+ find a univariate, real valued function <code>f</code> which for a given set of pairs
+ (<code>x<sub>i</sub></code>,<code>y<sub>i</sub></code>) yields
+ <code>f(x<sub>i</sub>)=y<sub>i</sub></code> to the best accuracy possible. Currently,
+ only an interpolator for generating natural cubic splines is available. There is
+ no interpolator factory, mainly because the interpolation algorithm is more determined
+ by the kind of the interpolated function rather than the set of points to interpolate.
+ There aren't currently any accuracy controls either.
+ </p>
+ <p>Typical usage:</p>
+ <source>double x[]={ 0.0, 1.0, 2.0 };
+double y[]={ 1.0, -1.0, 2.0);
+UnivariateRealInterpolator interpolator = SplineInterpolator();
+UnivariateRealFunction function = interpolator.interpolate();
+double x=0.5;
+double y=function.evaluate(x);
+System.out println("f("+x+")="+y);</source>
+ <p>
+ A natural spline is a function consisting of a polynominal of third degree for each interval.
+ A function interpolating <code>N</code> value pairs consists of <code>N-1</code> polynominals.
+ The function is continuous, smooth and can be derived two times. The second derivative
+ is continuous but not smooth. The curvature at the first and the last point is zero (that's
+ the "natural" part, coming from the flexible rulers used in construction drawing). The x values
+ passed to the interpolator must be ordered in ascendig order (there is no such restriction for
+ the y values, of course). It is not valid to evaluate the function for values outside the range
+ <code>x<sub>0</sub></code>..<code>x<sub>N</sub></code>. Currently, the original array for the
+ x values is referenced by the generated function, which is probably a bad idea.
+ </p>
</subsection>
</section>
</body>
1.5 +22 -28 jakarta-commons/math/xdocs/userguide/linear.xml
Index: linear.xml
===================================================================
RCS file: /home/cvs/jakarta-commons/math/xdocs/userguide/linear.xml,v
retrieving revision 1.4
retrieving revision 1.5
diff -u -r1.4 -r1.5
--- linear.xml 15 Nov 2003 18:38:16 -0000 1.4
+++ linear.xml 27 Dec 2003 15:22:34 -0000 1.5
@@ -3,34 +3,28 @@
<!-- $Revision$ $Date$ -->
<document url="linear.html">
-<properties>
+ <properties>
<title>The Commons Math User Guide - Linear Algebra</title>
<author email="phil@steitz.com">Phil Steitz</author>
-</properties>
+ </properties>
-<body>
-
-<section name="3 Linear Algebra">
-
-<subsection name="3.1 Overview" href="overview">
- <p>
- This is yet to be written. Any contributions will be gratefully accepted!
- </p>
-</subsection>
-
-<subsection name="3.2 Real matrices" href="real_matrices">
- <p>
- This is yet to be written. Any contributions will be gratefully accepted!
- </p>
-</subsection>
-
-<subsection name="3.3 Solving linear systems" href="solve">
- <p>
- This is yet to be written. Any contributions will be gratefully accepted!
- </p>
-</subsection>
-
-</section>
-
-</body>
+ <body>
+ <section name="3 Linear Algebra">
+ <subsection name="3.1 Overview" href="overview">
+ <p>
+ This is yet to be written. Any contributions will be gratefully accepted!
+ </p>
+ </subsection>
+ <subsection name="3.2 Real matrices" href="real_matrices">
+ <p>
+ This is yet to be written. Any contributions will be gratefully accepted!
+ </p>
+ </subsection>
+ <subsection name="3.3 Solving linear systems" href="solve">
+ <p>
+ This is yet to be written. Any contributions will be gratefully accepted!
+ </p>
+ </subsection>
+ </section>
+ </body>
</document>
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org