Ordinary Differential Equations¶ ↑
This chapter describes functions for solving ordinary differential equation (ODE) initial value problems. The library provides a variety of lowlevel methods, such as RungeKutta and BulirschStoer routines, and higherlevel components for adaptive stepsize control. The components can be combined by the user to achieve the desired solution, with full access to any intermediate steps.
Contents:
Classes for ODE solver¶ ↑

GSL::Odeiv::System

GSL::Odeiv::Step

GSL::Odeiv::Control

GSL::Odeiv::Evolve
These are GSL structure wrappers.

GSL::Odeiv::Solver
Another higherlevel interface to ODE system classes.
Class Descriptions¶ ↑
System¶ ↑

GSL::Odeiv::System.alloc(func, jac, dim)

GSL::Odeiv::System.alloc(func, dim)
Constructor. This defines a general ODE system with the dimension
dim
.# 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) } # 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 } sys = GSL:Odeiv::System.alloc(func, jac, dim) # for "BSIMP" algorithm
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 as
sys = GSL:Odeiv::System.alloc(func, nil, dim) # for others, replaced by nil sys = GSL:Odeiv::System.alloc(func, dim) # or omit

GSL::Odeiv::System#set(func, jac, parameters…)
This method sets or resets the procedures to evaluate the function and jacobian, and the constant parameters.

GSL::Odeiv::System#set_params(…)
Set the constant parameters of the function.

GSL::Odeiv::System#function

GSL::Odeiv::System#func

GSL::Odeiv::System#jacobian

GSL::Odeiv::System#jac
Return Proc objects

GSL::Odeiv::System#dimension

GSL::Odeiv::System#dim
Step¶ ↑
The lowest level components are the stepping functions which advance a
solution from time t
to t+h
for a fixed stepsize
h
and estimate the resulting local error.

GSL::Odeiv::Step.alloc(T, dim)
Constructor for a stepping function of an algorithm type
T
for a system of dimensiondim
. The algorithms are specified by one of the constants under theGSL::Odeiv::Step
class, as
GSL::Odeiv::Step::RK2
, Embedded 2nd order RungeKutta with 3rd order error estimate. 
GSL::Odeiv::Step::RK4
, 4th order (classical) RungeKutta. 
GSL::Odeiv::Step::RKF45
, Embedded 4th order RungeKuttaFehlberg method with 5th order error estimate. This method is a good generalpurpose integrator. 
GSL::Odeiv::Step::RKCK
, Embedded 4th order RungeKutta CashKarp method with 5th order error estimate. 
GSL::Odeiv::Step::RK8PD
, Embedded 8th order RungeKutta PrinceDormand method with 9th order error estimate. 
GSL::Odeiv::Step::RK2IMP
, Implicit 2nd order RungeKutta at Gaussian points 
GSL::Odeiv::Step::RK4IMP
, Implicit 4th order RungeKutta at Gaussian points 
GSL::Odeiv::Step::BSIMP
, Implicit BulirschStoer method of Bader and Deuflhard. This algorithm requires the Jacobian. 
GSL::Odeiv::Step::GEAR1
, M=1 implicit Gear method 
GSL::Odeiv::Step::GEAR2
, M=2 implicit Gear method 
GSL::Odeiv::Step::RK2SIMP
(GSL1.6)

Ex:
step = Odeiv::Step.alloc(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
” (GSL1.6)

Ex:
step = Odeiv::Step.alloc("bsimp", 4) step2 = Odeiv::Step.alloc("gsl_odeiv_step_rkck", 3)


GSL::Odeiv::Step#reset
This method resets the stepper. It should be used whenever the next use of s will not be a continuation of a previous step.

GSL::Odeiv::Step#name
Returns the name of the stepper as a String. For example,
require("gsl") include Odeiv s = Step.alloc(Step::RK4, 2) printf("step method is '%s'\n", s.name)
would print something like step method is 'rk4'.

GSL::Odeiv::Step#order
Returns the order of the stepper on the previous step. This order can vary if the stepper itself is adaptive.

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)
This method applies the stepper to the system of equations defined by
dydt
, using the step sizeh
to advance the system from timet
and statey
to timet+h
. The new state of the system is stored iny
on output, with an estimate of the absolute error in each component stored inyerr
. If the argumentdydt_in
is notnil
it should be a GSL::Vector object containing the derivatives for the system at timet
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 timet+h
will be stored indydt_out
if it is not nil.
Control¶ ↑

GSL::Odeiv::Control.standard_new(epsabs, epsrel, a_y, a_dydt)

GSL::Odeiv::Control.alloc(epsabs, epsrel, a_y, a_dydt)
The standard control object is a four parameter heuristic based on absolute and relative errors
epsabs
andepsrel
, and scaling factorsa_y
anda_dydt
for the system statey(t)
and derivativesy'(t)
respectively.

GSL::Odeiv::Control.y_new(epsabs, epsrel)
This method creates a new control object which will keep the local error on each step within an absolute error of
epsabs
and relative error ofepsrel
with respect to the solutiony_i(t)
. This is equivalent to the standard control object witha_y=1
anda_dydt=0
.

GSL::Odeiv::Control.yp_new(epsabs, epsrel)
This method creates a new control object which will keep the local error on each step within an absolute error of
epsabs
and relative error ofepsrel
with respect to the derivatives of the solutiony'_i(t)
. This is equivalent to the standard control object witha_y=0
anda_dydt=1
.

GSL::Odeiv::Control.alloc(epsabs, epsrel, a_y, a_dydt, vscale, dim)
This creates a new control object which uses the same algorithm as
GSL::Odeiv::Control.standard_new
but with an absolute error which is scaled for each component by theGSL::Vector
objectvscale
.

GSL::Odeiv::Control#init(epsabs, epsrel, a_y, a_dydt)
This method initializes the controler with the parameters
epsabs
(absolute error),epsrel
(relative error),a_y
(scaling factor for y) anda_dydt
(scaling factor for derivatives).

GSL::Odeiv::Control#name

GSL::Odeiv::Control#hadjust(step, y0, yerr, dydt, h)
This method adjusts the stepsize
h
using the control function object, and the current values ofy
,yerr
anddydt
. The stepping functionstep
is also needed to determine the order of the method. On output, an array of two elements [hadj, status
] is returned: If the error in the yvaluesyerr
is found to be too large then the stepsizeh
is reduced and the method returns [hadj, status
=GSL::ODEIV::HADJ_DEC
]. If the error is sufficiently small thenh
may be increased and [hadj, status
=GSL::ODEIV::HADJ_INC
] is returned. The method returns [hadj, status
=GSL::ODEIV::HADJ_NIL
] if the stepsize is unchanged. The goal of the method is to estimate the largest stepsize which satisfies the userspecified accuracy requirements for the current point.
Evolve¶ ↑
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 stepsize should be decreased the GSL::Evolve
object backs out of the current step and tries the proposed smaller
stepsize. This process is continued until an acceptable stepsize is
found.

GSL::Odeiv::Evolve.alloc(dim)
These create a
GSL::Evolve
object for a system ofdim
dimensions.

GSL::Odeiv::Evolve#reset
This method resets the GSL::Evolve object. It should be used whenever the next use of e will not be a continuation of a previous step.

GSL::Odeiv::Evolve#apply(evolve, control, step, sys, t, t1, h, y)
This method advances the system
sys
from timet
and positiony
using the stepping functionstep
. The initial stepsize is taken ash
. The maximum timet1
is guaranteed not to be exceeded by the timestep. On output, an array of three elements is returned, [tnext, hnext, status
], wheretnext
is the time advanced,hnext
is the stepsize for the next step, andstatus
is an error code retunred bygsl_odeiv_evolve_apply()
function. On the final timestep the value oftnext
will be set tot1
exactly.

GSL::Odeiv::Evolve#count
Solver¶ ↑
This is the highest level interface to solve ODE system, which contains System, Step, Control, and Evolve classes.

GSL::Odeiv::Solver.alloc(T, cary, fac, dim)
This creates a ODE solver with the algorithm type
T
for the system of dimentiondim
. Herecary
is an array as an argument for theGSL::Odeive:Control
constructor.
Ex1
solver = Odeiv::Solver.alloc(Odeiv::Step::RKF45, [1e6, 0.0], func, dim)

Type: RKF45,

Control: epsabs = 1e6, epsrel = 0.0, a_y = 1, a_dydt = 0

System: function =
func
, jacobian =nil

Dimension: dim


Ex2:
solver = Odeiv::Solver.alloc(Odeiv::Step::BSIMP, [1e6, 0.0, 1, 0], func, jac, dim)

Type: BSIMP,

Control: epsabs = 1e6, epsrel = 0.0, a_y = 1, a_dydt = 0

System: function =
func
, jacobian =jac

Dimension: dim



GSL::Odeiv:::Solver#reset
Reset the solver elements (step, evolve)

GSL::Odeiv:::Solver#step

GSL::Odeiv:::Solver#control

GSL::Odeiv:::Solver#evolve

GSL::Odeiv:::Solver#system
Access to the solver elements.

GSL::Odeiv::System#set_params(…)
Set the constant parameters of the function to solve.

GSL::Odeiv:::Solver#apply(t, t1, h, y)
This method advances the system from time
t
and positiony
(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 timet1
is guaranteed not to be exceeded by the timestep. On the final timestep the value oftnext
will be set tot1
exactly.
Example¶ ↑
The following program solves the secondorder 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),

x' = y

y' = x + mu y (1x^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.alloc(Step::RKF45, [1e6, 0.0], func, dim) mu = 10.0 solver.set_params(mu) t = 0.0 # initial time t1 = 100.0 # end time h = 1e6 # initial step y = Vector.alloc([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