This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: multidimensional root-finding problem :-(


-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On Friday 23 May 2003 04:34, Slaven Peles wrote:
> Hi,
>
> I was working with the example in GSL Reference Manual for multidimensional
> root finding using Newton's method:
> http://sources.redhat.com/gsl/ref/gsl-ref_33.html#SEC445
> I wanted to have both, equations and the Jacobi matrix, in the same
> function in order to optimize calculation, so I left functions rosenbrock_f
> and rosenbrock_df empty, and put everything into rosenbrock_fdf (see
> below). According to the Manual this should work, however when I run this
> program my initial value for x remains unchanged with every iteration. If
> anyone on the list has experience with GSL multidimensional root-finding
> functions please help. My code is below.
>
> Many thanks,
> Slaven
>
> /**************************************************************************
>******/ #include <stdlib.h>
> #include <stdio.h>
> #include <gsl/gsl_vector.h>
> #include <gsl/gsl_multiroots.h>
>
> struct rparams
>   {
>     double a;
>     double b;
>   };
>
> int
> rosenbrock_f (const gsl_vector * x, void *params,
>               gsl_vector * f)
> {
>   return GSL_SUCCESS;
> }
>
> int
> rosenbrock_df (const gsl_vector * x, void *params,
>                gsl_matrix * J)
> {
>   return GSL_SUCCESS;
> }
>
> int
> rosenbrock_fdf (const gsl_vector * x, void *params,
>                 gsl_vector * f, gsl_matrix * J)
> {
>   double a = ((struct rparams *) params)->a;
>   double b = ((struct rparams *) params)->b;
>
>   double x0 = gsl_vector_get (x, 0);
>   double x1 = gsl_vector_get (x, 1);
>
>   double y0 = a * (1 - x0);
>   double y1 = b * (x1 - x0 * x0);
>
>   gsl_vector_set (f, 0, y0);
>   gsl_vector_set (f, 1, y1);
>
>   gsl_matrix_set (J, 0, 0, -a);
>   gsl_matrix_set (J, 0, 1, 0);
>   gsl_matrix_set (J, 1, 0,  -2 * b  * x0);
>   gsl_matrix_set (J, 1, 1, b);
>
>   return GSL_SUCCESS;
> }
>
>
> int
> print_state (size_t iter, gsl_multiroot_fdfsolver * s)
> {
>   printf ("iter = %3u x = % .3f % .3f "
>           "f(x) = % .3e % .3e\n",
>           iter,
>           gsl_vector_get (s->x, 0),
>           gsl_vector_get (s->x, 1),
>           gsl_vector_get (s->f, 0),
>           gsl_vector_get (s->f, 1));
> }
>
>
> int
> main (void)
> {
>   const gsl_multiroot_fdfsolver_type *T;
>   gsl_multiroot_fdfsolver *s;
>
>   int status;
>   size_t iter = 0;
>
>   const size_t n = 2;
>   struct rparams p = {1.0, 10.0};
>   gsl_multiroot_function_fdf f = {&rosenbrock_f,
>                                   &rosenbrock_df,
>                                   &rosenbrock_fdf,
>                                   n, &p};
>
>   double x_init[2] = {-10.0, -5.0};
>   gsl_vector *x = gsl_vector_alloc (n);
>
>   gsl_vector_set (x, 0, x_init[0]);
>   gsl_vector_set (x, 1, x_init[1]);
>
>   T = gsl_multiroot_fdfsolver_gnewton;
>   s = gsl_multiroot_fdfsolver_alloc (T, n);
>   gsl_multiroot_fdfsolver_set (s, &f, x);
>
>   print_state (iter, s);
>
>   do
>     {
>       iter++;
>
>       status = gsl_multiroot_fdfsolver_iterate (s);
>
>       print_state (iter, s);
>
>       if (status)
>         break;
>
>       status = gsl_multiroot_test_residual (s->f, 1e-7);
>     }
>   while (status == GSL_CONTINUE && iter < 1000);
>
>   printf ("status = %s\n", gsl_strerror (status));
>
>   gsl_multiroot_fdfsolver_free (s);
>   gsl_vector_free (x);
>   return 0;
> }
i think there is a missconception here
you need all three (f, df, and fdf) perform what they are supposed to do. your 
iteration fails because your f is just a dummy function.
the "additional" function fdf is used in cases where both f and df are needed 
at the same time. you often can code it to be much faster that separate calls 
to f and df. so the way it is in the manual is actually correct. 
- -- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
Paul Scherrer Institute             fax: +41 56 310 2624
CH-5232 Villigen                   http://solar.web.psi.ch
Switzerland                   gnupg key: 0x515E30C7
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.0 (GNU/Linux)

iD8DBQE+zifeAd7CE1FeMMcRAiS6AJ9+ZPXznDOuAgmKQh1FCO4Jj7+NLACgrn17
FK6c5DE++DVjAA/BdKNFkbI=
=1Q/q
-----END PGP SIGNATURE-----


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]