Mailing-List: contact user-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Users List" <user@commons.apache.org>
Received-SPF: pass (nike.apache.org: domain of phil.steitz@gmail.com
 designates 209.85.212.43 as permitted sender)
DomainKey-Signature: a=rsa-sha1; c=nofws;
        d=gmail.com; s=gamma;
        h=message-id:date:from:user-agent:mime-version:to:subject:references
         :in-reply-to:x-enigmail-version:content-type
         :content-transfer-encoding;
        b=dDWlMZ8JLAavcXHEP1nqVFLoIx6x9MK6ZGYxaS9ew9FqGQzDbLszIF16Ll9jtdeWzi
         esq5azeO9PpVA19ifP4mdz2gbpFXh/3lTjvc0BxDixiYgCLpOqlO3eLYoZTkHXnWs8mP
         wMw7TuAEFJ3dJvu5d3u6MpiS+HooAoUh2q6OQ=
Message-ID: <4C67635A.9080008@gmail.com>
Date: Sat, 14 Aug 2010 23:47:38 -0400
From: Phil Steitz <phil.steitz@gmail.com>
User-Agent: Thunderbird 2.0.0.24 (Macintosh/20100228)
MIME-Version: 1.0
To: Commons Users List <user@commons.apache.org>
Subject: Re: [math] Accuracy in root finding with newBrentSolver
References: <AANLkTimoToDZShWz=BF5n6MieKg7oM49LC0vgs9Z5tQO@mail.gmail.com>
In-Reply-To: <AANLkTimoToDZShWz=BF5n6MieKg7oM49LC0vgs9Z5tQO@mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit

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