Last Update: Oct. 4, 2007
/***********************************************/
/* T.Kouya's GSL sample program collection */
/* Solving Nonlinear System of Equations */
/* Written by Tomonori Kouya */
/* */
/* Version 0.01: 2007-10-04 */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multiroots.h>
/* Dimension of Matrix and Vectors */
#define DIM 3
/* Limit of Iterative Times */
#define MAXTIMES 100
/* func(x) = 0 */
/* x1 + x2 + x3 - 6 = 0 */
/* x1 * x2 * x3 - 6 = 0 */
/* x1^2 + x2^2 + x3^2 - 14 = 0 */
int func_m(const gsl_vector *x, void *param, gsl_vector *ret)
{
double x1, x2, x3;
x1 = gsl_vector_get(x, 0);
x2 = gsl_vector_get(x, 1);
x3 = gsl_vector_get(x, 2);
gsl_vector_set(ret, 0, x1 + x2 + x3 - 6.0);
gsl_vector_set(ret, 1, x1 * x2 * x3 - 6.0);
gsl_vector_set(ret, 2, x1*x1 + x2*x2 + x3*x3 - 14.0);
return GSL_SUCCESS;
}
/* J = [ 1, 1, 1 ] */
/* [x2 * x3, x1 * x3, x1 * x2 ] */
/* [2 * x1 , 2 * x2, 2 * x3 ] */
int dfunc_m(const gsl_vector *x, void *param, gsl_matrix *ret)
{
double x1, x2, x3;
x1 = gsl_vector_get(x, 0);
x2 = gsl_vector_get(x, 1);
x3 = gsl_vector_get(x, 2);
gsl_matrix_set(ret, 0, 0, 1.0);
gsl_matrix_set(ret, 0, 1, 1.0);
gsl_matrix_set(ret, 0, 2, 1.0);
gsl_matrix_set(ret, 1, 0, x2 * x3);
gsl_matrix_set(ret, 1, 1, x1 * x3);
gsl_matrix_set(ret, 1, 2, x1 * x2);
gsl_matrix_set(ret, 2, 0, 2.0 * x1);
gsl_matrix_set(ret, 2, 1, 2.0 * x2);
gsl_matrix_set(ret, 2, 2, 2.0 * x3);
return GSL_SUCCESS;
}
/* simultaneous eval of f and df */
int fdfunc_m(const gsl_vector *x, void *param, gsl_vector *ret_vec, gsl_matrix *ret_mat)
{
func_m(x, param, ret_vec);
dfunc_m(x, param, ret_mat);
return GSL_SUCCESS;
}
int main(void)
{
int i, times, status;
gsl_multiroot_function f;
gsl_multiroot_fsolver *workspace_f;
gsl_multiroot_function_fdf fdf;
gsl_multiroot_fdfsolver *workspace;
gsl_vector *x;
/* allocate a, x, b */
x = gsl_vector_alloc(DIM);
/* Fsolver */
/* Define Solver */
workspace_f = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrid, DIM);
printf("F solver: %s\n", gsl_multiroot_fsolver_name(workspace_f));
f.f = &func_m;
f.n = DIM;
f.params = 0;
/* set initial value */
for(i = 0; i < DIM; i++)
gsl_vector_set(x, i, 1.0);
/* set solver */
gsl_multiroot_fsolver_set(workspace_f, &f, x);
/* main loop */
for(times = 0; times < MAXTIMES; times++)
{
status = gsl_multiroot_fsolver_iterate(workspace_f);
printf("%d times: ", times);
for(i = 0; i < DIM; i++)
printf("%10.3e ", gsl_vector_get(workspace_f->x, i));
printf("\n");
if((status == GSL_EBADFUNC) || (status == GSL_ENOPROG))
{
printf("Status: %s\n", gsl_strerror(status));
break;
}
}
/* print answer */
for(i = 0; i < DIM; i++)
printf("%3d %25.17e\n", i, gsl_vector_get(workspace_f->x, i));
gsl_multiroot_fsolver_free(workspace_f);
/* FDFsolver */
/* Define Solver */
workspace = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybridsj, DIM);
printf("FDF solver: %s\n", gsl_multiroot_fdfsolver_name(workspace));
fdf.f = &func_m;
fdf.df = &dfunc_m;
fdf.fdf = &fdfunc_m;
fdf.n = DIM;
fdf.params = 0;
/* set initial value */
for(i = 0; i < DIM; i++)
gsl_vector_set(x, i, 1.0);
/* set solver */
gsl_multiroot_fdfsolver_set(workspace, &fdf, x);
/* main loop */
for(times = 0; times < MAXTIMES; times++)
{
status = gsl_multiroot_fdfsolver_iterate(workspace);
printf("%d times: ", times);
for(i = 0; i < DIM; i++)
printf("%10.3e ", gsl_vector_get(workspace->x, i));
printf("\n");
if((status == GSL_EBADFUNC) || (status == GSL_ENOPROG))
{
printf("Status: %s\n", gsl_strerror(status));
break;
}
}
/* print answer */
for(i = 0; i < DIM; i++)
printf("%3d %25.17e\n", i, gsl_vector_get(workspace->x, i));
gsl_multiroot_fdfsolver_free(workspace);
/* free x */
gsl_vector_free(x);
return 0;
}
$ ./nonlinear_system F solver: hybrid 0 times: 1.000e+00 1.000e+00 1.000e+00 1 times: 1.000e+00 1.000e+00 1.000e+00 2 times: 1.000e+00 1.000e+00 1.000e+00 3 times: 1.000e+00 1.000e+00 1.000e+00 4 times: 2.006e+00 1.994e+00 2.000e+00 5 times: 2.006e+00 1.994e+00 2.000e+00 6 times: 2.081e+00 3.180e+00 7.391e-01 7 times: 2.056e+00 2.785e+00 1.158e+00 8 times: 2.066e+00 2.944e+00 9.899e-01 9 times: 2.068e+00 2.965e+00 9.669e-01 10 times: 2.068e+00 2.963e+00 9.691e-01 11 times: 2.068e+00 2.963e+00 9.691e-01 12 times: 2.068e+00 2.963e+00 9.691e-01 13 times: 1.999e+00 3.002e+00 9.986e-01 14 times: 2.000e+00 3.000e+00 9.999e-01 15 times: 2.000e+00 3.000e+00 1.000e-00 16 times: 2.000e+00 3.000e+00 1.000e-00 17 times: 2.000e+00 3.000e+00 1.000e-00 18 times: 2.000e+00 3.000e+00 1.000e+00 19 times: 2.000e+00 3.000e+00 1.000e-00 20 times: 2.000e+00 3.000e+00 1.000e-00 21 times: 2.000e+00 3.000e+00 1.000e-00 22 times: 2.000e+00 3.000e+00 1.000e-00 23 times: 2.000e+00 3.000e+00 1.000e-00 24 times: 2.000e+00 3.000e+00 1.000e-00 25 times: 2.000e+00 3.000e+00 1.000e-00 26 times: 2.000e+00 3.000e+00 1.000e-00 27 times: 2.000e+00 3.000e+00 1.000e-00 28 times: 2.000e+00 3.000e+00 1.000e-00 29 times: 2.000e+00 3.000e+00 1.000e-00 Status: iteration is not making progress towards solution 0 2.00000000000000133e+00 1 2.99999999999999911e+00 2 9.99999999999999778e-01 FDF solver: hybridsj 0 times: 1.000e+00 1.000e+00 1.000e+00 1 times: 1.000e+00 1.000e+00 1.000e+00 2 times: 1.000e+00 1.000e+00 1.000e+00 3 times: 1.000e+00 1.000e+00 1.000e+00 4 times: 2.006e+00 1.994e+00 2.000e+00 5 times: 2.006e+00 1.994e+00 2.000e+00 6 times: 2.081e+00 3.180e+00 7.391e-01 7 times: 2.056e+00 2.785e+00 1.158e+00 8 times: 2.066e+00 2.944e+00 9.899e-01 9 times: 2.068e+00 2.965e+00 9.669e-01 10 times: 2.068e+00 2.963e+00 9.691e-01 11 times: 2.068e+00 2.963e+00 9.691e-01 12 times: 2.068e+00 2.963e+00 9.691e-01 13 times: 2.028e+00 2.987e+00 9.884e-01 14 times: 2.000e+00 3.001e+00 9.994e-01 15 times: 2.000e+00 3.000e+00 1.000e-00 16 times: 2.000e+00 3.000e+00 1.000e-00 17 times: 2.000e+00 3.000e+00 1.000e-00 18 times: 2.000e+00 3.000e+00 1.000e-00 19 times: 2.000e+00 3.000e+00 1.000e-00 20 times: 2.000e+00 3.000e+00 1.000e+00 21 times: 2.000e+00 3.000e+00 1.000e+00 22 times: 2.000e+00 3.000e+00 1.000e+00 23 times: 2.000e+00 3.000e+00 1.000e+00 24 times: 2.000e+00 3.000e+00 1.000e+00 25 times: 2.000e+00 3.000e+00 1.000e+00 26 times: 2.000e+00 3.000e+00 1.000e+00 27 times: 2.000e+00 3.000e+00 1.000e+00 28 times: 2.000e+00 3.000e+00 1.000e+00 29 times: 2.000e+00 3.000e+00 1.000e+00 30 times: 2.000e+00 3.000e+00 1.000e+00 Status: iteration is not making progress towards solution 0 1.99999999999999556e+00 1 3.00000000000000222e+00 2 1.00000000000000155e+00