# GSL sample: Nonlinear System of Equations

Last Update: Oct. 4, 2007

### Sample code

nonlinear_system.c

```/***********************************************/
/* 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;
}
}

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;
}
}

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;
}```

### Example execution

```\$ ./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```

<- Go back