/* * The program solves the following system of four * first order differential equations, which describe * the motion of a glider * * v' = -sin(theta) - v^2 / R * theta' = - cos(theta)/v + v * x' = v cos(theta) * y' = v sin(theta) * * Here v is the dimensionless speed of the glider, * theta is the angle that the velocity direction * makes with the horizontal, x and y are * dimensionless cartesian coordinates of the glider. * * The step-size of the integrator is automatically * adjusted by the controller to maintain the * requested accuracy */ #include #include #include #include int func (double t, const double y[], double f[], void *params); int main (void) { size_t neqs = 4; /* number of equations */ double eps_abs = 1.e-8, eps_rel = 0.; /* desired precision */ double stepsize = 1e-6; /* initial integration step */ double R = 5.; /* the aerodynamic efficiency */ double t = 0., t1 = 120.; /* time interval */ int status; /* * Initial conditions */ double y[4] = { 1.3, -0.1, 0., 2. }; /* for res1 */ //double y[4] = { 2.3, -0.1, 0., 2. }; /* for res2 */ //double y[4] = { 3.3, -0.1, 0., 2. }; /* for res3 */ /* * Explicit embedded Runge-Kutta-Fehlberg (4,5) method. * This method is a good general-purpose integrator. */ gsl_odeiv2_step *s = gsl_odeiv2_step_alloc (gsl_odeiv2_step_rkf45, neqs); gsl_odeiv2_control *c = gsl_odeiv2_control_y_new (eps_abs, eps_rel); gsl_odeiv2_evolve *e = gsl_odeiv2_evolve_alloc (neqs); gsl_odeiv2_system sys = {func, NULL, neqs, &R}; /* * Evolution loop */ while ( (t < t1) && (y[3] > 0) ) { status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &stepsize, y); if (status != GSL_SUCCESS) { printf ("Troubles: % .5e % .5e % .5e % .5e % .5e\n", t, y[0], y[1], y[2], y[3]); break; } printf ("% .5e % .5e % .5e % .5e % .5e\n", t, y[0], y[1], y[2], y[3]); } gsl_odeiv2_evolve_free (e); gsl_odeiv2_control_free (c); gsl_odeiv2_step_free (s); return 0; }