home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / liboctave / ODE.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  11KB  |  599 lines
  1. // ODE.cc                                                -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <math.h>
  29. #include <float.h>
  30. #include <iostream.h>
  31.  
  32. #include "ODE.h"
  33. #include "f77-uscore.h"
  34. #include "lo-error.h"
  35.  
  36. extern "C"
  37. {
  38.   int F77_FCN (lsode) (int (*)(int*, double*, double*, double*, int*),
  39.                int *, double *, double *, double *,
  40.                int *, double *, double *, int *, int *, int *,
  41.                double *, int *, int *, int *,
  42.                int (*)(int*, double*, double*, int*, int*,
  43.                    double*, int*), int *);
  44. }
  45.  
  46. static ODEFunc::ODERHSFunc user_fun;
  47. static ODEFunc::ODEJacFunc user_jac;
  48. static ColumnVector *tmp_x;
  49.  
  50. ODE::ODE (void)
  51. {
  52.   n = 0;
  53.   t = 0.0;
  54.  
  55.   stop_time_set = 0;
  56.   stop_time = 0.0;
  57.  
  58.   integration_error = 0;
  59.   restart = 1;
  60.  
  61.   istate = 1;
  62.   itol = 1;
  63.   itask = 1;
  64.   iopt = 0;
  65.  
  66.   liw = 20 + n;
  67.   lrw = 22 + n * (9 + n);
  68.  
  69.   iwork = new int [liw];
  70.   rwork = new double [lrw];
  71.   for (int i = 4; i < 9; i++)
  72.     {
  73.       iwork[i] = 0;
  74.       rwork[i] = 0.0;
  75.     }
  76.  
  77.   fun = 0;
  78.   jac = 0;
  79. }
  80.  
  81. ODE::ODE (int size)
  82. {
  83.   n = size;
  84.   t = 0.0;
  85.  
  86.   stop_time_set = 0;
  87.   stop_time = 0.0;
  88.  
  89.   integration_error = 0;
  90.   restart = 1;
  91.  
  92.   istate = 1;
  93.   itol = 1;
  94.   itask = 1;
  95.   iopt = 0;
  96.  
  97.   liw = 20 + n;
  98.   lrw = 22 + n * (9 + n);
  99.  
  100.   iwork = new int [liw];
  101.   rwork = new double [lrw];
  102.   for (int i = 4; i < 9; i++)
  103.     {
  104.       iwork[i] = 0;
  105.       rwork[i] = 0.0;
  106.     }
  107.  
  108.   fun = 0;
  109.   jac = 0;
  110. }
  111.  
  112. ODE::ODE (const ColumnVector& state, double time, const ODEFunc& f)
  113. {
  114.   n = state.capacity ();
  115.   t = time;
  116.   x = state;
  117.  
  118.   stop_time_set = 0;
  119.   stop_time = 0.0;
  120.  
  121.   integration_error = 0;
  122.   restart = 1;
  123.  
  124.   istate = 1;
  125.   itol = 1;
  126.   itask = 1;
  127.   iopt = 1;
  128.  
  129.   liw = 20 + n;
  130.   lrw = 22 + n * (9 + n);
  131.  
  132.   iwork = new int [liw];
  133.   rwork = new double [lrw];
  134.   for (int i = 4; i < 9; i++)
  135.     {
  136.       iwork[i] = 0;
  137.       rwork[i] = 0.0;
  138.     }
  139.  
  140.   fun = f.function ();
  141.   jac = f.jacobian_function ();
  142. }
  143.  
  144. ODE::~ODE (void)
  145. {
  146.   delete [] rwork;
  147.   delete [] iwork;
  148. }
  149.  
  150. int
  151. lsode_f (int *neq, double *time, double *state, double *deriv, int *ierr)
  152. {
  153.   int nn = *neq;
  154.   ColumnVector tmp_deriv (nn);
  155.  
  156.   /*
  157.    * NOTE: this won't work if LSODE passes copies of the state vector.
  158.    *       In that case we have to create a temporary vector object
  159.    *       and copy.
  160.    */
  161.   tmp_deriv = (*user_fun) (*tmp_x, *time);
  162.  
  163.   if (tmp_deriv.length () == 0)
  164.     *ierr = -1;
  165.   else
  166.     {
  167.       for (int i = 0; i < nn; i++)
  168.     deriv [i] = tmp_deriv.elem (i);
  169.     }
  170.  
  171.   return 0;
  172. }
  173.  
  174. int
  175. lsode_j (int *neq, double *time, double *state, int *ml, int *mu,
  176.          double *pd, int *nrowpd)
  177. {
  178.   int nn = *neq;
  179.   Matrix tmp_jac (nn, nn);
  180.  
  181.   /*
  182.    * NOTE: this won't work if LSODE passes copies of the state vector.
  183.    *       In that case we have to create a temporary vector object
  184.    *       and copy.
  185.    */
  186.   tmp_jac = (*user_jac) (*tmp_x, *time);
  187.  
  188.   for (int j = 0; j < nn; j++)
  189.     for (int i = 0; i < nn; i++)
  190.       pd [*nrowpd * j + i] = tmp_jac (i, j);
  191.  
  192.   return 0;
  193. }
  194.  
  195. ColumnVector
  196. ODE::integrate (double tout)
  197. {
  198.   if (jac)
  199.     method_flag = 21;
  200.   else
  201.     method_flag = 22;
  202.  
  203.   integration_error = 0;
  204.  
  205.   double *xp = x.fortran_vec ();
  206.  
  207. // NOTE: this won't work if LSODE passes copies of the state vector.
  208. //       In that case we have to create a temporary vector object
  209. //       and copy.
  210.  
  211.   tmp_x = &x;
  212.   user_fun = fun;
  213.   user_jac = jac;
  214.  
  215. // Try 5000 steps before giving up.
  216.  
  217.   iwork[5] = 5000;
  218.   int working_too_hard = 0;
  219.  
  220.   if (stop_time_set)
  221.     {
  222.       itask = 4;
  223.       rwork [0] = stop_time;
  224.     }
  225.   else
  226.     {
  227.       itask = 1;
  228.     }
  229.  
  230.   double abs_tol = absolute_tolerance ();
  231.   double rel_tol = relative_tolerance ();
  232.  
  233.   rwork[4] = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0;
  234.   rwork[5] = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0;
  235.   rwork[6] = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0;
  236.  
  237.   if (restart)
  238.     {
  239.       restart = 0;
  240.       istate = 1;
  241.     }
  242.  
  243.  again:
  244.  
  245.   (void) F77_FCN (lsode) (lsode_f, &n, xp, &t, &tout, &itol,
  246.               &rel_tol, &abs_tol, &itask, &istate, &iopt,
  247.               rwork, &lrw, iwork, &liw, lsode_j,
  248.               &method_flag);
  249.  
  250.   switch (istate)
  251.     {
  252.     case -13: // Return requested in user-supplied function.
  253.     case -6: // error weight became zero during problem. (solution
  254.          // component i vanished, and atol or atol(i) = 0.)
  255.     case -5: // repeated convergence failures (perhaps bad jacobian
  256.          // supplied or wrong choice of mf or tolerances).
  257.     case -4: // repeated error test failures (check all inputs).
  258.     case -3: // illegal input detected (see printed message).
  259.     case -2: // excess accuracy requested (tolerances too small).
  260.       integration_error = 1;
  261.       return ColumnVector ();
  262.       break;
  263.     case -1: // excess work done on this call (perhaps wrong mf).
  264.       working_too_hard++;
  265.       if (working_too_hard > 20)
  266.     {
  267.       (*current_liboctave_error_handler)
  268.         ("giving up after more than %d steps attempted in lsode",
  269.          iwork[5] * 20);
  270.       integration_error = 1;
  271.       return ColumnVector ();
  272.     }
  273.       else
  274.     {
  275.       istate = 2;
  276.       goto again;
  277.     }
  278.       break;
  279.     case 2: // lsode was successful
  280.       break;
  281.     default:
  282.       // Error?
  283.       break;
  284.     }
  285.  
  286.   t = tout;
  287.  
  288.   return x;
  289. }
  290.  
  291. void
  292. ODE::integrate (int nsteps, double tstep, ostream& s)
  293. {
  294.   int time_to_quit = 0;
  295.   double tout = t;
  296.  
  297.   s << t << " " << x << "\n";
  298.  
  299.   for (int i = 0; i < nsteps; i++)
  300.     {
  301.       tout += tstep;
  302.       if (stop_time_set && tout > stop_time)
  303.     {
  304.       tout = stop_time;
  305.       time_to_quit = 1;
  306.     }
  307.  
  308.       x = integrate (tout);
  309.  
  310.       s << t << " " << x << "\n";
  311.  
  312.       if (time_to_quit)
  313.     return;
  314.     }
  315. }
  316.  
  317. Matrix
  318. ODE::integrate (const ColumnVector& tout)
  319. {
  320.   Matrix retval;
  321.   int n_out = tout.capacity ();
  322.  
  323.   if (n_out > 0 && n > 0)
  324.     {
  325.       retval.resize (n_out, n);
  326.  
  327.       for (int i = 0; i < n; i++)
  328.     retval.elem (0, i) = x.elem (i);
  329.  
  330.       for (int j = 1; j < n_out; j++)
  331.     {
  332.       ColumnVector x_next = integrate (tout.elem (j));
  333.  
  334.       if (integration_error)
  335.         return retval;
  336.  
  337.       for (i = 0; i < n; i++)
  338.         retval.elem (j, i) = x_next.elem (i);
  339.     }
  340.     }
  341.  
  342.   return retval;
  343. }
  344.  
  345. Matrix
  346. ODE::integrate (const ColumnVector& tout, const ColumnVector& tcrit)
  347. {
  348.   Matrix retval;
  349.   int n_out = tout.capacity ();
  350.  
  351.   if (n_out > 0 && n > 0)
  352.     {
  353.       retval.resize (n_out, n);
  354.  
  355.       for (int i = 0; i < n; i++)
  356.     retval.elem (0, i) = x.elem (i);
  357.  
  358.       int n_crit = tcrit.capacity ();
  359.  
  360.       if (n_crit > 0)
  361.     {
  362.       int i_crit = 0;
  363.       int i_out = 1;
  364.       double next_crit = tcrit.elem (0);
  365.       double next_out;
  366.       while (i_out < n_out)
  367.         {
  368.           int do_restart = 0;
  369.  
  370.           next_out = tout.elem (i_out);
  371.           if (i_crit < n_crit)
  372.         next_crit = tcrit.elem (i_crit);
  373.  
  374.           int save_output;
  375.           double t_out;
  376.  
  377.           if (next_crit == next_out)
  378.         {
  379.           set_stop_time (next_crit);
  380.           t_out = next_out;
  381.           save_output = 1;
  382.           i_out++;
  383.           i_crit++;
  384.           do_restart = 1;
  385.         }
  386.           else if (next_crit < next_out)
  387.         {
  388.           if (i_crit < n_crit)
  389.             {
  390.               set_stop_time (next_crit);
  391.               t_out = next_crit;
  392.               save_output = 0;
  393.               i_crit++;
  394.               do_restart = 1;
  395.             }
  396.           else
  397.             {
  398.               clear_stop_time ();
  399.               t_out = next_out;
  400.               save_output = 1;
  401.               i_out++;
  402.             }
  403.         }
  404.           else
  405.         {
  406.           set_stop_time (next_crit);
  407.           t_out = next_out;
  408.           save_output = 1;
  409.           i_out++;
  410.         }
  411.  
  412.           ColumnVector x_next = integrate (t_out);
  413.  
  414.           if (integration_error)
  415.         return retval;
  416.  
  417.           if (save_output)
  418.         {
  419.           for (i = 0; i < n; i++)
  420.             retval.elem (i_out-1, i) = x_next.elem (i);
  421.         }
  422.  
  423.           if (do_restart)
  424.         force_restart ();
  425.         }
  426.     }
  427.       else
  428.     {
  429.       retval = integrate (tout);
  430.  
  431.       if (integration_error)
  432.         return retval;
  433.     }
  434.     }
  435.  
  436.   return retval;
  437. }
  438.  
  439. int
  440. ODE::size (void) const
  441. {
  442.   return n;
  443. }
  444.  
  445. ColumnVector
  446. ODE::state (void) const
  447. {
  448.   return x;
  449. }
  450.  
  451. double ODE::time (void) const
  452. {
  453.   return t;
  454. }
  455.  
  456. void
  457. ODE::force_restart (void)
  458. {
  459.   restart = 1;
  460. }
  461.  
  462. void
  463. ODE::initialize (const ColumnVector& state, double time)
  464. {
  465.   restart = 1;
  466.   x = state;
  467.   t = time;
  468. }
  469.  
  470. void
  471. ODE::set_stop_time (double time)
  472. {
  473.   stop_time_set = 1;
  474.   stop_time = time;
  475. }
  476.  
  477. void
  478. ODE::clear_stop_time (void)
  479. {
  480.   stop_time_set = 0;
  481. }
  482.  
  483. ODE_options::ODE_options (void)
  484. {
  485.   init ();
  486. }
  487.  
  488. ODE_options::ODE_options (const ODE_options& opt)
  489. {
  490.   copy (opt);
  491. }
  492.  
  493. ODE_options&
  494. ODE_options::operator = (const ODE_options& opt)
  495. {
  496.   if (this != &opt)
  497.     copy (opt);
  498.  
  499.   return *this;
  500. }
  501.  
  502. ODE_options::~ODE_options (void)
  503. {
  504. }
  505.  
  506. void
  507. ODE_options::init (void)
  508. {
  509.   double sqrt_eps = sqrt (DBL_EPSILON);
  510.   x_absolute_tolerance = sqrt_eps;
  511.   x_initial_step_size = -1.0;
  512.   x_maximum_step_size = -1.0;
  513.   x_minimum_step_size = 0.0;
  514.   x_relative_tolerance = sqrt_eps;
  515. }
  516.  
  517. void
  518. ODE_options::copy (const ODE_options& opt)
  519. {
  520.   x_absolute_tolerance = opt.x_absolute_tolerance;
  521.   x_initial_step_size = opt.x_initial_step_size;
  522.   x_maximum_step_size = opt.x_maximum_step_size;
  523.   x_minimum_step_size = opt.x_minimum_step_size;
  524.   x_relative_tolerance = opt.x_relative_tolerance;
  525. }
  526.  
  527. void
  528. ODE_options::set_default_options (void)
  529. {
  530.   init ();
  531. }
  532.  
  533. void
  534. ODE_options::set_absolute_tolerance (double val)
  535. {
  536.   x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  537. }
  538.  
  539. void
  540. ODE_options::set_initial_step_size (double val)
  541. {
  542.   x_initial_step_size = (val >= 0.0) ? val : -1.0;
  543. }
  544.  
  545. void
  546. ODE_options::set_maximum_step_size (double val)
  547. {
  548.   x_maximum_step_size = (val >= 0.0) ? val : -1.0;
  549. }
  550.  
  551. void
  552. ODE_options::set_minimum_step_size (double val)
  553. {
  554.   x_minimum_step_size = (val >= 0.0) ? val : 0.0;
  555. }
  556.  
  557. void
  558. ODE_options::set_relative_tolerance (double val)
  559. {
  560.   x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  561. }
  562.  
  563. double
  564. ODE_options::absolute_tolerance (void)
  565. {
  566.   return x_absolute_tolerance;
  567. }
  568.  
  569. double
  570. ODE_options::initial_step_size (void)
  571. {
  572.   return x_initial_step_size;
  573. }
  574.  
  575. double
  576. ODE_options::maximum_step_size (void)
  577. {
  578.   return x_maximum_step_size;
  579. }
  580.  
  581. double
  582. ODE_options::minimum_step_size (void)
  583. {
  584.   return x_minimum_step_size;
  585. }
  586.  
  587. double
  588. ODE_options::relative_tolerance (void)
  589. {
  590.   return x_relative_tolerance;
  591. }
  592.  
  593. /*
  594. ;;; Local Variables: ***
  595. ;;; mode: C++ ***
  596. ;;; page-delimiter: "^/\\*" ***
  597. ;;; End: ***
  598. */
  599.