# # = 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] # #