#
# = Ordinary Differential Equations
# This chapter describes functions for solving ordinary differential equation
# (ODE) initial value problems. The library provides a variety of low-level
# methods, such as Runge-Kutta and Bulirsch-Stoer routines, and higher-level
# components for adaptive step-size control. The components can be combined
# by the user to achieve the desired solution, with full access to any
# intermediate steps.
#
#
# Contents:
# 1. {Classes for ODE solver}[link:rdoc/odeiv_rdoc.html#label-Classes+for+ODE+solver]
# 1. {Class Descriptions}[link:rdoc/odeiv_rdoc.html#label-Class+Descriptions]
#    1. {GSL::Odeiv::System : Defining the ODE System}[link:rdoc/odeiv_rdoc.html#label-System]
#    1. {GSL::Odeiv::Step : Stepping Algorithms}[link:rdoc/odeiv_rdoc.html#label-Step]
#    1. {GSL::Odeiv::Control : Adaptive Step-size Control}[link:rdoc/odeiv_rdoc.html#label-Control]
#    1. {GSL::Odeiv::Evolve : Evolution}[link:rdoc/odeiv_rdoc.html#label-Evolve]
#    1. {GSL::Odeiv::Solver : Higher level interface}[link:rdoc/odeiv_rdoc.html#label-Solver]
# 1. {Examples}[link:rdoc/odeiv_rdoc.html#label-Example]
#
# == 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 higher-level 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 <tt>dim</tt>.
#
#       # 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 <tt>t</tt> to <tt>t+h</tt> for a fixed step-size <tt>h</tt> and estimate the resulting local error.
#
# ---
# * GSL::Odeiv::Step.alloc(T, dim)
#
#   Constructor for a stepping function of an algorithm type <tt>T</tt> for a system of
#   dimension <tt>dim</tt>. The algorithms are specified by one of the constants under the
#   <tt>GSL::Odeiv::Step</tt> class, as
#
#   1. <tt>GSL::Odeiv::Step::RK2</tt>, Embedded 2nd order Runge-Kutta with 3rd order error estimate.
#   1. <tt>GSL::Odeiv::Step::RK4</tt>, 4th order (classical) Runge-Kutta.
#   1. <tt>GSL::Odeiv::Step::RKF45</tt>, Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error estimate. This method is a good general-purpose integrator.
#   1. <tt>GSL::Odeiv::Step::RKCK</tt>, Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error estimate.
#   1. <tt>GSL::Odeiv::Step::RK8PD</tt>, Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimate.
#   1. <tt>GSL::Odeiv::Step::RK2IMP</tt>, Implicit 2nd order Runge-Kutta at Gaussian points
#   1. <tt>GSL::Odeiv::Step::RK4IMP</tt>, Implicit 4th order Runge-Kutta at Gaussian points
#   1. <tt>GSL::Odeiv::Step::BSIMP</tt>, Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.
#   1. <tt>GSL::Odeiv::Step::GEAR1</tt>, M=1 implicit Gear method
#   1. <tt>GSL::Odeiv::Step::GEAR2</tt>, M=2 implicit Gear method
#   1. <tt>GSL::Odeiv::Step::RK2SIMP</tt> (GSL-1.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,
#
#   1. "<tt>rk2</tt>" or "<tt>gsl_odeiv_step_rk2</tt>"
#   1. "<tt>rk4</tt>" or "<tt>gsl_odeiv_step_rk4</tt>"
#   1. "<tt>rkf45</tt>" or "<tt>gsl_odeiv_step_rkf45</tt>"
#   1. "<tt>rkck</tt>" or "<tt>gsl_odeiv_step_rkck</tt>"
#   1. "<tt>rk8pd</tt>" or "<tt>gsl_odeiv_step_rk8pd</tt>"
#   1. "<tt>rk2imp</tt>" or "<tt>gsl_odeiv_step_rk2imp</tt>"
#   1. "<tt>rk4imp</tt>" or "<tt>gsl_odeiv_step_rk4imp</tt>"
#   1. "<tt>bsimp</tt>" or "<tt>gsl_odeiv_step_bsimp</tt>"
#   1. "<tt>gear1</tt>" or "<tt>gsl_odeiv_step_gear1</tt>"
#   1. "<tt>gear2</tt>" or "<tt>gsl_odeiv_step_gear2</tt>"
#   1. "<tt>rk2simp</tt>" or "<tt>gsl_odeiv_step_rk2simp</tt>" (GSL-1.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
#   <tt>dydt</tt>, using the step size <tt>h</tt> to advance the system from time
#   <tt>t</tt> and state <tt>y</tt> to time <tt>t+h</tt>. The new state of the system
#   is stored in <tt>y</tt> on output, with an estimate of the absolute error in
#   each component stored in <tt>yerr</tt>. If the argument <tt>dydt_in</tt> is not
#   <tt>nil</tt> it should be a {GSL::Vector}[link:rdoc/vector_rdoc.html] object
#   containing the derivatives for the system at time <tt>t</tt> 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 <tt>t+h</tt> will be
#   stored in <tt>dydt_out</tt> 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 <tt>epsabs</tt> and <tt>epsrel</tt>, and
#   scaling factors <tt>a_y</tt> and <tt>a_dydt</tt> for the system state
#   <tt>y(t)</tt> and derivatives <tt>y'(t)</tt> 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 <tt>epsabs</tt> and relative error
#   of <tt>epsrel</tt> with respect to the solution <tt>y_i(t)</tt>.
#   This is equivalent to the standard control object with <tt>a_y=1</tt>
#   and <tt>a_dydt=0</tt>.
#
# ---
# * 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 <tt>epsabs</tt> and
#   relative error of <tt>epsrel</tt> with respect to the derivatives of the
#   solution <tt>y'_i(t)</tt>.
#   This is equivalent to the standard control object with <tt>a_y=0</tt>
#   and <tt>a_dydt=1</tt>.
#
# ---
# * GSL::Odeiv::Control.alloc(epsabs, epsrel, a_y, a_dydt, vscale, dim)
#
#   This creates a new control object which uses the same algorithm as
#   <tt>GSL::Odeiv::Control.standard_new</tt> but with an absolute error which
#   is scaled for each component by the <tt>GSL::Vector</tt> object <tt>vscale</tt>.
#
# ---
# * GSL::Odeiv::Control#init(epsabs, epsrel, a_y, a_dydt)
#
#   This method initializes the controler with the parameters <tt>epsabs</tt>
#   (absolute error),  <tt>epsrel</tt> (relative error), <tt>a_y</tt>
#   (scaling factor for y) and <tt>a_dydt</tt> (scaling factor for derivatives).
#
# ---
# * GSL::Odeiv::Control#name
# * GSL::Odeiv::Control#hadjust(step, y0, yerr, dydt, h)
#
#   This method adjusts the step-size <tt>h</tt> using the control function
#   object, and the current values of <tt>y</tt>,  <tt>yerr</tt> and <tt>dydt</tt>.
#   The stepping function <tt>step</tt> is also needed to determine the order
#   of the method. On output, an array of two elements [<tt>hadj, status</tt>]
#   is returned: If the error in the y-values  <tt>yerr</tt> is found to be
#   too large then the step-size <tt>h</tt> is reduced and the method returns
#   [<tt>hadj, status</tt>=<tt>GSL::ODEIV::HADJ_DEC</tt>].
#   If the error is sufficiently small then
#   <tt>h</tt> may be increased and [<tt>hadj, status</tt>=<tt>GSL::ODEIV::HADJ_INC</tt>]
#   is returned.
#   The method returns [<tt>hadj, status</tt>=<tt>GSL::ODEIV::HADJ_NIL</tt>] 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.
#
# === Evolve
# The higher level of the system is the <tt>GSL::Evolve</tt> class which combines the
# results of a stepper and controler to reliably advance the solution forward
# over an interval <tt>(t_0, t_1)</tt>. If the controler signals that the step-size
# should be decreased the <tt>GSL::Evolve</tt> 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.alloc(dim)
#
#   These create a <tt>GSL::Evolve</tt> object for a system of <tt>dim</tt> 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 <tt>sys</tt> from time <tt>t</tt> and position
#   <tt>y</tt> using the stepping function <tt>step</tt>. The initial step-size is
#   taken as <tt>h</tt>. The maximum time <tt>t1</tt> is guaranteed not to be exceeded by
#   the time-step. On output, an array of three elements is returned,
#   [<tt>tnext, hnext, status</tt>], where <tt>tnext</tt> is the time advanced,
#   <tt>hnext</tt> is the step-size
#   for the next step, and <tt>status</tt> is an error code retunred by <tt>gsl_odeiv_evolve_apply()</tt> function.
#   On the final time-step the value of <tt>tnext</tt> will be set to
#   <tt>t1</tt> 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 <tt>T</tt> for the system of dimention <tt>dim</tt>. Here <tt>cary</tt> is an array as an argument for the <tt>GSL::Odeive:Control</tt> constructor.
#
#   * Ex1
#       solver = Odeiv::Solver.alloc(Odeiv::Step::RKF45, [1e-6, 0.0], func, dim)
#
#     * Type: RKF45,
#     * Control: epsabs = 1e-6, epsrel = 0.0, a_y = 1, a_dydt = 0
#     * System: function = <tt>func</tt>, jacobian = <tt>nil</tt>
#     * Dimension: dim
#
#   * Ex2:
#       solver = Odeiv::Solver.alloc(Odeiv::Step::BSIMP, [1e-6, 0.0, 1, 0], func, jac, dim)
#
#     * Type: BSIMP,
#     * Control: epsabs = 1e-6, epsrel = 0.0, a_y = 1, a_dydt = 0
#     * System: function = <tt>func</tt>, jacobian = <tt>jac</tt>
#     * 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 <tt>t</tt> and position <tt>y</tt> (<tt>GSL::Vector</tt> object) using the stepping function. On output, the new time and position are returned as an array [<tt>tnext, hnext, status</tt>], i.e. <tt>t, y</tt> themselves are not modified by this method. The maximum time <tt>t1</tt> is guaranteed not to be exceeded by the time-step. On the final time-step the value of <tt>tnext</tt> will be set to <tt>t1</tt> exactly.
#
# == Example
#
# 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),
#
# * x' = y
# * 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.alloc(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.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
#
#
# {prev}[link:rdoc/siman_rdoc.html]
# {next}[link:rdoc/interp_rdoc.html]
#
# {Reference index}[link:rdoc/ref_rdoc.html]
# {top}[link:index.html]
#
#