GSL求根

来源：互联网发布：win7网络图标有个红叉编辑：程序博客网时间：2024/05/16 15:57

demo_fn.h

#ifndef ___DEMO_FN_HHH_HH_H_#define ___DEMO_FN_HHH_HH_H_struct quadratic_params  {    double a, b, c;  };double quadratic (double x, void *params);double quadratic_deriv (double x, void *params);void quadratic_fdf (double x, void *params,                     double *y, double *dy);#endif

demo_fn.c

#include "demo_fn.h"doublequadratic (double x, void *params){  struct quadratic_params *p     = (struct quadratic_params *) params;  double a = p->a;  double b = p->b;  double c = p->c;  return (a * x + b) * x + c;}doublequadratic_deriv (double x, void *params){  struct quadratic_params *p     = (struct quadratic_params *) params;  double a = p->a;  double b = p->b;  double c = p->c;  return 2.0 * a * x + b;}voidquadratic_fdf (double x, void *params,                double *y, double *dy){  struct quadratic_params *p     = (struct quadratic_params *) params;  double a = p->a;  double b = p->b;  double c = p->c;  *y = (a * x + b) * x + c;  *dy = 2.0 * a * x + b;}

main.c

#include <stdio.h>#include <gsl/gsl_errno.h>#include <gsl/gsl_math.h>#include <gsl/gsl_roots.h>#include "demo_fn.h"// #include "demo_fn.c"#pragma comment(lib, "libgsl_d.lib")#pragma comment(lib, "libgslcblas_d.lib")intmain (void){  int status;  int iter = 0, max_iter = 100;  const gsl_root_fsolver_type *T;  gsl_root_fsolver *s;  double r = 0, r_expected = sqrt (5.0);  double x_lo = 0.0, x_hi = 5.0;  gsl_function F;  struct quadratic_params params = {1.0, 0.0, -5.0};  F.function = &quadratic;  F.params = ¶ms;  T = gsl_root_fsolver_brent;  s = gsl_root_fsolver_alloc (T);  gsl_root_fsolver_set (s, &F, x_lo, x_hi);  printf ("using %s method\n",           gsl_root_fsolver_name (s));  printf ("%5s [%9s, %9s] %9s %10s %9s\n",          "iter", "lower", "upper", "root",           "err", "err(est)");  do    {      iter++;      status = gsl_root_fsolver_iterate (s);      r = gsl_root_fsolver_root (s);      x_lo = gsl_root_fsolver_x_lower (s);      x_hi = gsl_root_fsolver_x_upper (s);      status = gsl_root_test_interval (x_lo, x_hi,                                       0, 0.001);      if (status == GSL_SUCCESS)        printf ("Converged:\n");      printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",              iter, x_lo, x_hi,              r, r - r_expected,               x_hi - x_lo);    }  while (status == GSL_CONTINUE && iter < max_iter);  return status;}