/* * The program solves the second-order nonlinear * Van der Pol oscillator equation, * * u''(t) + \mu u'(t) (u(t)^2 - 1) + u(t) = 0 * * This can be converted into a first order system * by introducing a separate variable for the velocity, * v = u'(t), * * u' = v * v' = -u + \mu v (1-u^2) * * The program begins by defining functions for these * derivatives. The main function uses * driver level functions to solve the problem. The program * evolves the solution from (u, v) = (1, 0) at t=0 to t=100. * The step-size h is automatically adjusted by the controller * to maintain an absolute accuracy of 10^{-6} in the function * values (u, v). */ #include #include #include #include int func (double t, const double y[], double f[], void *params) { double mu = *(double *) params; f[0] = y[1]; f[1] = -y[0] + mu*y[1]*(1. - y[0]*y[0]); return GSL_SUCCESS; } int main (void) { size_t neqs = 2; /* number of equations */ double eps_abs = 1.e-6, eps_rel = 0.; /* desired precision */ double stepsize = 1e-6; /* initial integration step size */ double mu = 0.1; /* the parameter of the oscillator */ double t = 0., t1 = 100.; /* time interval of the evolution */ int status; /* * Initial conditions */ double y[2] = { 1., 0. }; /* * 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, &mu}; /* * Evolution loop */ while (t < t1) { status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &stepsize, y); if (status != GSL_SUCCESS) { printf ("Troubles at % .5e % .5e % .5e\n", t, y[0], y[1]); break; } printf ("% .5e % .5e % .5e\n", t, y[0], y[1]); } gsl_odeiv2_evolve_free (e); gsl_odeiv2_control_free (c); gsl_odeiv2_step_free (s); return 0; }