This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: ODEs - runge kutta


On Thursday 21 February 2002 11:19 am, Sam Halliday wrote:
> hello,
>
> i am currently trying to numerically solve the lane-emden equation, which
> is a 2nd order differential equation, and experience has shown that the 4th
> order runge kutta is the best way to solve this.
>
> i was interrested in using the GSL for solving this problem, although the
> examples section for ODE's scared me too much... must i really define so
> much, (jacobians and so forth), is there not a way in which i can just send
> my function, parameter (which actually is the time) and time step length to
> a library and let it do all the work? does anyone have a simple example i
> can follow so that i can use this library for this task? teh one in the
> documenattion seems to be too crypic, and tries to do too much for my
> needs...
>

It looks scary, indeed. But, at this point it seems to me that results 
obtained with gsl are much more accurate than those obtained using Numerical 
Recipes, for example. So I guess it is worth to put a little effort in 
learning how to use gsl. Below is a working example for solving Van der Pol 
equation using RK4 fixed stepsize, without defining Jacobian. This is the 
best what I could do at the moment, so if anybody has some better example, 
please send it to the list. 

Cheers,
Slaven



#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

int
func (double t, const double y[], double f[],
      void *params)
{
  double mu = *(double *)params;
  f[0] = y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, 
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat 
    = gsl_matrix_view_array (dfdy=NULL, 2, 2);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

int
main (void)
{
  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rk4;

  gsl_odeiv_step * s 
    = gsl_odeiv_step_alloc (T, 2);

  double mu = 10;
  gsl_odeiv_system sys = {func, jac, 2, &mu};

  double t = 0.0, t1 = 100.0;
  double h = 1e-2;
  double y[2] = { 1.0, 0.0 }, y_err[2], dfdy[4], dydt_in[2], dydt_out[2];

  /* initialise dydt_in */
  GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in);

  while (t < t1)
    {
      int status = gsl_odeiv_step_apply (s, t, h, y, y_err, 
                                         dydt_in, dydt_out, &sys);

      if (status != GSL_SUCCESS)
          break;

      dydt_in[0] = dydt_out[0];
      dydt_in[1] = dydt_out[1];

      t += h;

      printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

  gsl_odeiv_step_free(s);
  return 0;
}


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]