Preliminary checkin of SoC code.
Contributed by: Xiaogang Zhang git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@239940 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Joerg Pietschmann
parent
fd07147f86
commit
69d29afc48
|
|
@ -0,0 +1,119 @@
|
|||
/*
|
||||
* Copyright 2003-2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import java.io.Serializable;
|
||||
import org.apache.commons.math.MathException;
|
||||
|
||||
/**
|
||||
* Implements the <a href="
|
||||
* "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
|
||||
* Divided Difference Algorithm</a> for interpolation of real univariate
|
||||
* functions. For reference, see <b>Introduction to Numerical Analysis</b>,
|
||||
* ISBN 038795452X, chapter 2.
|
||||
* <p>
|
||||
* The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
|
||||
* this class provides an easy-to-use interface to it.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
|
||||
Serializable {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = 107049519551235069L;
|
||||
|
||||
/**
|
||||
* Computes an interpolating function for the data set.
|
||||
*
|
||||
* @param x the interpolating points array
|
||||
* @param y the interpolating values array
|
||||
* @return a function which interpolates the data set
|
||||
* @throws MathException if arguments are invalid
|
||||
*/
|
||||
public UnivariateRealFunction interpolate(double x[], double y[]) throws
|
||||
MathException {
|
||||
|
||||
/**
|
||||
* a[] and c[] are defined in the general formula of Newton form:
|
||||
* p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
|
||||
* a[n](x-c[0])(x-c[1])...(x-c[n-1])
|
||||
*/
|
||||
double a[], c[];
|
||||
|
||||
PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
|
||||
|
||||
/**
|
||||
* When used for interpolation, the Newton form formula becomes
|
||||
* p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
|
||||
* f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
|
||||
* Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
|
||||
* <p>
|
||||
* Note x[], y[], a[] have the same length but c[]'s size is one less.
|
||||
*/
|
||||
c = new double[x.length-1];
|
||||
for (int i = 0; i < c.length; i++) {
|
||||
c[i] = x[i];
|
||||
}
|
||||
a = computeDividedDifference(x, y);
|
||||
|
||||
PolynomialFunctionNewtonForm p;
|
||||
p = new PolynomialFunctionNewtonForm(a, c);
|
||||
return p;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the divided difference array defined recursively by
|
||||
* f[x0] = f(x0)
|
||||
* f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
|
||||
* <p>
|
||||
* The computational complexity is O(N^2).
|
||||
*
|
||||
* @return a fresh copy of the divided difference array
|
||||
* @throws MathException if any abscissas coincide
|
||||
*/
|
||||
protected static double[] computeDividedDifference(double x[], double y[])
|
||||
throws MathException {
|
||||
|
||||
int i, j, n;
|
||||
double divdiff[], a[], denominator;
|
||||
|
||||
PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
|
||||
|
||||
n = x.length;
|
||||
divdiff = new double[n];
|
||||
for (i = 0; i < n; i++) {
|
||||
divdiff[i] = y[i]; // initialization
|
||||
}
|
||||
|
||||
a = new double [n];
|
||||
a[0] = divdiff[0];
|
||||
for (i = 1; i < n; i++) {
|
||||
for (j = 0; j < n-i; j++) {
|
||||
denominator = x[j+i] - x[j];
|
||||
if (denominator == 0.0) {
|
||||
// This happens only when two abscissas are identical.
|
||||
throw new MathException
|
||||
("Identical abscissas cause division by zero.");
|
||||
}
|
||||
divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
|
||||
}
|
||||
a[i] = divdiff[0];
|
||||
}
|
||||
|
||||
return a;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,319 @@
|
|||
/*
|
||||
* Copyright 2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import org.apache.commons.math.ConvergenceException;
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.complex.*;
|
||||
|
||||
/**
|
||||
* Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
|
||||
* Laguerre's Method</a> for root finding of real coefficient polynomials.
|
||||
* For reference, see <b>A First Course in Numerical Analysis</b>,
|
||||
* ISBN 048641454X, chapter 8.
|
||||
* <p>
|
||||
* Laguerre's method is global in the sense that it can start with any initial
|
||||
* approximation and be able to solve all roots from that point.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class LaguerreSolver extends UnivariateRealSolverImpl {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = 5287689975005870178L;
|
||||
|
||||
/** polynomial function to solve */
|
||||
private PolynomialFunction p;
|
||||
|
||||
/**
|
||||
* Construct a solver for the given function.
|
||||
*
|
||||
* @param f function to solve
|
||||
* @throws IllegalArgumentException if function is not polynomial
|
||||
*/
|
||||
public LaguerreSolver(UnivariateRealFunction f) throws
|
||||
IllegalArgumentException {
|
||||
|
||||
super(f, 100, 1E-6);
|
||||
if (f instanceof PolynomialFunction) {
|
||||
p = (PolynomialFunction)f;
|
||||
} else {
|
||||
throw new IllegalArgumentException("Function is not polynomial.");
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the polynomial function.
|
||||
*
|
||||
* @return a fresh copy of the polynomial function
|
||||
*/
|
||||
public PolynomialFunction getPolynomialFunction() {
|
||||
return new PolynomialFunction(p.getCoefficients());
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a real root in the given interval with initial value.
|
||||
* <p>
|
||||
* Requires bracketing condition.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max, double initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (p.value(min) == 0.0) { return min; }
|
||||
if (p.value(max) == 0.0) { return max; }
|
||||
if (p.value(initial) == 0.0) { return initial; }
|
||||
|
||||
verifyBracketing(min, max, p);
|
||||
verifySequence(min, initial, max);
|
||||
if (isBracketing(min, initial, p)) {
|
||||
return solve(min, initial);
|
||||
} else {
|
||||
return solve(initial, max);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a real root in the given interval.
|
||||
* <p>
|
||||
* Despite the bracketing condition, the root returned by solve(Complex[],
|
||||
* Complex) may not be a real zero inside [min, max]. For example,
|
||||
* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
|
||||
* another initial value, or, as we did here, call solveAll() to obtain
|
||||
* all roots and pick up the one that we're looking for.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max) throws ConvergenceException,
|
||||
FunctionEvaluationException {
|
||||
|
||||
Complex c[], root[], initial, z;
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (p.value(min) == 0.0) { return min; }
|
||||
if (p.value(max) == 0.0) { return max; }
|
||||
verifyBracketing(min, max, p);
|
||||
|
||||
double coefficients[] = p.getCoefficients();
|
||||
c = new Complex[coefficients.length];
|
||||
for (int i = 0; i < coefficients.length; i++) {
|
||||
c[i] = new Complex(coefficients[i], 0.0);
|
||||
}
|
||||
initial = new Complex(0.5 * (min + max), 0.0);
|
||||
z = solve(c, initial);
|
||||
if (isRootOK(min, max, z)) {
|
||||
setResult(z.getReal(), iterationCount);
|
||||
return result;
|
||||
}
|
||||
|
||||
// solve all roots and select the one we're seeking
|
||||
root = solveAll(c, initial);
|
||||
for (int i = 0; i < root.length; i++) {
|
||||
if (isRootOK(min, max, root[i])) {
|
||||
setResult(root[i].getReal(), iterationCount);
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
// should never happen
|
||||
throw new ConvergenceException("Convergence failed.");
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns true iff the given complex root is actually a real zero
|
||||
* in the given interval, within the solver tolerance level.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @param z the complex root
|
||||
* @return true iff z is the sought-after real zero
|
||||
*/
|
||||
protected boolean isRootOK(double min, double max, Complex z) {
|
||||
double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
|
||||
return (isSequence(min, z.getReal(), max)) &&
|
||||
(Math.abs(z.getImaginary()) <= tolerance ||
|
||||
z.abs() <= functionValueAccuracy);
|
||||
}
|
||||
|
||||
/**
|
||||
* Find all complex roots for the polynomial with the given coefficients,
|
||||
* starting from the given initial value.
|
||||
*
|
||||
* @param coefficients the polynomial coefficients array
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public Complex[] solveAll(double coefficients[], double initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
Complex c[] = new Complex[coefficients.length];
|
||||
Complex z = new Complex(initial, 0.0);
|
||||
for (int i = 0; i < c.length; i++) {
|
||||
c[i] = new Complex(coefficients[i], 0.0);
|
||||
}
|
||||
return solveAll(c, z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Find all complex roots for the polynomial with the given coefficients,
|
||||
* starting from the given initial value.
|
||||
*
|
||||
* @param coefficients the polynomial coefficients array
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public Complex[] solveAll(Complex coefficients[], Complex initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
int i, j, n, iterationCount = 0;
|
||||
Complex root[], c[], subarray[], oldc, newc;
|
||||
|
||||
n = coefficients.length - 1;
|
||||
if (n < 1) {
|
||||
throw new IllegalArgumentException
|
||||
("Polynomial degree must be positive: degree=" + n);
|
||||
}
|
||||
c = new Complex[n+1]; // coefficients for deflated polynomial
|
||||
for (i = 0; i <= n; i++) {
|
||||
c[i] = coefficients[i];
|
||||
}
|
||||
|
||||
// solve individual root successively
|
||||
root = new Complex[n];
|
||||
for (i = 0; i < n; i++) {
|
||||
subarray = new Complex[n-i+1];
|
||||
System.arraycopy(c, 0, subarray, 0, subarray.length);
|
||||
root[i] = solve(subarray, initial);
|
||||
// polynomial deflation using synthetic division
|
||||
newc = c[n-i];
|
||||
for (j = n-i-1; j >= 0; j--) {
|
||||
oldc = c[j];
|
||||
c[j] = newc;
|
||||
newc = oldc.add(newc.multiply(root[i]));
|
||||
}
|
||||
iterationCount += this.iterationCount;
|
||||
}
|
||||
|
||||
resultComputed = true;
|
||||
this.iterationCount = iterationCount;
|
||||
return root;
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a complex root for the polynomial with the given coefficients,
|
||||
* starting from the given initial value.
|
||||
*
|
||||
* @param coefficients the polynomial coefficients array
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public Complex solve(Complex coefficients[], Complex initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
Complex z = initial, oldz, pv, dv, d2v, G, G2, H, delta, denominator;
|
||||
|
||||
int n = coefficients.length - 1;
|
||||
if (n < 1) {
|
||||
throw new IllegalArgumentException
|
||||
("Polynomial degree must be positive: degree=" + n);
|
||||
}
|
||||
Complex N = new Complex((double)n, 0.0);
|
||||
Complex N1 = new Complex((double)(n-1), 0.0);
|
||||
|
||||
int i = 1;
|
||||
oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
|
||||
while (i <= maximalIterationCount) {
|
||||
// Compute pv (polynomial value), dv (derivative value), and
|
||||
// d2v (second derivative value) simultaneously.
|
||||
pv = coefficients[n];
|
||||
dv = d2v = new Complex(0.0, 0.0);
|
||||
for (int j = n-1; j >= 0; j--) {
|
||||
d2v = dv.add(z.multiply(d2v));
|
||||
dv = pv.add(z.multiply(dv));
|
||||
pv = coefficients[j].add(z.multiply(pv));
|
||||
}
|
||||
d2v = d2v.multiply(new Complex(2.0, 0.0));
|
||||
|
||||
// check for convergence
|
||||
double tolerance = Math.max(relativeAccuracy * z.abs(),
|
||||
absoluteAccuracy);
|
||||
if ((z.subtract(oldz)).abs() <= tolerance) {
|
||||
resultComputed = true;
|
||||
iterationCount = i;
|
||||
return z;
|
||||
}
|
||||
if (pv.abs() <= functionValueAccuracy) {
|
||||
resultComputed = true;
|
||||
iterationCount = i;
|
||||
return z;
|
||||
}
|
||||
|
||||
// now pv != 0, calculate the new approximation
|
||||
G = dv.divide(pv);
|
||||
G2 = G.multiply(G);
|
||||
H = G2.subtract(d2v.divide(pv));
|
||||
delta = N1.multiply((N.multiply(H)).subtract(G2));
|
||||
// choose a denominator larger in magnitude
|
||||
Complex dplus = G.add(ComplexUtils.sqrt(delta));
|
||||
Complex dminus = G.subtract(ComplexUtils.sqrt(delta));
|
||||
denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
|
||||
// Perturb z if denominator is zero, for instance,
|
||||
// p(x) = x^3 + 1, z = 0.
|
||||
if (denominator.equals(new Complex(0.0, 0.0))) {
|
||||
z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
|
||||
oldz = new Complex(Double.POSITIVE_INFINITY,
|
||||
Double.POSITIVE_INFINITY);
|
||||
} else {
|
||||
oldz = z;
|
||||
z = z.subtract(N.divide(denominator));
|
||||
}
|
||||
i++;
|
||||
}
|
||||
throw new ConvergenceException("Maximum number of iterations exceeded.");
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,277 @@
|
|||
/*
|
||||
* Copyright 2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import org.apache.commons.math.ConvergenceException;
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
* Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
|
||||
* Muller's Method</a> for root finding of real univariate functions. For
|
||||
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
|
||||
* chapter 3.
|
||||
* <p>
|
||||
* Muller's method applies to both real and complex functions, but here we
|
||||
* restrict ourselves to real functions. Methods solve() and solve2() find
|
||||
* real zeros, using different ways to bypass complex arithmetics.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class MullerSolver extends UnivariateRealSolverImpl {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = 2619993603551148137L;
|
||||
|
||||
/**
|
||||
* Construct a solver for the given function.
|
||||
*
|
||||
* @param f function to solve
|
||||
*/
|
||||
public MullerSolver(UnivariateRealFunction f) {
|
||||
super(f, 100, 1E-6);
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a real root in the given interval with initial value.
|
||||
* <p>
|
||||
* Requires bracketing condition.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max, double initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (f.value(min) == 0.0) { return min; }
|
||||
if (f.value(max) == 0.0) { return max; }
|
||||
if (f.value(initial) == 0.0) { return initial; }
|
||||
|
||||
verifyBracketing(min, max, f);
|
||||
verifySequence(min, initial, max);
|
||||
if (isBracketing(min, initial, f)) {
|
||||
return solve(min, initial);
|
||||
} else {
|
||||
return solve(initial, max);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a real root in the given interval.
|
||||
* <p>
|
||||
* Original Muller's method would have function evaluation at complex point.
|
||||
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
|
||||
* condition is one way to go: by requiring bracketing in every iteration,
|
||||
* the newly computed approximation is guaranteed to be real.
|
||||
* <p>
|
||||
* Normally Muller's method converges quadratically in the vicinity of a
|
||||
* zero, however it may be very slow in regions far away from zeros. For
|
||||
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
|
||||
* bisection as a safety backup if it performs very poorly.
|
||||
* <p>
|
||||
* The formulas here use divided differences directly.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max) throws ConvergenceException,
|
||||
FunctionEvaluationException {
|
||||
|
||||
// [x0, x2] is the bracketing interval in each iteration
|
||||
// x1 is the last approximation and an interpolation point in (x0, x2)
|
||||
// x is the new root approximation and new x1 for next round
|
||||
// d01, d12, d012 are divided differences
|
||||
double x0, x1, x2, x, oldx, y0, y1, y2, y;
|
||||
double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
|
||||
|
||||
x0 = min; y0 = f.value(x0);
|
||||
x2 = max; y2 = f.value(x2);
|
||||
x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (y0 == 0.0) { return min; }
|
||||
if (y2 == 0.0) { return max; }
|
||||
verifyBracketing(min, max, f);
|
||||
|
||||
int i = 1;
|
||||
oldx = Double.POSITIVE_INFINITY;
|
||||
while (i <= maximalIterationCount) {
|
||||
// Muller's method employs quadratic interpolation through
|
||||
// x0, x1, x2 and x is the zero of the interpolating parabola.
|
||||
// Due to bracketing condition, this parabola must have two
|
||||
// real roots and we choose one in [x0, x2] to be x.
|
||||
d01 = (y1 - y0) / (x1 - x0);
|
||||
d12 = (y2 - y1) / (x2 - x1);
|
||||
d012 = (d12 - d01) / (x2 - x0);
|
||||
c1 = d01 + (x1 - x0) * d012;
|
||||
delta = c1 * c1 - 4 * y1 * d012;
|
||||
xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
|
||||
xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
|
||||
// xplus and xminus are two roots of parabola and at least
|
||||
// one of them should lie in (x0, x2)
|
||||
x = isSequence(x0, xplus, x2) ? xplus : xminus;
|
||||
y = f.value(x);
|
||||
|
||||
// check for convergence
|
||||
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
|
||||
if (Math.abs(x - oldx) <= tolerance) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
if (Math.abs(y) <= functionValueAccuracy) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
|
||||
// Bisect if convergence is too slow. Bisection would waste
|
||||
// our calculation of x, hopefully it won't happen often.
|
||||
boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
|
||||
(x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
|
||||
(x == x1);
|
||||
// prepare the new bracketing interval for next iteration
|
||||
if (!bisect) {
|
||||
x0 = x < x1 ? x0 : x1;
|
||||
y0 = x < x1 ? y0 : y1;
|
||||
x2 = x > x1 ? x2 : x1;
|
||||
y2 = x > x1 ? y2 : y1;
|
||||
x1 = x; y1 = y;
|
||||
oldx = x;
|
||||
} else {
|
||||
double xm = 0.5 * (x0 + x2);
|
||||
double ym = f.value(xm);
|
||||
if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
|
||||
x2 = xm; y2 = ym;
|
||||
} else {
|
||||
x0 = xm; y0 = ym;
|
||||
}
|
||||
x1 = 0.5 * (x0 + x2);
|
||||
y1 = f.value(x1);
|
||||
oldx = Double.POSITIVE_INFINITY;
|
||||
}
|
||||
i++;
|
||||
}
|
||||
throw new ConvergenceException("Maximum number of iterations exceeded.");
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a real root in the given interval.
|
||||
* <p>
|
||||
* solve2() differs from solve() in the way it avoids complex operations.
|
||||
* Except for the initial [min, max], solve2() does not require bracketing
|
||||
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
|
||||
* number arises in the computation, we simply use its modulus as real
|
||||
* approximation.
|
||||
* <p>
|
||||
* Because the interval may not be bracketing, bisection alternative is
|
||||
* not applicable here. However in practice our treatment usually works
|
||||
* well, especially near real zeros where the imaginary part of complex
|
||||
* approximation is often negligible.
|
||||
* <p>
|
||||
* The formulas here do not use divided differences directly.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve2(double min, double max) throws ConvergenceException,
|
||||
FunctionEvaluationException {
|
||||
|
||||
// x2 is the last root approximation
|
||||
// x is the new approximation and new x2 for next round
|
||||
// x0 < x1 < x2 does not hold here
|
||||
double x0, x1, x2, x, oldx, y0, y1, y2, y;
|
||||
double q, A, B, C, delta, denominator, tolerance;
|
||||
|
||||
x0 = min; y0 = f.value(x0);
|
||||
x1 = max; y1 = f.value(x1);
|
||||
x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (y0 == 0.0) { return min; }
|
||||
if (y1 == 0.0) { return max; }
|
||||
verifyBracketing(min, max, f);
|
||||
|
||||
int i = 1;
|
||||
oldx = Double.POSITIVE_INFINITY;
|
||||
while (i <= maximalIterationCount) {
|
||||
// quadratic interpolation through x0, x1, x2
|
||||
q = (x2 - x1) / (x1 - x0);
|
||||
A = q * (y2 - (1 + q) * y1 + q * y0);
|
||||
B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
|
||||
C = (1 + q) * y2;
|
||||
delta = B * B - 4 * A * C;
|
||||
if (delta >= 0.0) {
|
||||
// choose a denominator larger in magnitude
|
||||
double dplus = B + Math.sqrt(delta);
|
||||
double dminus = B - Math.sqrt(delta);
|
||||
denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
|
||||
} else {
|
||||
// take the modulus of (B +/- Math.sqrt(delta))
|
||||
denominator = Math.sqrt(B * B - delta);
|
||||
}
|
||||
if (denominator != 0) {
|
||||
x = x2 - 2.0 * C * (x2 - x1) / denominator;
|
||||
// perturb x if it coincides with x1 or x2
|
||||
while (x == x1 || x == x2) {
|
||||
x += absoluteAccuracy;
|
||||
}
|
||||
} else {
|
||||
// extremely rare case, get a random number to skip it
|
||||
x = min + Math.random() * (max - min);
|
||||
oldx = Double.POSITIVE_INFINITY;
|
||||
}
|
||||
y = f.value(x);
|
||||
|
||||
// check for convergence
|
||||
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
|
||||
if (Math.abs(x - oldx) <= tolerance) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
if (Math.abs(y) <= functionValueAccuracy) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
|
||||
// prepare the next iteration
|
||||
x0 = x1; y0 = y1;
|
||||
x1 = x2; y1 = y2;
|
||||
x2 = x; y2 = y;
|
||||
oldx = x;
|
||||
i++;
|
||||
}
|
||||
throw new ConvergenceException("Maximum number of iterations exceeded.");
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,53 @@
|
|||
/*
|
||||
* Copyright 2003-2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import java.io.Serializable;
|
||||
import org.apache.commons.math.MathException;
|
||||
|
||||
/**
|
||||
* Implements the <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
|
||||
* Neville's Algorithm</a> for interpolation of real univariate functions. For
|
||||
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
|
||||
* chapter 2.
|
||||
* <p>
|
||||
* The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
|
||||
* this class provides an easy-to-use interface to it.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class NevilleInterpolator implements UnivariateRealInterpolator,
|
||||
Serializable {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = 3003707660147873733L;
|
||||
|
||||
/**
|
||||
* Computes an interpolating function for the data set.
|
||||
*
|
||||
* @param x the interpolating points array
|
||||
* @param y the interpolating values array
|
||||
* @return a function which interpolates the data set
|
||||
* @throws MathException if arguments are invalid
|
||||
*/
|
||||
public UnivariateRealFunction interpolate(double x[], double y[]) throws
|
||||
MathException {
|
||||
|
||||
PolynomialFunctionLagrangeForm p;
|
||||
p = new PolynomialFunctionLagrangeForm(x, y);
|
||||
return p;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,290 @@
|
|||
/*
|
||||
* Copyright 2003-2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import java.io.Serializable;
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
|
||||
/**
|
||||
* Implements the representation of a real polynomial function in
|
||||
* <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
|
||||
* Lagrange Form</a>. For reference, see <b>Introduction to Numerical
|
||||
* Analysis</b>, ISBN 038795452X, chapter 2.
|
||||
* <p>
|
||||
* The approximated function should be smooth enough for Lagrange polynomial
|
||||
* to work well. Otherwise, consider using splines instead.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
|
||||
Serializable {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = -3965199246151093920L;
|
||||
|
||||
/**
|
||||
* The coefficients of the polynomial, ordered by degree -- i.e.
|
||||
* coefficients[0] is the constant term and coefficients[n] is the
|
||||
* coefficient of x^n where n is the degree of the polynomial.
|
||||
*/
|
||||
private double coefficients[];
|
||||
|
||||
/**
|
||||
* Interpolating points (abscissas) and the function values at these points.
|
||||
*/
|
||||
private double x[], y[];
|
||||
|
||||
/**
|
||||
* Whether the polynomial coefficients are available.
|
||||
*/
|
||||
private boolean coefficientsComputed;
|
||||
|
||||
/**
|
||||
* Construct a Lagrange polynomial with the given abscissas and function
|
||||
* values. The order of interpolating points are not important.
|
||||
* <p>
|
||||
* The constructor makes copy of the input arrays and assigns them.
|
||||
*
|
||||
* @param x interpolating points
|
||||
* @param y function values at interpolating points
|
||||
* @throws IllegalArgumentException if input arrays are not valid
|
||||
*/
|
||||
PolynomialFunctionLagrangeForm(double x[], double y[]) throws
|
||||
IllegalArgumentException {
|
||||
|
||||
verifyInterpolationArray(x, y);
|
||||
this.x = new double[x.length];
|
||||
this.y = new double[y.length];
|
||||
System.arraycopy(x, 0, this.x, 0, x.length);
|
||||
System.arraycopy(y, 0, this.y, 0, y.length);
|
||||
coefficientsComputed = false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate the function value at the given point.
|
||||
*
|
||||
* @param z the point at which the function value is to be computed
|
||||
* @return the function value
|
||||
* @throws FunctionEvaluationException if a runtime error occurs
|
||||
* @see UnivariateRealFunction#value(double)
|
||||
*/
|
||||
public double value(double z) throws FunctionEvaluationException {
|
||||
return evaluate(x, y, z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the degree of the polynomial.
|
||||
*
|
||||
* @return the degree of the polynomial
|
||||
*/
|
||||
public int degree() {
|
||||
return x.length - 1;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the interpolating points array.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of the interpolating points array
|
||||
*/
|
||||
public double[] getInterpolatingPoints() {
|
||||
double[] out = new double[x.length];
|
||||
System.arraycopy(x, 0, out, 0, x.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the interpolating values array.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of the interpolating values array
|
||||
*/
|
||||
public double[] getInterpolatingValues() {
|
||||
double[] out = new double[y.length];
|
||||
System.arraycopy(y, 0, out, 0, y.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the coefficients array.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of the coefficients array
|
||||
*/
|
||||
public double[] getCoefficients() {
|
||||
if (!coefficientsComputed) {
|
||||
computeCoefficients();
|
||||
}
|
||||
double[] out = new double[coefficients.length];
|
||||
System.arraycopy(coefficients, 0, out, 0, coefficients.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Evaluate the Lagrange polynomial using
|
||||
* <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
|
||||
* Neville's Algorithm</a>. It takes O(N^2) time.
|
||||
* <p>
|
||||
* This function is made public static so that users can call it directly
|
||||
* without instantiating PolynomialFunctionLagrangeForm object.
|
||||
*
|
||||
* @param x the interpolating points array
|
||||
* @param y the interpolating values array
|
||||
* @param z the point at which the function value is to be computed
|
||||
* @return the function value
|
||||
* @throws FunctionEvaluationException if a runtime error occurs
|
||||
* @throws IllegalArgumentException if inputs are not valid
|
||||
*/
|
||||
public static double evaluate(double x[], double y[], double z) throws
|
||||
FunctionEvaluationException, IllegalArgumentException {
|
||||
|
||||
int i, j, n, nearest = 0;
|
||||
double value, c[], d[], tc, td, divider, w, dist, min_dist;
|
||||
|
||||
verifyInterpolationArray(x, y);
|
||||
|
||||
n = x.length;
|
||||
c = new double[n];
|
||||
d = new double[n];
|
||||
min_dist = Double.POSITIVE_INFINITY;
|
||||
for (i = 0; i < n; i++) {
|
||||
// initialize the difference arrays
|
||||
c[i] = y[i];
|
||||
d[i] = y[i];
|
||||
// find out the abscissa closest to z
|
||||
dist = Math.abs(z - x[i]);
|
||||
if (dist < min_dist) {
|
||||
nearest = i;
|
||||
min_dist = dist;
|
||||
}
|
||||
}
|
||||
|
||||
// initial approximation to the function value at z
|
||||
value = y[nearest];
|
||||
|
||||
for (i = 1; i < n; i++) {
|
||||
for (j = 0; j < n-i; j++) {
|
||||
tc = x[j] - z;
|
||||
td = x[i+j] - z;
|
||||
divider = x[j] - x[i+j];
|
||||
if (divider == 0.0) {
|
||||
// This happens only when two abscissas are identical.
|
||||
throw new FunctionEvaluationException(z,
|
||||
"Identical abscissas cause division by zero: x[" +
|
||||
i + "] = x[" + (i+j) + "] = " + x[i]);
|
||||
}
|
||||
// update the difference arrays
|
||||
w = (c[j+1] - d[j]) / divider;
|
||||
c[j] = tc * w;
|
||||
d[j] = td * w;
|
||||
}
|
||||
// sum up the difference terms to get the final value
|
||||
if (nearest < 0.5*(n-i+1)) {
|
||||
value += c[nearest]; // fork down
|
||||
} else {
|
||||
nearest--;
|
||||
value += d[nearest]; // fork up
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate the coefficients of Lagrange polynomial from the
|
||||
* interpolation data. It takes O(N^2) time.
|
||||
* <p>
|
||||
* Note this computation can be ill-conditioned. Use with caution
|
||||
* and only when it is necessary.
|
||||
*
|
||||
* @throws ArithmeticException if any abscissas coincide
|
||||
*/
|
||||
protected void computeCoefficients() throws ArithmeticException {
|
||||
int i, j, n;
|
||||
double c[], tc[], d, t;
|
||||
|
||||
n = degree() + 1;
|
||||
coefficients = new double[n];
|
||||
for (i = 0; i < n; i++) {
|
||||
coefficients[i] = 0.0;
|
||||
}
|
||||
|
||||
// c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
|
||||
c = new double[n+1];
|
||||
c[0] = 1.0;
|
||||
for (i = 0; i < n; i++) {
|
||||
for (j = i; j > 0; j--) {
|
||||
c[j] = c[j-1] - c[j] * x[i];
|
||||
}
|
||||
c[0] *= (-x[i]);
|
||||
c[i+1] = 1;
|
||||
}
|
||||
|
||||
tc = new double[n];
|
||||
for (i = 0; i < n; i++) {
|
||||
// d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
|
||||
d = 1;
|
||||
for (j = 0; j < n; j++) {
|
||||
if (i != j) {
|
||||
d *= (x[i] - x[j]);
|
||||
}
|
||||
}
|
||||
if (d == 0.0) {
|
||||
// This happens only when two abscissas are identical.
|
||||
throw new ArithmeticException
|
||||
("Identical abscissas cause division by zero.");
|
||||
}
|
||||
t = y[i] / d;
|
||||
// Lagrange polynomial is the sum of n terms, each of which is a
|
||||
// polynomial of degree n-1. tc[] are the coefficients of the i-th
|
||||
// numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
|
||||
tc[n-1] = c[n]; // actually c[n] = 1
|
||||
coefficients[n-1] += t * tc[n-1];
|
||||
for (j = n-2; j >= 0; j--) {
|
||||
tc[j] = c[j+1] + tc[j+1] * x[i];
|
||||
coefficients[j] += t * tc[j];
|
||||
}
|
||||
}
|
||||
|
||||
coefficientsComputed = true;
|
||||
}
|
||||
|
||||
/**
|
||||
* Verifies that the interpolation arrays are valid.
|
||||
* <p>
|
||||
* The interpolating points must be distinct. However it is not
|
||||
* verified here, it is checked in evaluate() and computeCoefficients().
|
||||
*
|
||||
* @throws IllegalArgumentException if not valid
|
||||
* @see #evaluate(double[], double[], double)
|
||||
* @see #computeCoefficients()
|
||||
*/
|
||||
protected static void verifyInterpolationArray(double x[], double y[]) throws
|
||||
IllegalArgumentException {
|
||||
|
||||
if (x.length < 2 || y.length < 2) {
|
||||
throw new IllegalArgumentException
|
||||
("Interpolation requires at least two points.");
|
||||
}
|
||||
if (x.length != y.length) {
|
||||
throw new IllegalArgumentException
|
||||
("Abscissa and value arrays must have the same length.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,214 @@
|
|||
/*
|
||||
* Copyright 2003-2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import java.io.Serializable;
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
|
||||
/**
|
||||
* Implements the representation of a real polynomial function in
|
||||
* Newton Form. For reference, see <b>Elementary Numerical Analysis</b>,
|
||||
* ISBN 0070124477, chapter 2.
|
||||
* <p>
|
||||
* The formula of polynomial in Newton form is
|
||||
* p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
|
||||
* a[n](x-c[0])(x-c[1])...(x-c[n-1])
|
||||
* Note that the length of a[] is one more than the length of c[]
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class PolynomialFunctionNewtonForm implements UnivariateRealFunction,
|
||||
Serializable {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = -3353896576191389897L;
|
||||
|
||||
/**
|
||||
* The coefficients of the polynomial, ordered by degree -- i.e.
|
||||
* coefficients[0] is the constant term and coefficients[n] is the
|
||||
* coefficient of x^n where n is the degree of the polynomial.
|
||||
*/
|
||||
private double coefficients[];
|
||||
|
||||
/**
|
||||
* Members of c[] are called centers of the Newton polynomial.
|
||||
* When all c[i] = 0, a[] becomes normal polynomial coefficients,
|
||||
* i.e. a[i] = coefficients[i].
|
||||
*/
|
||||
private double a[], c[];
|
||||
|
||||
/**
|
||||
* Whether the polynomial coefficients are available.
|
||||
*/
|
||||
private boolean coefficientsComputed;
|
||||
|
||||
/**
|
||||
* Construct a Newton polynomial with the given a[] and c[]. The order of
|
||||
* centers are important in that if c[] shuffle, then values of a[] would
|
||||
* completely change, not just a permutation of old a[].
|
||||
* <p>
|
||||
* The constructor makes copy of the input arrays and assigns them.
|
||||
*
|
||||
* @param a the coefficients in Newton form formula
|
||||
* @param c the centers
|
||||
* @throws IllegalArgumentException if input arrays are not valid
|
||||
*/
|
||||
PolynomialFunctionNewtonForm(double a[], double c[]) throws
|
||||
IllegalArgumentException {
|
||||
|
||||
verifyInputArray(a, c);
|
||||
this.a = new double[a.length];
|
||||
this.c = new double[c.length];
|
||||
System.arraycopy(a, 0, this.a, 0, a.length);
|
||||
System.arraycopy(c, 0, this.c, 0, c.length);
|
||||
coefficientsComputed = false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate the function value at the given point.
|
||||
*
|
||||
* @param z the point at which the function value is to be computed
|
||||
* @return the function value
|
||||
* @throws FunctionEvaluationException if a runtime error occurs
|
||||
* @see UnivariateRealFunction#value(double)
|
||||
*/
|
||||
public double value(double z) throws FunctionEvaluationException {
|
||||
return evaluate(a, c, z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the degree of the polynomial.
|
||||
*
|
||||
* @return the degree of the polynomial
|
||||
*/
|
||||
public int degree() {
|
||||
return c.length;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of coefficients in Newton form formula.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of coefficients in Newton form formula
|
||||
*/
|
||||
public double[] getNewtonCoefficients() {
|
||||
double[] out = new double[a.length];
|
||||
System.arraycopy(a, 0, out, 0, a.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the centers array.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of the centers array
|
||||
*/
|
||||
public double[] getCenters() {
|
||||
double[] out = new double[c.length];
|
||||
System.arraycopy(c, 0, out, 0, c.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the coefficients array.
|
||||
* <p>
|
||||
* Changes made to the returned copy will not affect the polynomial.
|
||||
*
|
||||
* @return a fresh copy of the coefficients array
|
||||
*/
|
||||
public double[] getCoefficients() {
|
||||
if (!coefficientsComputed) {
|
||||
computeCoefficients();
|
||||
}
|
||||
double[] out = new double[coefficients.length];
|
||||
System.arraycopy(coefficients, 0, out, 0, coefficients.length);
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* Evaluate the Newton polynomial using nested multiplication. It is
|
||||
* also called <a href="http://mathworld.wolfram.com/HornersRule.html">
|
||||
* Horner's Rule</a> and takes O(N) time.
|
||||
*
|
||||
* @param a the coefficients in Newton form formula
|
||||
* @param c the centers
|
||||
* @param z the point at which the function value is to be computed
|
||||
* @return the function value
|
||||
* @throws FunctionEvaluationException if a runtime error occurs
|
||||
* @throws IllegalArgumentException if inputs are not valid
|
||||
*/
|
||||
public static double evaluate(double a[], double c[], double z) throws
|
||||
FunctionEvaluationException, IllegalArgumentException {
|
||||
|
||||
verifyInputArray(a, c);
|
||||
|
||||
int n = c.length;
|
||||
double value = a[n];
|
||||
for (int i = n-1; i >= 0; i--) {
|
||||
value = a[i] + (z - c[i]) * value;
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate the normal polynomial coefficients given the Newton form.
|
||||
* It also uses nested multiplication but takes O(N^2) time.
|
||||
*/
|
||||
protected void computeCoefficients() {
|
||||
int i, j, n = degree();
|
||||
|
||||
coefficients = new double[n+1];
|
||||
for (i = 0; i <= n; i++) {
|
||||
coefficients[i] = 0.0;
|
||||
}
|
||||
|
||||
coefficients[0] = a[n];
|
||||
for (i = n-1; i >= 0; i--) {
|
||||
for (j = n-i; j > 0; j--) {
|
||||
coefficients[j] = coefficients[j-1] - c[i] * coefficients[j];
|
||||
}
|
||||
coefficients[0] = a[i] - c[i] * coefficients[0];
|
||||
}
|
||||
|
||||
coefficientsComputed = true;
|
||||
}
|
||||
|
||||
/**
|
||||
* Verifies that the input arrays are valid.
|
||||
* <p>
|
||||
* The centers must be distinct for interpolation purposes, but not
|
||||
* for general use. Thus it is not verified here.
|
||||
*
|
||||
* @throws IllegalArgumentException if not valid
|
||||
* @see DividedDifferenceInterpolator#computeDividedDifference(double[],
|
||||
* double[])
|
||||
*/
|
||||
protected static void verifyInputArray(double a[], double c[]) throws
|
||||
IllegalArgumentException {
|
||||
|
||||
if (a.length < 1 || c.length < 1) {
|
||||
throw new IllegalArgumentException
|
||||
("Input arrays must not be empty.");
|
||||
}
|
||||
if (a.length != c.length + 1) {
|
||||
throw new IllegalArgumentException
|
||||
("Bad input array sizes, should have difference 1.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,158 @@
|
|||
/*
|
||||
* Copyright 2005 The Apache Software Foundation.
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.analysis;
|
||||
|
||||
import org.apache.commons.math.ConvergenceException;
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
* Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
|
||||
* Ridders' Method</a> for root finding of real univariate functions. For
|
||||
* reference, see C. Ridders, <i>A new algorithm for computing a single root
|
||||
* of a real continuous function </i>, IEEE Transactions on Circuits and
|
||||
* Systems, 26 (1979), 979 - 980.
|
||||
* <p>
|
||||
* The function should be continuous but not necessarily smooth.
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class RiddersSolver extends UnivariateRealSolverImpl {
|
||||
|
||||
/** serializable version identifier */
|
||||
static final long serialVersionUID = -4703139035737911735L;
|
||||
|
||||
/**
|
||||
* Construct a solver for the given function.
|
||||
*
|
||||
* @param f function to solve
|
||||
*/
|
||||
public RiddersSolver(UnivariateRealFunction f) {
|
||||
super(f, 100, 1E-6);
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a root in the given interval with initial value.
|
||||
* <p>
|
||||
* Requires bracketing condition.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @param initial the start value to use
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max, double initial) throws
|
||||
ConvergenceException, FunctionEvaluationException {
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (f.value(min) == 0.0) { return min; }
|
||||
if (f.value(max) == 0.0) { return max; }
|
||||
if (f.value(initial) == 0.0) { return initial; }
|
||||
|
||||
verifyBracketing(min, max, f);
|
||||
verifySequence(min, initial, max);
|
||||
if (isBracketing(min, initial, f)) {
|
||||
return solve(min, initial);
|
||||
} else {
|
||||
return solve(initial, max);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Find a root in the given interval.
|
||||
* <p>
|
||||
* Requires bracketing condition.
|
||||
*
|
||||
* @param min the lower bound for the interval
|
||||
* @param max the upper bound for the interval
|
||||
* @return the point at which the function value is zero
|
||||
* @throws ConvergenceException if the maximum iteration count is exceeded
|
||||
* or the solver detects convergence problems otherwise
|
||||
* @throws FunctionEvaluationException if an error occurs evaluating the
|
||||
* function
|
||||
* @throws IllegalArgumentException if any parameters are invalid
|
||||
*/
|
||||
public double solve(double min, double max) throws ConvergenceException,
|
||||
FunctionEvaluationException {
|
||||
|
||||
// [x1, x2] is the bracketing interval in each iteration
|
||||
// x3 is the midpoint of [x1, x2]
|
||||
// x is the new root approximation and an endpoint of the new interval
|
||||
double x1, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance;
|
||||
|
||||
x1 = min; y1 = f.value(x1);
|
||||
x2 = max; y2 = f.value(x2);
|
||||
|
||||
// check for zeros before verifying bracketing
|
||||
if (y1 == 0.0) { return min; }
|
||||
if (y2 == 0.0) { return max; }
|
||||
verifyBracketing(min, max, f);
|
||||
|
||||
int i = 1;
|
||||
oldx = Double.POSITIVE_INFINITY;
|
||||
while (i <= maximalIterationCount) {
|
||||
// calculate the new root approximation
|
||||
x3 = 0.5 * (x1 + x2);
|
||||
y3 = f.value(x3);
|
||||
if (Math.abs(y3) <= functionValueAccuracy) {
|
||||
setResult(x3, i);
|
||||
return result;
|
||||
}
|
||||
delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
|
||||
correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
|
||||
(x3 - x1) / Math.sqrt(delta);
|
||||
x = x3 - correction; // correction != 0
|
||||
y = f.value(x);
|
||||
|
||||
// check for convergence
|
||||
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
|
||||
if (Math.abs(x - oldx) <= tolerance) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
if (Math.abs(y) <= functionValueAccuracy) {
|
||||
setResult(x, i);
|
||||
return result;
|
||||
}
|
||||
|
||||
// prepare the new interval for next iteration
|
||||
// Ridders' method guarantees x1 < x < x2
|
||||
if (correction > 0.0) { // x1 < x < x3
|
||||
if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
|
||||
x2 = x; y2 = y;
|
||||
} else {
|
||||
x1 = x; x2 = x3;
|
||||
y1 = y; y2 = y3;
|
||||
}
|
||||
} else { // x3 < x < x2
|
||||
if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
|
||||
x1 = x; y1 = y;
|
||||
} else {
|
||||
x1 = x3; x2 = x;
|
||||
y1 = y3; y2 = y;
|
||||
}
|
||||
}
|
||||
oldx = x;
|
||||
i++;
|
||||
}
|
||||
throw new ConvergenceException("Maximum number of iterations exceeded.");
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue