ode-initval.texi

开放gsl矩阵运算

relative error of @var{eps_rel} with respect to the derivatives of the
solution @math{y'_i(t)} .  This is equivalent to the standard control
object with @var{a_y}=0 and @var{a_dydt}=1.
@end deftypefun

@deftypefun {gsl_odeiv_control *} gsl_odeiv_control_alloc (const gsl_odeiv_control_type * @var{T})
This function returns a pointer to a newly allocated instance of a
control function of type @var{T}.  This function is only needed for
defining new types of control functions.  For most purposes the standard
control functions described above should be sufficient. 
@end deftypefun

@deftypefun int gsl_odeiv_control_init (gsl_odeiv_control * @var{c}, double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt})
This function initializes the control function @var{c} with the
parameters @var{eps_abs} (absolute error), @var{eps_rel} (relative
error), @var{a_y} (scaling factor for y) and @var{a_dydt} (scaling
factor for derivatives).
@end deftypefun

@deftypefun void gsl_odeiv_control_free (gsl_odeiv_control * @var{c})
This function frees all the memory associated with the control function
@var{c}.
@end deftypefun

@deftypefun int gsl_odeiv_control_hadjust (gsl_odeiv_control * @var{c}, gsl_odeiv_step * @var{s}, const double y0[], const double yerr[], const double dydt[], double * @var{h})
This function adjusts the step-size @var{h} using the control function
@var{c}, and the current values of @var{y}, @var{yerr} and @var{dydt}.
The stepping function @var{step} is also needed to determine the order
of the method.  If the error in the y-values @var{yerr} is found to be
too large then the step-size @var{h} is reduced and the function returns
@code{GSL_ODEIV_HADJ_DEC}.  If the error is sufficiently small then
@var{h} may be increased and @code{GSL_ODEIV_HADJ_INC} is returned.  The
function returns @code{GSL_ODEIV_HADJ_NIL} if the step-size is
unchanged.  The goal of the function is to estimate the largest
step-size which satisfies the user-specified accuracy requirements for
the current point.
@end deftypefun

@deftypefun {const char *} gsl_odeiv_control_name (const gsl_odeiv_control * @var{c})
This function returns a pointer to the name of the control function.
For example,

@example
printf("control method is '%s'\n", 
       gsl_odeiv_control_name (c));
@end example

@noindent
would print something like @code{control method is 'standard'}
@end deftypefun


@node Evolution
@section Evolution

The highest level of the system is the evolution function which combines
the results of a stepping function and control function to reliably
advance the solution forward over an interval @math{(t_0, t_1)}.  If the
control function signals that the step-size should be decreased the
evolution function backs out of the current step and tries the proposed
smaller step-size.  This is process is continued until an acceptable
step-size is found.

@deftypefun {gsl_odeiv_evolve *} gsl_odeiv_evolve_alloc (size_t @var{dim})
This function returns a pointer to a newly allocated instance of an
evolution function for a system of @var{dim} dimensions.
@end deftypefun

@deftypefun int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * @var{e}, gsl_odeiv_control * @var{con}, gsl_odeiv_step * @var{step}, const gsl_odeiv_system * @var{dydt}, double * @var{t}, double @var{t1}, double * @var{h}, double y[])
This function advances the system (@var{e}, @var{dydt}) from time
@var{t} and position @var{y} using the stepping function @var{step}.
The new time and position are stored in @var{t} and @var{y} on output.
The initial step-size is taken as @var{h}, but this will be modified
using the control function @var{c} to achieve the appropriate error
bound if necessary.  The routine may make several calls to @var{step} in
order to determine the optimum step-size. If the step-size has been
changed the value of @var{h} will be modified on output.  The maximum
time @var{t1} is guaranteed not to be exceeded by the time-step.  On the
final time-step the value of @var{t} will be set to @var{t1} exactly.
@end deftypefun

@deftypefun int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * @var{e})
This function resets the evolution function @var{e}.  It should be used
whenever the next use of @var{e} will not be a continuation of a
previous step.
@end deftypefun

@deftypefun void gsl_odeiv_evolve_free (gsl_odeiv_evolve * @var{e})
This function frees all the memory associated with the evolution function
@var{e}.
@end deftypefun

@node ODE Example programs
@section Examples
@cindex Van der Pol oscillator, example
The following program solves the second-order nonlinear Van der Pol
oscillator equation,

@tex
\beforedisplay
$$
x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
$$
\afterdisplay
@end tex
@ifinfo
@example
x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
@end example
@end ifinfo
@noindent
This can be converted into a first order system suitable for use with
the library by introducing a separate variable for the velocity, @math{y
= x'(t)},

@tex
\beforedisplay
$$
\eqalign{
x' &= y\cr
y' &= -x + \mu y (1-x^2)
}
$$
\afterdisplay
@end tex
@ifinfo
@example
x' = y
y' = -x + \mu y (1-x^2)
@end example
@end ifinfo
@noindent
The program begins by defining functions for these derivatives and
their Jacobian,

@example
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

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]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
@}

int
jac (double t, const double y[], double *dfdy, 
     double dfdt[], void *params)
@{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat 
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix; 
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
@}

int
main (void)
@{
  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rk8pd;

  gsl_odeiv_step * s 
    = gsl_odeiv_step_alloc (T, 2);
  gsl_odeiv_control * c 
    = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e 
    = gsl_odeiv_evolve_alloc (2);

  double mu = 10;
  gsl_odeiv_system sys = @{func, jac, 2, &mu@};

  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y[2] = @{ 1.0, 0.0 @};

  gsl_ieee_env_setup();

  while (t < t1)
    @{
      int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys, 
                                           &t, t1,
                                           &h, y);

      if (status != GSL_SUCCESS)
          break;

      printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    @}

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);
  return 0;
@}
@end example
@noindent
The main loop of the program evolves the solution from @math{(y, y') =
(1, 0)} at @math{t=0} to @math{t=100}.  The step-size @math{h} is
automatically adjusted by the controller to maintain an absolute
accuracy of @c{$10^{-6}$} 
@math{10^@{-6@}} in the function values @var{y}.

@iftex
@sp 1
@center @image{vdp}
@center Numerical solution of the Van der Pol oscillator equation 
@center using Prince-Dormand 8th order Runge-Kutta.
@end iftex

@noindent
To obtain the values at regular intervals, rather than the variable
spacings chosen by the control function, the main loop can be modified
to advance the solution from one point to the next.  For example, the
following main loop prints the solution at the fixed points @math{t = 0,
1, 2, \dots, 100},

@example
  for (i = 1; i <= 100; i++)
    @{
      double ti = i * t1 / 100.0;

      while (t < ti)
        @{
          gsl_odeiv_evolve_apply (e, c, s, 
                                  &sys, 
                                  &t, ti, &h,
                                  y);
        @}
 
      printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    @}
@end example


@node ODE References and Further Reading
@section References and Further Reading
@noindent
Many of the the basic Runge-Kutta formulas can be found in the Handbook of
Mathematical Functions,

@itemize @asis
@item
Abramowitz & Stegun (eds.), @cite{Handbook of Mathematical Functions},
Section 25.5.
@end itemize
@noindent
The implicit Bulirsch-Stoer algorithm @code{bsimp} is described in the
following paper,

@itemize @asis
@item
G. Bader and P. Deuflhard, ``A Semi-Implicit Mid-Point Rule for Stiff
Systems of Ordinary Differential Equations.'', Numer. Math. 41, 373-398,
1983.
@end itemize