htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/ode-initval.c

   1 #include <stdio.h>
   2 #include <gsl/gsl_errno.h>
   3 #include <gsl/gsl_matrix.h>
   4 #include <gsl/gsl_odeiv.h>
   5
   6 int
   7 func (double t, const double y[], double f[],
   8       void *params)
   9 {
  10   double mu = *(double *)params;
  11   f[0] = y[1];
  12   f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  13   return GSL_SUCCESS;
  14 }
  15
  16 int
  17 jac (double t, const double y[], double *dfdy,
  18      double dfdt[], void *params)
  19 {
  20   double mu = *(double *)params;
  21   gsl_matrix_view dfdy_mat
  22     = gsl_matrix_view_array (dfdy, 2, 2);
  23   gsl_matrix * m = &dfdy_mat.matrix;
  24   gsl_matrix_set (m, 0, 0, 0.0);
  25   gsl_matrix_set (m, 0, 1, 1.0);
  26   gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  27   gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  28   dfdt[0] = 0.0;
  29   dfdt[1] = 0.0;
  30   return GSL_SUCCESS;
  31 }
  32
  33 int
  34 main (void)
  35 {
  36   const gsl_odeiv_step_type * T
  37     = gsl_odeiv_step_rk8pd;
  38
  39   gsl_odeiv_step * s
  40     = gsl_odeiv_step_alloc (T, 2);
  41   gsl_odeiv_control * c
  42     = gsl_odeiv_control_y_new (1e-6, 0.0);
  43   gsl_odeiv_evolve * e
  44     = gsl_odeiv_evolve_alloc (2);
  45
  46   double mu = 10;
  47   gsl_odeiv_system sys = {func, jac, 2, &mu};
  48
  49   double t = 0.0, t1 = 100.0;
  50   double h = 1e-6;
  51   double y[2] = { 1.0, 0.0 };
  52
  53   while (t < t1)
  54     {
  55       int status = gsl_odeiv_evolve_apply (e, c, s,
  56                                            &sys,
  57                                            &t, t1,
  58                                            &h, y);
  59
  60       if (status != GSL_SUCCESS)
  61           break;
  62
  63       printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
  64     }
  65
  66   gsl_odeiv_evolve_free (e);
  67   gsl_odeiv_control_free (c);
  68   gsl_odeiv_step_free (s);
  69   return 0;
  70 }