[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #40092] false position root finding requires too many fun
From: |
Patrick Alken |
Subject: |
[Bug-gsl] [bug #40092] false position root finding requires too many function evals |
Date: |
Mon, 23 Sep 2013 15:49:25 +0000 |
User-agent: |
Mozilla/5.0 (X11; Linux i686 on x86_64; rv:23.0) Gecko/20100101 Firefox/23.0 |
URL:
<http://savannah.gnu.org/bugs/?40092>
Summary: false position root finding requires too many
function evals
Project: GNU Scientific Library
Submitted by: psa
Submitted on: Mon 23 Sep 2013 03:49:24 PM GMT
Category: Performance
Severity: 3 - Normal
Operating System:
Status: None
Assigned to: None
Open/Closed: Open
Release:
Discussion Lock: Any
_______________________________________________________
Details:
>From address@hidden:
Hi,
I have noticed that the GSL implementation of the false position method
for one-dimensional root finding may require twice more function
evaluations than the bisection method for an almost linear function that
should not be considered as a degenerate case.
Since the reference states:
> Its convergence is linear, but it is usually faster than bisection.
I think it worths a remark in the reference.
I suppose that it is a design flaw that user supplied tolerance is not
available in the iteration functions. The implementation of the Brent
method uses DBL_EPSILON, but in general user might request significantly
lower precision. As a result it deteriorates programs performance.
The problem appears if at a certain step the edge of a bracketing
interval almost coincides the root. Further linear interpolation gives
the same point (either due to rounding error or withing the required
precision). Function evaluation should be skipped for this point before
bisection step.
The following program demonstrates the issue. The bisection method
converges after ~50 function evaluations, the false position method
requires ~100 function evaluations for the same precision. I would
stress again that it is a well behaviored almost linear function.
Maxim Nikulin
/* Excessive function evaluations for an almost linear function
* in the current implementation of the regula falsi method
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
double fpol2 (double x, void *params)
{
double *p = params;
return (x*p[2] + p[1])*x + p[0];
}
struct f_count_wrap_par {
gsl_function *F;
int n;
};
double f_count_wrap (double x, void *params)
{
struct f_count_wrap_par *f = params;
++f->n;
return GSL_FN_EVAL(f->F, x);
}
int main (void)
{
int i, max_iter = 50;
double x_i, x_l = 0.5, x_h = 1.1;
double tol_abs = 5*GSL_DBL_EPSILON;
double tol_rel = 5*GSL_DBL_EPSILON;
gsl_root_fsolver *s;
double par[] = { -1, 1, -1e-10 }; /* f(x) = -1e-10*x*x + (x - 1) */
gsl_function F = {fpol2, par}, Fwrap;
struct f_count_wrap_par wrap_par = {&F, 0};
int status;
double r = -2*par[0]/(sqrt(par[1]*par[1]-4*par[0]*par[2])+par[1]);
Fwrap.function = f_count_wrap;
Fwrap.params = &wrap_par;
s = gsl_root_fsolver_alloc (gsl_root_fsolver_falsepos);
/* s = gsl_root_fsolver_alloc (gsl_root_fsolver_bisection); */
gsl_root_fsolver_set (s, &Fwrap, x_l, x_h);
printf ("# i N_f x_l-r x_i-r x_h-r x_h-x_l\n");
for (i = 0, status = GSL_CONTINUE ;
status == GSL_CONTINUE && i < max_iter ; ++i) {
status = gsl_root_fsolver_iterate (s);
x_i = gsl_root_fsolver_root (s);
x_l = gsl_root_fsolver_x_lower (s);
x_h = gsl_root_fsolver_x_upper (s);
if (status != GSL_SUCCESS)
break;
status = gsl_root_test_interval (x_l, x_h, tol_abs, tol_rel);
printf ("%3d %5d % .3e % .3e % .3e % .3e\n", i, wrap_par.n,
x_l - r, x_i - r, x_h - r, x_h - x_l);
}
if (status == GSL_SUCCESS) {
printf ("# Converged\n");
printf ("# f(x_l) = % .4e, f(x_h) = % .4e\n",
GSL_FN_EVAL(&F, x_l), GSL_FN_EVAL(&F, x_h));
} else if (i == max_iter)
printf ("# Max iteration number exceeded\n");
gsl_root_fsolver_free (s);
exit (EXIT_SUCCESS);
}
Output: iteration, number of function evaluations, shifts from the root
of the left interval edge, the current root estimation, and the right
edge, and inally the bracketing interval length
# i N_f x_l-r x_i-r x_h-r x_h-x_l
0 4 -2.000e-01 5.000e-12 5.000e-12 2.000e-01
1 6 -1.000e-01 0.000e+00 0.000e+00 1.000e-01
2 8 -5.000e-02 0.000e+00 0.000e+00 5.000e-02
3 10 -2.500e-02 0.000e+00 0.000e+00 2.500e-02
4 12 -1.250e-02 0.000e+00 0.000e+00 1.250e-02
5 14 -6.250e-03 0.000e+00 0.000e+00 6.250e-03
6 16 -3.125e-03 0.000e+00 0.000e+00 3.125e-03
7 18 -1.563e-03 0.000e+00 0.000e+00 1.563e-03
8 20 -7.813e-04 0.000e+00 0.000e+00 7.813e-04
9 22 -3.906e-04 0.000e+00 0.000e+00 3.906e-04
10 24 -1.953e-04 0.000e+00 0.000e+00 1.953e-04
11 26 -9.766e-05 0.000e+00 0.000e+00 9.766e-05
12 28 -4.883e-05 0.000e+00 0.000e+00 4.883e-05
13 30 -2.441e-05 0.000e+00 0.000e+00 2.441e-05
14 32 -1.221e-05 0.000e+00 0.000e+00 1.221e-05
15 34 -6.104e-06 0.000e+00 0.000e+00 6.104e-06
16 36 -3.052e-06 0.000e+00 0.000e+00 3.052e-06
17 38 -1.526e-06 0.000e+00 0.000e+00 1.526e-06
18 40 -7.629e-07 0.000e+00 0.000e+00 7.629e-07
19 42 -3.815e-07 0.000e+00 0.000e+00 3.815e-07
20 44 -1.907e-07 0.000e+00 0.000e+00 1.907e-07
21 46 -9.537e-08 0.000e+00 0.000e+00 9.537e-08
22 48 -4.768e-08 0.000e+00 0.000e+00 4.768e-08
23 50 -2.384e-08 0.000e+00 0.000e+00 2.384e-08
24 52 -1.192e-08 0.000e+00 0.000e+00 1.192e-08
25 54 -5.960e-09 0.000e+00 0.000e+00 5.960e-09
26 56 -2.980e-09 0.000e+00 0.000e+00 2.980e-09
27 58 -1.490e-09 0.000e+00 0.000e+00 1.490e-09
28 60 -7.451e-10 0.000e+00 0.000e+00 7.451e-10
29 62 -3.725e-10 0.000e+00 0.000e+00 3.725e-10
30 64 -1.863e-10 0.000e+00 0.000e+00 1.863e-10
31 66 -9.313e-11 0.000e+00 0.000e+00 9.313e-11
32 68 -4.657e-11 0.000e+00 0.000e+00 4.657e-11
33 70 -2.328e-11 0.000e+00 0.000e+00 2.328e-11
34 72 -1.164e-11 0.000e+00 0.000e+00 1.164e-11
35 74 -5.821e-12 0.000e+00 0.000e+00 5.821e-12
36 76 -2.910e-12 0.000e+00 0.000e+00 2.910e-12
37 78 -1.455e-12 0.000e+00 0.000e+00 1.455e-12
38 80 -7.276e-13 0.000e+00 0.000e+00 7.276e-13
39 82 -3.637e-13 0.000e+00 0.000e+00 3.637e-13
40 84 -1.819e-13 0.000e+00 0.000e+00 1.819e-13
41 86 -9.104e-14 0.000e+00 0.000e+00 9.104e-14
42 88 -4.552e-14 0.000e+00 0.000e+00 4.552e-14
43 90 -2.265e-14 0.000e+00 0.000e+00 2.265e-14
44 92 -1.132e-14 0.000e+00 0.000e+00 1.132e-14
45 94 -5.773e-15 0.000e+00 0.000e+00 5.773e-15
46 96 -2.887e-15 0.000e+00 0.000e+00 2.887e-15
47 98 -1.332e-15 0.000e+00 0.000e+00 1.332e-15
# Converged
# f(x_l) = -1.3240e-15, f(x_h) = 8.2399e-18
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?40092>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [bug #40092] false position root finding requires too many function evals,
Patrick Alken <=