GSL::Odeiv::System
GSL::Odeiv::Step
GSL::Odeiv::Control
GSL::Odeiv::Evolve
GSL::Odeiv::Solver
The lowest level components are the stepping functions which advance a solution from time t to t+h for a fixed step-size h and estimate the resulting local error.
GSL::Odeiv::Step.new(T, dim)
Constructor for a stepping function of an algorithm type T for a system of dimension dim. The algorithms are specified by one of the constants under the GSL::Odeiv::Step
class, as
GSL::Odeiv::Step::RK2
, Embedded 2nd order Runge-Kutta with 3rd order error estimate.GSL::Odeiv::Step::RK4
, 4th order (classical) Runge-Kutta.GSL::Odeiv::Step::RKF45
, Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error estimate. This method is a good general-purpose integrator.GSL::Odeiv::Step::RKCK
, Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error estimate.GSL::Odeiv::Step::RK8PD
, Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimate.GSL::Odeiv::Step::RK2IMP
, Implicit 2nd order Runge-Kutta at Gaussian pointsGSL::Odeiv::Step::RK4IMP
, Implicit 4th order Runge-Kutta at Gaussian pointsGSL::Odeiv::Step::BSIMP
, Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.GSL::Odeiv::Step::GEAR1
, M=1 implicit Gear methodGSL::Odeiv::Step::GEAR2
, M=2 implicit Gear methodGSL::Odeiv::Step::RK2SIMP
(GSL-1.6)Ex:
step = Odeiv::Step.new(Odeiv::Step::RKF45, 2)
The algorithm types can also be given by a String, same as the C struct name,
rk2
" or "gsl_odeiv_step_rk2
"rk4
" or "gsl_odeiv_step_rk4
"rkf45
" or "gsl_odeiv_step_rkf45
"rkck
" or "gsl_odeiv_step_rkck
"rk8pd
" or "gsl_odeiv_step_rk8pd
"rk2imp
" or "gsl_odeiv_step_rk2imp
"rk4imp
" or "gsl_odeiv_step_rk4imp
"bsimp
" or "gsl_odeiv_step_bsimp
"gear1
" or "gsl_odeiv_step_gear1
"gear2
" or "gsl_odeiv_step_gear2
"rk2simp
" or "gsl_odeiv_step_rk2simp
" (GSL-1.6)Ex:
step = Odeiv::Step.new("bsimp", 4) step2 = Odeiv::Step.new("gsl_odeiv_step_rkck", 3)
GSL::Odeiv::Step#reset
GSL::Odeiv::Step#name
Returns the name of the stepper as a String. For example,
require("gsl") include Odeiv s = Step.new(Step::RK4, 2) printf("step method is '%s'\n", s.name)
would print something like step method is 'rk4'.
GSL::Odeiv::Step#order
GSL::Odeiv::Step#apply(t, h, y, yerr, dydt_in, dydt_out, sys)
GSL::Odeiv::Step#apply(t, h, y, yerr, dydt_in, sys)
GSL::Odeiv::Step#apply(t, h, y, yerr, sys)
nil
it should be a GSL::Vector object containing the derivatives for the system at time t on input. This is optional as the derivatives will be computed internally if they are not provided, but allows the reuse of existing derivative information. On output the new derivatives of the system at time t+h will be stored in dydt_out if it is not nil.GSL::Odeiv::System.new(func, jac, dim)
GSL::Odeiv::System.new(func, dim)
Constructor. This defines a general ODE system with the dimension dim.
Ex:
# t: variable (scalar) # y: vector # dydt: vector # params: scalar or an array func = Proc.new { |t, y, dydt, params| mu = params dydt[0] = y[1] dydt[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1.0) }
The third argument jac is also a Proc objects, to calculate jacobian.
Ex:
# t: scalar # y: vector # dfdy: matrix, jacobian # dfdt: vector # params: scalar of an array jac = Proc.new { |t, y, dfdy, dfdt, params| mu = params dfdy.set(0, 0, 0.0) dfdy.set(0, 1, 1.0) dfdy.set(1, 0, -2*mu*y[0]*y[1] - 1.0) dfdy.set(1, 1, -mu*(y[0]*y[0] - 1.0)) dfdt[0] = 0.0 dfdt[1] = 0.0 }
Note that some of the simpler solver algorithms do not make use of the Jacobian matrix, so it is not always strictly necessary to provide it. Thus the constructor is called as
sys = GSL:Odeiv::System.new(func, jac, dim) # for "BSIMP" algorithm sys = GSL:Odeiv::System.new(func, nil, dim) # for others, replaced by nil sys = GSL:Odeiv::System.new(func, dim) # or omit
GSL::Odeiv::System#set(func, jac, parameters...)
GSL::Odeiv::System#set_params(...)
GSL::Odeiv::System#function
GSL::Odeiv::System#func
GSL::Odeiv::System#jacobian
GSL::Odeiv::System#jac
GSL::Odeiv::System#dimension
GSL::Odeiv::System#dim
GSL::Odeiv::Control.standard_new(epsabs, epsrel, a_y, a_dydt)
GSL::Odeiv::Control.new(epsabs, epsrel, a_y, a_dydt)
GSL::Odeiv::Control.y_new(epsabs, epsrel)
GSL::Odeiv::Control.yp_new(epsabs, epsrel)
GSL::Odeiv::Control.new(epsabs, epsrel, a_y, a_dydt, vscale, dim)
GSL::Odeiv::Control.standard_new
but with an absolute error which is scaled for each component by the GSL::Vector
object vscale.GSL::Odeiv::Control#init(epsabs, epsrel, a_y, a_dydt)
GSL::Odeiv::Control#name
GSL::Odeiv::Control#hadjust(step, y0, yerr, dydt, h)
GSL::ODEIV::HADJ_DEC
]. If the error is sufficiently small then h may be increased and [hadj, status=GSL::ODEIV::HADJ_INC
] is returned. The method returns [hadj, status=GSL::ODEIV::HADJ_NIL
] if the step-size is unchanged. The goal of the method is to estimate the largest step-size which satisfies the user-specified accuracy requirements for the current point.The higher level of the system is the GSL::Evolve
class which combines the results of a stepper and controler to reliably advance the solution forward over an interval (t_0, t_1)
. If the controler signals that the step-size should be decreased the GSL::Evolve
object backs out of the current step and tries the proposed smaller step-size. This process is continued until an acceptable step-size is found.
GSL::Odeiv::Evolve.new(dim)
GSL::Odeiv::Evolve.alloc(dim)
GSL::Evolve
object for a system of dim dimensions.GSL::Odeiv::Evolve#reset
GSL::Odeiv::Evolve#apply(evolve, control, step, sys, t, t1, h, y)
gsl_odeiv_evolve_apply()
function. On the final time-step the value of tnext will be set to t1 exactly.GSL::Odeiv::Evolve#count
This is the highest level interface to solve ODE system, which contains System, Step, Control, and Evolve classes.
GSL::Odeiv::Solver.new(T, cary, fac, jac, dim)
GSL::Odeiv::Solver.new(T, cary, fac, dim)
GSL::Odeive:Control
constructor.
Ex1
solver = Odeiv::Solver.new(Odeiv::Step::RKF45, [1e-6, 0.0], func, dim)
Ex2:
solver = Odeiv::Solver.new(Odeiv::Step::BSIMP, [1e-6, 0.0, 1, 0], func, jac, dim)
GSL::Odeiv:::Solver#reset
GSL::Odeiv:::Solver#step
GSL::Odeiv:::Solver#control
GSL::Odeiv:::Solver#evolve
GSL::Odeiv:::Solver#system
GSL::Odeiv::System#set_params(...)
GSL::Odeiv:::Solver#apply(t, t1, h, y)
GSL::Vector
object) using the stepping function. On output, the new time and position are returned as an array [tnext, hnext, status], i.e. t, y themselves are not modified by this method. The maximum time t1 is guaranteed not to be exceeded by the time-step. On the final time-step the value of tnext will be set to t1 exactly.The following program solves the second-order nonlinear Van der Pol oscillator equation, as found in the GSL manual, x"(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0,
This can be converted into a first order system suitable for use with the routines described in this chapter by introducing a separate variable for the velocity, y = x'(t),
y' = -x + \mu y (1-x^2)
require("gsl") include Odeiv dim = 2 # dimension of the system # Proc object to calculate the derivatives func = Proc.new { |t, y, dydt, mu| dydt[0] = y[1] dydt[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1.0) } # Create the solver solver = Solver.new(Step::RKF45, [1e-6, 0.0], func, dim) mu = 10.0 solver.set_params(mu) t = 0.0 # initial time t1 = 100.0 # end time h = 1e-6 # initial step y = Vector.new([1.0, 0.0]) # initial value while t < t1 t, h, status = solver.apply(t, t1, h, y) break if status != GSL::SUCCESS printf("%.5e %.5e %.5e %.5e\n", t, y[0], y[1], h) end