This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: multidimensional root-finding problem :-(
- From: Ivo Alxneit <ivo dot alxneit at psi dot ch>
- To: Slaven Peles <peles at cns dot physics dot gatech dot edu>
- Cc: gsl-discuss at sources dot redhat dot com
- Date: Fri, 23 May 2003 15:53:24 +0200
- Subject: Re: multidimensional root-finding problem :-(
- Organization: PSI / ENE / HTS
- References: <200305222234.16084.peles@cns.physics.gatech.edu>
-----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-----