commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <phil.ste...@gmail.com>
Subject Re: [math] Accuracy in root finding with newBrentSolver
Date Sun, 15 Aug 2010 12:47:29 GMT
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


Mime
View raw message