commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Mikkel Meyer Andersen <m...@mikl.dk>
Subject Re: [math] Accuracy in root finding with newBrentSolver
Date Sun, 15 Aug 2010 13:29:52 GMT
Hi,

Isn't it possible to use BigDecimal or similar (e.g.
http://www.apfloat.org)? I know that some code might need to be
rewritten.

Cheers, Mikkel.

2010/8/15 Phil Steitz <phil.steitz@gmail.com>:
> Jbb wrote:
>> Thanks Phil (and Luc).
>>
>> I understand the point, but is there any known workaround?
>> Is there anyway to get more precision on the root finder solver taken into account
that calculation time (speed of the algorithm) is not important in this case?
>> Anyway to exchange speed for accuracy?
>
> The only thing that I can think of is to recast the function so that
> its values near the root are distinguishable from zero.  If the
> function is returning values indistinguishable from zero near the
> root, there is nothing that the solver can do to distinguish the one
> you are looking for.  I did try changing the function definition to
> eliminate the pow and sqrt, but that did not help much.  If you have
> access to R, you can see the problem by executing the following:
>
> x <- seq( -2*pi - .00005, -2*pi + .00005, length = 100)
> y <- x * x - ((1 - cos(x)) * (1 - cos(x))) - ((2*pi - sin(x)) *
> (2*pi - sin(x)))
> plot(x, y, axes=T, type="b", pch="*", xlab=" ", ylab=" ")
>
> The y vector is your function with sqrt and pow removed.  You can
> see from the plot and also by just examining the y values that the
> function is not too well-behaved near the root and appears constant
> at zero over (-2pi - .00001, -2pi + .00001).
>
> You can see the same thing by tabulating values using your Java code:
>
> AlfaFunctionClass alfaFunction = new AlfaFunctionClass(x, y, z);
> final double answer = -2d * Math.PI;
> final double diff = .00005;
> final int count = 100;
> final double delta = diff / count;
> double root = answer - diff / 2;
> for (int i = 0; i < count; i++) {
>  System.out.println(
>   "point: " + root +
>   " value: " + alfaFunction.value(root) +
>   " arg error: " + (root + 2*Math.PI) +
>   " close to zero: " + (Math.abs(alfaFunction.value(root)) <1E-300)
>         );
>   root += delta;
> }
>
> You will see that the above displays "true" for "close to zero"
> throughout the interval described above. Replacing "<1E-300" with
> "==0d" yields the same result.  So bottom line is that from the
> perspective of the JVM, the function vanishes identically over the
> an interval of length approximately 1E-5 around -2pi. So the only
> way you are going to be able to distinguish the "true" root
> numerically is to transform the function so that values near the
> root are distinguishable from zero.
>
> Phil
>
>
>>
>> Thanks again and best regards,
>>
>> Jbb
>>
>> El 15/08/2010, a las 05:47, Phil Steitz <phil.steitz@gmail.com> escribió:
>>
>>> Juan Barandiaran wrote:
>>>> Hi, I'm working in a theoretical physics problem in which I have to find
the
>>>> roots of a function for every point in x,y,z.
>>>> I can bracket quite precisely where every root will be found and I know the
>>>> exact solution for many points:
>>>> For every point (x= 1, y= n*2*PI, z= 0) the root is in -n*2*PI.
>>>>
>>>> I have been using the newBrentSolver function, but to my surprise, even when
>>>> I set the Accuracy settings to very
>>>> small values and the MaxIterationCount to Integer.MAX_VALUE I get an error
>>>> of 1.427E-5 (too much).
>>>>
>>>> The function to be solved is alfa in
>>>> : alfa+Math.sqrt(Math.pow(x-Math.cos(alfa),2.0)+Math.pow(y-Math.sin(alfa),2.0)+z*z);
>>>>
>>>> Any hint of how could I get more precision on a root finder solver?
>>>> I attach the code of a short Test Case program that shows the error I'm
>>>> getting.
>>> The problem is in the function evaluation.  The function vanishes
>>> numerically (i.e., returns a value with absolute value less than
>>> 1E-90) throughout a neighborhood of width 1E-5 around the root you
>>> are trying to find.  The solver will (correctly, given its contract)
>>> return the first value that it finds within this interval.
>>>
>>> Phil
>>>> Thanks and best regards,
>>>>
>>>> JBB
>>>>
>>>> /*****************************************************************************************************************************************
>>>>
>>>>   ELECTRON DENSITY OF ENERGY: ALFA Retarded Position Calculation
>>>>   SpinningParticles Model
>>>>
>>>> @author Martin Rivas (SpinningParticles theory) & Juan Barandiaran
>>>> (Numerical Analysis & Calculation Program)
>>>> @web http://spinningparticles.com/
>>>>
>>>> Test case to try different solvers from the Apache Commons-Math library,
>>>> showing 1.427E-5 (too much) error in the resulting roots.
>>>>
>>>> Compiled in JAVA with Eclipse SDK 3.6:
>>>> http://www.eclipse.org/downloads/packages/eclipse-classic-360/heliosr
>>>> and with the Apache Commons-Math library 2.1  <
>>>> http://commons.apache.org/math/>
>>>>
>>>> ******************************************************************************************************************************************/
>>>>
>>>>
>>>> package testCaseRootFindingSolversPackage;
>>>>
>>>> import org.apache.commons.math.analysis.UnivariateRealFunction;
>>>> import org.apache.commons.math.analysis.solvers.*;
>>>> import org.apache.commons.math.ConvergenceException;
>>>> import org.apache.commons.math.FunctionEvaluationException;
>>>>
>>>> public class TestCaseRootFindingSolvers  {
>>>>
>>>> public static class AlfaFunctionClass implements UnivariateRealFunction {
>>>> //Alfa is the angle of the point dependent retarded position where the
>>>> electron charge was when it produced the field that influences the evaluated
>>>> x,y,z point
>>>> //To calculate it we must find the root of its equation between
>>>> -sqr((x^2+y^2+z^2))-1 and -sqr((x^2+y^2+z^2))+1 as the correct root must
be
>>>> negative (retarded)
>>>> //and bracketed around the module of the position vector
>>>>
>>>> private double x;
>>>> private double y;
>>>> private double z;
>>>>
>>>> public AlfaFunctionClass(double x, double y, double z) {
>>>> this.x  = x;
>>>> this.y = y;
>>>> this.z = z;
>>>> }
>>>> // this is the method specified by interface UnivariateRealFunction
>>>> public double value(double alfa) {
>>>> // here, we have to evaluate the function that calculates the retarded
>>>> position angle alfa.
>>>> return
>>>> alfa+Math.sqrt(Math.pow(x-Math.cos(alfa),2.0)+Math.pow(y-Math.sin(alfa),2.0)+z*z);
>>>> }
>>>> }
>>>>
>>>> public static class EnergyDensityFunctionAlfaClass {
>>>> public double evaluate(double x, double y, double z) {
>>>> double alfa;
>>>> UnivariateRealFunction alfaFunction = new AlfaFunctionClass(x,y,z);
>>>> UnivariateRealSolverFactory factory =
>>>> UnivariateRealSolverFactory.newInstance();
>>>> UnivariateRealSolver solver = factory.newBrentSolver();  //Fails when the
>>>> result in the two bracketing points has the same sign
>>>> // UnivariateRealSolver solver = factory.newBisectionSolver();
>>>>
>>>> //Trying to get the maximal precision with the available settings
>>>> solver.setAbsoluteAccuracy(1E-90);
>>>> solver.setMaximalIterationCount(Integer.MAX_VALUE);
>>>> solver.setRelativeAccuracy(1E-90);
>>>> solver.setFunctionValueAccuracy(1E-90);
>>>> try {
>>>> alfa = solver.solve(alfaFunction, -1.0*Math.sqrt(x*x+y*y+z*z)-1.0,
>>>> -1.0*Math.sqrt(x*x+y*y+z*z)+1.0);
>>>> } catch (ConvergenceException e) {
>>>> // TODO Auto-generated catch block
>>>> e.printStackTrace();
>>>> alfa= 0.0;
>>>> } catch (FunctionEvaluationException e) {
>>>> // TODO Auto-generated catch block
>>>> e.printStackTrace();
>>>> alfa= 0.0;
>>>> }
>>>> System.out.println("getFunctionValue="+String.valueOf(solver.getFunctionValue())+"???
>>>> Can't be cero with 1,42E-5 error?????");
>>>> System.out.println("getResult="+String.valueOf(solver.getResult()));
>>>> System.out.println("getAbsoluteAccuracy="+String.valueOf(solver.getAbsoluteAccuracy()));
>>>>
>>>> System.out.println("getRelativeAccuracy="+String.valueOf(solver.getRelativeAccuracy()));
>>>> System.out.println("getFunctionValueAccuracy="+String.valueOf(solver.getFunctionValueAccuracy()));
>>>> return alfa;
>>>> }
>>>> }
>>>> public static void main(String args[]) {
>>>> // Variables
>>>> double alfa=0.0;
>>>> double error;
>>>> EnergyDensityFunctionAlfaClass energyDensityFunction= new
>>>> EnergyDensityFunctionAlfaClass();
>>>> //Example point, although we have to evaluate the function for every point
>>>> in space.
>>>> double x=1.0;
>>>> double n=1.0;
>>>> double y=n*2.0*Math.PI;
>>>> double z=0.0;
>>>> alfa= energyDensityFunction.evaluate(x, y, z);
>>>> System.out.println();
>>>> System.out.println("Alfa= Root calculated with Apache Solver=
>>>> "+String.valueOf(alfa));
>>>> error= alfa+n*2*Math.PI;
>>>> System.out.println();
>>>> System.out.println("Error= Difference between the Alfa Root calculated with
>>>> Apache Solver and n*2*PI (exact root for any integer n)=
>>>> "+String.valueOf(error));
>>>> } // End of main
>>>> } //End
>>>>
>>>
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: user-help@commons.apache.org
>>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>> For additional commands, e-mail: user-help@commons.apache.org
>>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
> For additional commands, e-mail: user-help@commons.apache.org
>
>

---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message