Hi Lukasz,
I'm not sure if I understand what you mean by "I cannot explicitly put
an array in the ODE system definition". I guess that you want the
right hand side of the ODE to depend not only on the time t, but also
on an array of doubles? You can do that using the 'params' argument of
the function which is of type 'void *' and can point to anything, see
the attached example.
Regards,
Frank
2008/4/23, [EMAIL PROTECTED] <[EMAIL PROTECTED]>:
> Does anyone know how can I represent pulsed stimuli in ODEs system
> definition. I need the driving term in the equation be the function of time,
> and have the form of a rectangular pulse. The problem is, I cannot explicitly
> put an array in the ODE system definition. Is there any way to do it?
> Best,
> Lukasz
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
// The function calculating the rhs of the ODE
int
func (double t, const double y[], double dydt[],
void *params)
{
// params is casted to 'double *' and assigned to my_array.
// Now 'my_array' is identical to the variable 'array' in main ()
// as intended.
const double * my_array = (const double *) params;
// Evaluate the function
dydt [0] = t + my_array [0] + my_array [1] + my_array [2];
return GSL_SUCCESS;
}
int
main (void)
{
// This is the array that is used as a parameter of the function
// that calculates the rhs of the ODE
double array [3] = {2, 4, 6};
// Initialise ODE solver
const gsl_odeiv_step_type * T
= gsl_odeiv_step_rkf45;
gsl_odeiv_step * s
= gsl_odeiv_step_alloc (T, 1);
gsl_odeiv_control * c
= gsl_odeiv_control_y_new (1e-6, 0.0);
gsl_odeiv_evolve * e
= gsl_odeiv_evolve_alloc (1);
gsl_odeiv_system sys = {func, 0, 1, array};
double t = 0.0, t1 = 100.0;
double h = 1e-6;
double y [1] = { 0.0 };
while (t < t1)
{
int status = gsl_odeiv_evolve_apply (e, c, s,
&sys,
&t, t1,
&h, y);
if (status != GSL_SUCCESS)
break;
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
return 0;
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
