Hello, and thanks in advance :)
I'm currently working on solving an ODE, which I already know what it the solution supposed to "look like" . I wanted to try with GSL ODE solvers, but I'm having some trouble, because I'm getting something completely wrong!
I just wanted to see if the problem is in the GSL code, since I'm only "guessing" on how to use the functions- (I couldn't find more examples or so to see for myself) or if it is something else.
I'm guessing the problem is in the code, maybe I'm not using the functions properly? I think so, because I can't even get the right answer for a harmonic oscillator!
I searched for old post of GSL problems but the other post were too old, or not related at all.
/************************************************/
Ok, the equation is this:
x' = y
y' = 2/ x^3
and the code is something like this:
//(the derivates function)
func (double t, const double y[], double f[],void *params)
{
double mu = *(double *)params;
//ex=exp(-mu*t);
f[0] = y[1];
f[1] = (1/(y[0]*y[0]*y[0]))*2.0;
return GSL_SUCCESS;
}
//the jacobian
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, 2, 2);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, 0.0);
gsl_matrix_set (m, 0, 1, 1.0);
gsl_matrix_set (m, 1, 0, -6/(y[0]*y[0]*y[0]*y[0]));
gsl_matrix_set (m, 1, 1, 0.0);
dfdt[0] = 0.0;
dfdt[1] = -mu*K*t/(y[0]*y[0]*y[0]);
return GSL_SUCCESS;
}
//main (partial) content
const gsl_odeiv_step_type * T = gsl_odeiv_step_gear1;
gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
gsl_odeiv_control * c= gsl_odeiv_control_y_new (1e-3,1e-3);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);
double mu = 0.0;
gsl_odeiv_system sys = {func, jac, 2, &mu};
double t = 0.0, t1 = 1e11;
double h = 1e-6;
double y[2] = { 1.0, -5.1e-12 };As I said, any help is appreciated! What do you think of this?
(And sorry for my bad English ^^u )
