This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
ODE problem
- From: Anze Slosar <anze at mrao dot cam dot ac dot uk>
- To: <gsl-discuss at sources dot redhat dot com>
- Date: Wed, 23 Apr 2003 23:43:20 +0100 (BST)
- Subject: 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]