alvinalexander.com | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Commons Math example source code file (MullerSolver.java)

This example Commons Math source code file (MullerSolver.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Java - Commons Math tags/keywords

convergenceexception, deprecated, deprecated, functionevaluationexception, functionevaluationexception, maxiterationsexceededexception, maxiterationsexceededexception, mullersolver, mullersolver, univariaterealfunction, univariaterealsolverimpl

The Commons Math MullerSolver.java source code

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You 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.solvers;

import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
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, 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.</p>
 *
 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
 * @since 1.2
 */
public class MullerSolver extends UnivariateRealSolverImpl {

    /**
     * Construct a solver for the given function.
     *
     * @param f function to solve
     * @deprecated as of 2.0 the function to solve is passed as an argument
     * to the {@link #solve(UnivariateRealFunction, double, double)} or
     * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
     * method.
     */
    @Deprecated
    public MullerSolver(UnivariateRealFunction f) {
        super(f, 100, 1E-6);
    }

    /**
     * Construct a solver.
     */
    public MullerSolver() {
        super(100, 1E-6);
    }

    /** {@inheritDoc} */
    @Deprecated
    public double solve(final double min, final double max)
        throws ConvergenceException, FunctionEvaluationException {
        return solve(f, min, max);
    }

    /** {@inheritDoc} */
    @Deprecated
    public double solve(final double min, final double max, final double initial)
        throws ConvergenceException, FunctionEvaluationException {
        return solve(f, min, max, initial);
    }

    /**
     * Find a real root in the given interval with initial value.
     * <p>
     * Requires bracketing condition.</p>
     *
     * @param f the function to solve
     * @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 MaxIterationsExceededException 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(final UnivariateRealFunction f,
                        final double min, final double max, final double initial)
        throws MaxIterationsExceededException, 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(f, min, initial);
        } else {
            return solve(f, 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>
     * <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>
     * <p>
     * The formulas here use divided differences directly.</p>
     *
     * @param f the function to solve
     * @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 MaxIterationsExceededException 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(final UnivariateRealFunction f,
                        final double min, final double max)
        throws MaxIterationsExceededException, 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 = min;
        double y0 = f.value(x0);
        double x2 = max;
        double y2 = f.value(x2);
        double x1 = 0.5 * (x0 + x2);
        double 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);

        double oldx = Double.POSITIVE_INFINITY;
        for (int i = 1; i <= maximalIterationCount; ++i) {
            // 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.
            final double d01 = (y1 - y0) / (x1 - x0);
            final double d12 = (y2 - y1) / (x2 - x1);
            final double d012 = (d12 - d01) / (x2 - x0);
            final double c1 = d01 + (x1 - x0) * d012;
            final double delta = c1 * c1 - 4 * y1 * d012;
            final double xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
            final double 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)
            final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
            final double y = f.value(x);

            // check for convergence
            final double 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.
            // the real number equality test x == x1 is intentional and
            // completes the proximity tests above it
            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;
            }
        }
        throw new MaxIterationsExceededException(maximalIterationCount);
    }

    /**
     * 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>
     * <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>
     * <p>
     * The formulas here do not use divided differences directly.</p>
     *
     * @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 MaxIterationsExceededException 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
     * @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)}
     * since 2.0
     */
    @Deprecated
    public double solve2(final double min, final double max)
        throws MaxIterationsExceededException, FunctionEvaluationException {
        return solve2(f, min, max);
    }

    /**
     * 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>
     * <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>
     * <p>
     * The formulas here do not use divided differences directly.</p>
     *
     * @param f the function to solve
     * @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 MaxIterationsExceededException 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(final UnivariateRealFunction f,
                         final double min, final double max)
        throws MaxIterationsExceededException, 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 = min;
        double y0 = f.value(x0);
        double x1 = max;
        double y1 = f.value(x1);
        double x2 = 0.5 * (x0 + x1);
        double 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);

        double oldx = Double.POSITIVE_INFINITY;
        for (int i = 1; i <= maximalIterationCount; ++i) {
            // quadratic interpolation through x0, x1, x2
            final double q = (x2 - x1) / (x1 - x0);
            final double a = q * (y2 - (1 + q) * y1 + q * y0);
            final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
            final double c = (1 + q) * y2;
            final double delta = b * b - 4 * a * c;
            double x;
            final double denominator;
            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 exactly coincides with x1 or x2
                // the equality tests here are intentional
                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;
            }
            final double y = f.value(x);

            // check for convergence
            final double 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;
        }
        throw new MaxIterationsExceededException(maximalIterationCount);
    }
}

Other Commons Math examples (source code examples)

Here is a short list of links related to this Commons Math MullerSolver.java source code file:

... this post is sponsored by my books ...

#1 New Release!

FP Best Seller

 

new blog posts

 

Copyright 1998-2021 Alvin Alexander, alvinalexander.com
All Rights Reserved.

A percentage of advertising revenue from
pages under the /java/jwarehouse URI on this website is
paid back to open source projects.