001: /*
002: * Copyright 2003-2004 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.analysis;
017:
018: import org.apache.commons.math.FunctionEvaluationException;
019: import org.apache.commons.math.ConvergenceException;
020:
021: /**
022: * Utility routines for {@link UnivariateRealSolver} objects.
023: *
024: * @version $Revision: 155427 $ $Date: 2005-02-26 06:11:52 -0700 (Sat, 26 Feb 2005) $
025: */
026: public class UnivariateRealSolverUtils {
027: /**
028: * Default constructor.
029: */
030: private UnivariateRealSolverUtils() {
031: super ();
032: }
033:
034: /** Cached solver factory */
035: private static UnivariateRealSolverFactory factory = null;
036:
037: /**
038: * Convenience method to find a zero of a univariate real function. A default
039: * solver is used.
040: *
041: * @param f the function.
042: * @param x0 the lower bound for the interval.
043: * @param x1 the upper bound for the interval.
044: * @return a value where the function is zero.
045: * @throws ConvergenceException if the iteration count was exceeded
046: * @throws FunctionEvaluationException if an error occurs evaluating
047: * the function
048: * @throws IllegalArgumentException if f is null or the endpoints do not
049: * specify a valid interval
050: */
051: public static double solve(UnivariateRealFunction f, double x0,
052: double x1) throws ConvergenceException,
053: FunctionEvaluationException {
054: setup(f);
055: return factory.newDefaultSolver(f).solve(x0, x1);
056: }
057:
058: /**
059: * Convenience method to find a zero of a univariate real function. A default
060: * solver is used.
061: *
062: * @param f the function
063: * @param x0 the lower bound for the interval
064: * @param x1 the upper bound for the interval
065: * @param absoluteAccuracy the accuracy to be used by the solver
066: * @return a value where the function is zero
067: * @throws ConvergenceException if the iteration count is exceeded
068: * @throws FunctionEvaluationException if an error occurs evaluating the
069: * function
070: * @throws IllegalArgumentException if f is null, the endpoints do not
071: * specify a valid interval, or the absoluteAccuracy is not valid for the
072: * default solver
073: */
074: public static double solve(UnivariateRealFunction f, double x0,
075: double x1, double absoluteAccuracy)
076: throws ConvergenceException, FunctionEvaluationException {
077:
078: setup(f);
079: UnivariateRealSolver solver = factory.newDefaultSolver(f);
080: solver.setAbsoluteAccuracy(absoluteAccuracy);
081: return solver.solve(x0, x1);
082: }
083:
084: /**
085: * This method attempts to find two values a and b satisfying <ul>
086: * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
087: * <li> <code> f(a) * f(b) < 0 </code></li>
088: * </ul>
089: * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
090: * and <code>b</code> bracket a root of f.
091: * <p>
092: * The algorithm starts by setting
093: * <code>a := initial -1; b := initial +1,</code> examines the value of the
094: * function at <code>a</code> and <code>b</code> and keeps moving
095: * the endpoints out by one unit each time through a loop that terminates
096: * when one of the following happens: <ul>
097: * <li> <code> f(a) * f(b) < 0 </code> -- success!</li>
098: * <li> <code> a = lower </code> and <code> b = upper</code>
099: * -- ConvergenceException </li>
100: * <li> <code> Integer.MAX_VALUE</code> iterations elapse
101: * -- ConvergenceException </li>
102: * </ul>
103: * <p>
104: * <strong>Note: </strong> this method can take
105: * <code>Integer.MAX_VALUE</code> iterations to throw a
106: * <code>ConvergenceException.</code> Unless you are confident that there
107: * is a root between <code>lowerBound</code> and <code>upperBound</code>
108: * near <code>initial,</code> it is better to use
109: * {@link #bracket(UnivariateRealFunction, double, double, double, int)},
110: * explicitly specifying the maximum number of iterations.
111: *
112: * @param function the function
113: * @param initial initial midpoint of interval being expanded to
114: * bracket a root
115: * @param lowerBound lower bound (a is never lower than this value)
116: * @param upperBound upper bound (b never is greater than this
117: * value)
118: * @return a two element array holding {a, b}
119: * @throws ConvergenceException if a root can not be bracketted
120: * @throws FunctionEvaluationException if an error occurs evaluating the
121: * function
122: * @throws IllegalArgumentException if function is null, maximumIterations
123: * is not positive, or initial is not between lowerBound and upperBound
124: */
125: public static double[] bracket(UnivariateRealFunction function,
126: double initial, double lowerBound, double upperBound)
127: throws ConvergenceException, FunctionEvaluationException {
128: return bracket(function, initial, lowerBound, upperBound,
129: Integer.MAX_VALUE);
130: }
131:
132: /**
133: * This method attempts to find two values a and b satisfying <ul>
134: * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
135: * <li> <code> f(a) * f(b) < 0 </code> </li>
136: * </ul>
137: * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
138: * and <code>b</code> bracket a root of f.
139: * <p>
140: * The algorithm starts by setting
141: * <code>a := initial -1; b := initial +1,</code> examines the value of the
142: * function at <code>a</code> and <code>b</code> and keeps moving
143: * the endpoints out by one unit each time through a loop that terminates
144: * when one of the following happens: <ul>
145: * <li> <code> f(a) * f(b) < 0 </code> -- success!</li>
146: * <li> <code> a = lower </code> and <code> b = upper</code>
147: * -- ConvergenceException </li>
148: * <li> <code> maximumIterations</code> iterations elapse
149: * -- ConvergenceException </li></ul>
150: *
151: * @param function the function
152: * @param initial initial midpoint of interval being expanded to
153: * bracket a root
154: * @param lowerBound lower bound (a is never lower than this value)
155: * @param upperBound upper bound (b never is greater than this
156: * value)
157: * @param maximumIterations maximum number of iterations to perform
158: * @return a two element array holding {a, b}.
159: * @throws ConvergenceException if the algorithm fails to find a and b
160: * satisfying the desired conditions
161: * @throws FunctionEvaluationException if an error occurs evaluating the
162: * function
163: * @throws IllegalArgumentException if function is null, maximumIterations
164: * is not positive, or initial is not between lowerBound and upperBound
165: */
166: public static double[] bracket(UnivariateRealFunction function,
167: double initial, double lowerBound, double upperBound,
168: int maximumIterations) throws ConvergenceException,
169: FunctionEvaluationException {
170:
171: if (function == null) {
172: throw new IllegalArgumentException("function is null.");
173: }
174: if (maximumIterations <= 0) {
175: throw new IllegalArgumentException(
176: "bad value for maximumIterations: "
177: + maximumIterations);
178: }
179: if (initial < lowerBound || initial > upperBound
180: || lowerBound >= upperBound) {
181: throw new IllegalArgumentException(
182: "Invalid endpoint parameters: lowerBound="
183: + lowerBound + " initial=" + initial
184: + " upperBound=" + upperBound);
185: }
186: double a = initial;
187: double b = initial;
188: double fa;
189: double fb;
190: int numIterations = 0;
191:
192: do {
193: a = Math.max(a - 1.0, lowerBound);
194: b = Math.min(b + 1.0, upperBound);
195: fa = function.value(a);
196:
197: fb = function.value(b);
198: numIterations++;
199: } while ((fa * fb > 0.0) && (numIterations < maximumIterations)
200: && ((a > lowerBound) || (b < upperBound)));
201:
202: if (fa * fb >= 0.0) {
203: throw new ConvergenceException("Number of iterations= "
204: + numIterations + " maximum iterations= "
205: + maximumIterations + " initial= " + initial
206: + " lowerBound=" + lowerBound + " upperBound="
207: + upperBound + " final a value=" + a
208: + " final b value=" + b + " f(a)=" + fa + " f(b)="
209: + fb);
210: }
211:
212: return new double[] { a, b };
213: }
214:
215: /**
216: * Compute the midpoint of two values.
217: *
218: * @param a first value.
219: * @param b second value.
220: * @return the midpoint.
221: */
222: public static double midpoint(double a, double b) {
223: return (a + b) * .5;
224: }
225:
226: /**
227: * Checks to see if f is null, throwing IllegalArgumentException if so.
228: * Also initializes factory if factory is null.
229: *
230: * @param f input function
231: * @throws IllegalArgumentException if f is null
232: */
233: private static void setup(UnivariateRealFunction f) {
234:
235: if (f == null) {
236: throw new IllegalArgumentException(
237: "function can not be null.");
238: }
239:
240: if (factory == null) {
241: factory = UnivariateRealSolverFactory.newInstance();
242: }
243: }
244: }
|