Last Update: Oct. 4, 2007
/***********************************************/
/* T.Kouya's GSL sample program collection */
/* Solving Nonlinear Equation */
/* Written by Tomonori Kouya */
/* */
/* Version 0.01: 2007-10-04 */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
/* Limit of Iterative Times */
#define MAXTIMES 100
/* func(x) = 0 */
/* x^2 - 2 = 0 */
double func(double x, void *param)
{
return x * x - 2;
}
/* f'(x) */
double dfunc(double x, void *param)
{
return 2.0 * x;
}
/* simultaneous eval of f and df */
void fdfunc(double x, void *param, double *ret_f, double *ret_df)
{
*ret_f = func(x, param);
*ret_df = dfunc(x, param);
return;
}
int main(void)
{
int i, times, status;
gsl_function f;
gsl_root_fsolver *workspace_f;
gsl_function_fdf fdf;
gsl_root_fdfsolver *workspace;
double x, x_l, x_r;
/* Fsolver */
/* Define Solver */
workspace_f = gsl_root_fsolver_alloc(gsl_root_fsolver_bisection);
printf("F solver: %s\n", gsl_root_fsolver_name(workspace_f));
f.function = &func;
f.params = 0;
/* set initial interval */
x_l = 0.0;
x_r = 2 * M_PI;
/* set solver */
gsl_root_fsolver_set(workspace_f, &f, x_l, x_r);
/* main loop */
for(times = 0; times < MAXTIMES; times++)
{
status = gsl_root_fsolver_iterate(workspace_f);
x_l = gsl_root_fsolver_x_lower(workspace_f);
x_r = gsl_root_fsolver_x_upper(workspace_f);
printf("%d times: [%10.3e, %10.3e]\n", times, x_l, x_r);
status = gsl_root_test_interval(x_l, x_r, 1.0e-13, 1.0e-20);
if(status != GSL_CONTINUE)
{
printf("Status: %s\n", gsl_strerror(status));
printf("\n Root = [%25.17e, %25.17e]\n\n", x_l, x_r);
break;
}
}
/* free */
gsl_root_fsolver_free(workspace_f);
/* FDFsolver */
/* Define Solver */
workspace = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
printf("FDF solver: %s\n", gsl_root_fdfsolver_name(workspace));
fdf.f = &func;
fdf.df = &dfunc;
fdf.fdf = &fdfunc;
fdf.params = 0;
/* set initial value */
x = 2.0;
/* set solver */
gsl_root_fdfsolver_set(workspace, &fdf, x);
/* main loop */
for(times = 0; times < MAXTIMES; times++)
{
status = gsl_root_fdfsolver_iterate(workspace);
x = gsl_root_fdfsolver_root(workspace);
status = gsl_root_test_residual(func(x, NULL), 1.0e-5);
printf("%d times: %10.3e\n", times, x);
if((status != GSL_CONTINUE) || (status == GSL_EZERODIV) || (status == GSL_EZERODIV))
{
printf("Status: %s\n", gsl_strerror(status));
printf("\n Root = %25.17e\n\n", x);
break;
}
}
/* free */
gsl_root_fdfsolver_free(workspace);
return 0;
}
$ ./nonlinear F solver: bisection 0 times: [ 0.000e+00, 3.142e+00] 1 times: [ 0.000e+00, 1.571e+00] 2 times: [ 7.854e-01, 1.571e+00] 3 times: [ 1.178e+00, 1.571e+00] 4 times: [ 1.374e+00, 1.571e+00] 5 times: [ 1.374e+00, 1.473e+00] 6 times: [ 1.374e+00, 1.424e+00] 7 times: [ 1.399e+00, 1.424e+00] 8 times: [ 1.411e+00, 1.424e+00] 9 times: [ 1.411e+00, 1.417e+00] 10 times: [ 1.411e+00, 1.414e+00] 11 times: [ 1.413e+00, 1.414e+00] 12 times: [ 1.414e+00, 1.414e+00] 13 times: [ 1.414e+00, 1.414e+00] 14 times: [ 1.414e+00, 1.414e+00] 15 times: [ 1.414e+00, 1.414e+00] 16 times: [ 1.414e+00, 1.414e+00] 17 times: [ 1.414e+00, 1.414e+00] 18 times: [ 1.414e+00, 1.414e+00] 19 times: [ 1.414e+00, 1.414e+00] 20 times: [ 1.414e+00, 1.414e+00] 21 times: [ 1.414e+00, 1.414e+00] 22 times: [ 1.414e+00, 1.414e+00] 23 times: [ 1.414e+00, 1.414e+00] 24 times: [ 1.414e+00, 1.414e+00] 25 times: [ 1.414e+00, 1.414e+00] 26 times: [ 1.414e+00, 1.414e+00] 27 times: [ 1.414e+00, 1.414e+00] 28 times: [ 1.414e+00, 1.414e+00] 29 times: [ 1.414e+00, 1.414e+00] 30 times: [ 1.414e+00, 1.414e+00] 31 times: [ 1.414e+00, 1.414e+00] 32 times: [ 1.414e+00, 1.414e+00] 33 times: [ 1.414e+00, 1.414e+00] 34 times: [ 1.414e+00, 1.414e+00] 35 times: [ 1.414e+00, 1.414e+00] 36 times: [ 1.414e+00, 1.414e+00] 37 times: [ 1.414e+00, 1.414e+00] 38 times: [ 1.414e+00, 1.414e+00] 39 times: [ 1.414e+00, 1.414e+00] 40 times: [ 1.414e+00, 1.414e+00] 41 times: [ 1.414e+00, 1.414e+00] 42 times: [ 1.414e+00, 1.414e+00] 43 times: [ 1.414e+00, 1.414e+00] 44 times: [ 1.414e+00, 1.414e+00] 45 times: [ 1.414e+00, 1.414e+00] Status: success Root = [ 1.41421356237308871e+00, 1.41421356237317797e+00] FDF solver: newton 0 times: 1.500e+00 1 times: 1.417e+00 2 times: 1.414e+00 Status: success Root = 1.41421568627450989e+00