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]

ODE problem


Hello everybody!

I have big trouble with gsl solver on various fronts. One of them is the
following; consider this short  piece of code (I tried to make it as short
as possible):




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


int PerturbationFunc (double t, const double y[], double f[],
		       void *pars)
{

  f[0]=y[1];
  f[1]=-1e6*y[0];

  f[2]=y[3];
  f[3]=-1e6*y[2];
  return GSL_SUCCESS;
}




void main()
{

  double tau=-0.1;
  const gsl_odeiv_step_type * T =gsl_odeiv_step_rkf45;


  gsl_odeiv_step * s
    = gsl_odeiv_step_alloc (T, 4);
  gsl_odeiv_control * c
    = gsl_odeiv_control_y_new (0,1.e-15);
  gsl_odeiv_evolve * e
    = gsl_odeiv_evolve_alloc (4);



  gsl_odeiv_system sys = {PerturbationFunc, NULL, 1, NULL};

  double h = 1e-6;
  double y [4];
  int evolve=1000;
  double k=1000;

  y[0]=1.0/sqrt(2*k)*sin(k*tau);       // u
  y[1]=k*1.0/sqrt(2*k)*cos(k*tau);     // u'
  y[2]=1.0/sqrt(2*k)*sin(k*tau);       // u
  y[3]=k*1.0/sqrt(2*k)*cos(k*tau);     // u'
  double p;



  while (evolve)
    {
      int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys,
                                           &tau, 1.e200,
                                           &h, y);

      if (status != GSL_SUCCESS)
	{
	  return 0;
        }

      printf ("%g \n",y[2]-y[0]);
      evolve--;
    }

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);

}



Basically y[0] and y[1] is a SHM kind of differential equation with y[1]
being derivative of y[0]. Exactly the same is true for y[2] and y[3] with
*the same initial conditions.* Modulus numerical errors the ODE solver
should evolve them exactly the same. Ie. it should print just zeros! It
doesn't however, after a while they separate and soon the are different by
a non-negigible amount! You say it doesn't matter, but when I want to
follow a more complicated ODE with real and imag being out of phase by
pi/2 I soon loose quadrature and it is really annoying...

Any ideas or comments?



anze


--- Space isn't remote at all. It's only an hour's drive away,
    if your car could go straight upwards. [Fred Hoyle]











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