Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Formatting, minor edits / updates. 

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@141407 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz 2004-07-23 05:14:43 +00:00
parent 191ec1a14d
commit 207da9b291
1 changed files with 125 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

228
xdocs/userguide/analysis.xml
View File

@ -17,7 +17,7 @@
-->
<?xml-stylesheet type="text/xsl" href="./xdoc.xsl"?>
<!-- $Revision: 1.9 $ $Date: 2004/04/26 20:47:07 $ -->
<!-- $Revision: 1.10 $ $Date: 2004/07/23 05:14:43 $ -->
<document url="analysis.html">
<properties>
<title>The Commons Math User Guide - Numerical Analysis</title>
@ -26,68 +26,68 @@
<section name="4 Numerical Analysis">
<subsection name="4.1 Overview" href="overview">
<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.
The analysis package provides numerical root-finding and interpolation
implementations for real-valued functions of one real variable.
</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.
Possible future additions may include numerical differentation,
integration and optimization.
</p>
</subsection>
<subsection name="4.2 Root-finding" href="rootfinding">
<p>
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. 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.
A <a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealSolver.html">
org.apache.commons.math.analysis.UnivariateRealSolver.</a>
provides the means to find roots of univariate real-valued functions.
A root is the value where the function takes the value 0. Commons-Math
includes implementations of the following root-finding algorithms: <ul>
<li><a href="../apidocs/org/apache/commons/math/analysis/BisectionSolver.html">
Bisection</a></li>
<li><a href="../apidocs/org/apache/commons/math/analysis/BrentSolver.html">
Brent-Dekker</a></li>
<li><a href="../apidocs/org/apache/commons/math/analysis/NewtonSolver.html">
Newton's Method</a></li>
<li><a href="../apidocs/org/apache/commons/math/analysis/SecantSolver.html">
Secant Method</a></li>
</ul>
</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.
There are numerous non-obvious traps and pitfalls in root finding.
First, 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 will not 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, which may be 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
ill-conditioned problems are met. A way around this is to transform
the problem in order to get a better conditioned function. Proper
selection of a root-finding algorithm and its configuration parameters
requires knowledge of the analytical properties of the function under
analysis and numerical analysis techniques. Users are encouraged
to consult a numerical analysis text (or an numerical analyst) when
selecting and configuring a solver.
</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
encouraged that all solver object creation occurs via the <code>org.apache.commons.math.analysis.UnivariateRealSolverFactory</code>
class. <code>UnivariateRealSolverFactory</code> is a simple factory used to create all
of the solver objects supported by Commons-Math. The typical usage of <code>UnivariateRealSolverFactory</code>
In order to use the root-finding features, first a solver object must
be created. It is encouraged that all solver object creation occurs
via the <code>org.apache.commons.math.analysis.UnivariateRealSolverFactory</code>
class. <code>UnivariateRealSolverFactory</code> is a simple factory
used to create all 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>
<p>
The solvers that can be instantiated via the <code>UnivariateRealSolverFactory</code> are detailed below:
The solvers that can be instantiated via the
<code>UnivariateRealSolverFactory</code> are detailed below:
<table>
<tr><th>Solver</th><th>Factory Method</th><th>Notes on Use</th></tr>
<tr><td>Bisection</td><td>newBisectionSolver</td><td><div>Root must be bracketted.</div><div>Linear, guaranteed convergence</div></td></tr>
@ -99,7 +99,7 @@ UnivariateRealSolver solver = factory.newDefaultSolver(function);</source>
<p>
Using a solver object, roots of functions are easily found using the <code>solve</code>
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:
<code>max</code>, <code>solve</code> computes a value <code>c</code> such that:
<ul>
<li><code>f(c) = 0.0</code> (see "function value accuracy")</li>
<li><code>min &lt;= c &lt;= max</code></li>
@ -113,31 +113,39 @@ 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
The <code>BrentSolve</code> uses the Brent-Dekker algorithm which is
fast and robust. This algorithm is recommended for most users and the
<code>BrentSolver</code> is the default solver provided by the
<code>UnivariateRealSolverFactory</code>. 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.
The <code>SecantSolver</code> 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.
The <code>BisectionSolver</code> 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
adverse 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
part, these properties should not have to change from their default values to produce
quality results. In the circumstances where changing these property values is needed, it
is easily done through getter and setter methods on <code>UnivariateRealSolver</code>:
The <code>UnivariateRealSolver</code> interface exposes many
properties to control the convergence of a solver. For the most part,
these properties should not have to change from their default values
to produce good results. In the circumstances where changing these
property values is needed, it 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>
@ -148,13 +156,16 @@ double c = solver.solve(1.0, 5.0);</source>
<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.
The Absolute Accuracy is (estimated) maximal difference between
the computed root and the true root of the function. This is
what most people think of as "accuracy" intuitively. The default
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 default absolute
accuracy may never be reached because the given relative
accuracy is reached first.
</td>
</tr>
<tr>
@ -165,11 +176,13 @@ double c = solver.solve(1.0, 5.0);</source>
<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
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
better 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>
@ -181,11 +194,14 @@ double c = solver.solve(1.0, 5.0);</source>
<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.
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 sped up by setting this value
appropriately.
</td>
</tr>
<tr>
@ -196,11 +212,14 @@ double c = solver.solve(1.0, 5.0);</source>
<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
This is the maximal number of iterations the algorith will try.
If this number is exceeded, non-convergence is assumed and a
<code>ConvergenceException</code> exception is thrown. The
default value is 100, which should be plenty, given 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>
@ -209,14 +228,17 @@ double c = solver.solve(1.0, 5.0);</source>
</subsection>
<subsection name="4.3 Interpolation" href="interpolation">
<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
A <a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealInterpolator.html">
org.apache.commons.math.analysis.UnivariateRealInterpolator</a>
is used to find a univariate real-valued function <code>f</code> which
for a given set of ordered 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
<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.
There aren't currently any accuracy controls either, as interpolation
accuracy is in general determined by the algorithm.
</p>
<p>Typical usage:</p>
<source>double x[]={ 0.0, 1.0, 2.0 };
@ -227,15 +249,15 @@ 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.
A natural cubic spline is a function consisting of a polynominal of
third degree for each subinterval determined by the x-coordinates of the
interpolated points. A function interpolating <code>N</code>
value pairs consists of <code>N-1</code> polynominals. The function
is continuous, smooth and can be differentiated twice. The second
derivative is continuous but not smooth. The x values passed to the
interpolator must be ordered in ascending order. 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>.
</p>
</subsection>
</section>