[Help-gsl] ode problem

[Tomasz Samotyjak, PWSZ]

Hi,

I've tried to use gsl fo ode's. The jac and func is:

int
func (double t, const double y[], double f[], void *params)
{
  double* mu = (double *)params;
  double R1 = mu[0], R2 = R1;
  double L1 = mu[1], L2 = L1;
  double C1 = mu[2];

  f[0] = 1/L1 * ( 400*sin(314*t) - f[1] - R1 * f[0] );
  f[1] = 1/C1 * ( f[0] - f[2] );
  f[2] = 1/L2 * ( f[1] - R2 * f[2] );

  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
  double* mu = (double *)params;
  double R1 = mu[0], R2 = R1;
  double L1 = mu[1], L2 = L1;
  double C1 = mu[2];
  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 3,3);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, -R1/L1);
  gsl_matrix_set (m, 0, 1, -1/L1);
  gsl_matrix_set (m, 0, 2, 0);

  gsl_matrix_set (m, 1, 0, 1/C1);
  gsl_matrix_set (m, 1, 1, 0);
  gsl_matrix_set (m, 1, 2, -1/C1);

  gsl_matrix_set (m, 2, 0, 0);
  gsl_matrix_set (m, 2, 1, 1/L2);
  gsl_matrix_set (m, 2, 2, -R2/L2);

  for( int i=0 ; i<3 ; i++)  dfdt[i] = 0.0;
  return GSL_SUCCESS;
}

The main function is:

int
main (void)
{
  const gsl_odeiv_step_type * T = gsl_odeiv_step_bsimp;
  gsl_odeiv_step * s  = gsl_odeiv_step_alloc (T, 3);
  gsl_odeiv_control * c  = gsl_odeiv_control_yp_new (1e-3,0);
  gsl_odeiv_evolve * e  = gsl_odeiv_evolve_alloc (3);

  double mu[] = {1, 1e-3, 1e-3};
  gsl_odeiv_system sys = {func, jac, 3, mu};

  double t = 0.0, t1 = 0.20;
  double h = 1e-5;
  double y[] = { 0.0 , 0.0, 0.0};

  ofstream out("wynik.txt");

  while (t < t1)
    {
     int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys,
                                           &t, t1,
                                           &h, y);

      if (status != GSL_SUCCESS)
          break;

      cout << t << "\t" << y[0] << "\t" << y[1] << "\t" << y[2] << endl;
      out <<  t << "\t" << y[0] << "\t" << y[1] << "\t" << y[2] << endl;
    }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  out.close();
  return 0;
}

But program don't work and prints:

8.192e-015 0 0 0
4.9152e-014 0 0 0
2.53952e-013 0 0 0
1.27795e-012 0 0 0
6.39795e-012 0 0 0
3.1998e-011 0 0 0
1.59998e-010 0 0 0
7.99998e-010 0 0 0
4e-009 -1.#IND -1.#IND -1.#IND
2e-008 -1.#IND -1.#IND -1.#IND
1e-007 -1.#IND -1.#IND -1.#IND
5e-007 -1.#IND -1.#IND -1.#IND
2.5e-006 -1.#IND -1.#IND -1.#IND
1.25e-005 -1.#IND -1.#IND -1.#IND
6.25e-005 -1.#IND -1.#IND -1.#IND
0.0003125 -1.#IND -1.#IND -1.#IND
0.0015625 -1.#IND -1.#IND -1.#IND
0.0078125 -1.#IND -1.#IND -1.#IND
0.0390625 -1.#IND -1.#IND -1.#IND
0.195312 -1.#IND -1.#IND -1.#IND
0.2 -1.#IND -1.#IND -1.#IND

I've tried change abs error and other parameters, but the result is the
same. What to do?

Best
Tomasz Samotyjak