Mercurial > octave
diff src/lsode.cc @ 1:78fd87e624cb
[project @ 1993-08-08 01:13:40 by jwe]
Initial revision
author | jwe |
---|---|
date | Sun, 08 Aug 1993 01:13:40 +0000 |
parents | |
children | d68036bcad4c |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/lsode.cc Sun Aug 08 01:13:40 1993 +0000 @@ -0,0 +1,141 @@ +// tc-lsode.cc -*- C++ -*- +/* + +Copyright (C) 1993 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef __GNUG__ +#pragma implementation +#endif + +#include "ODE.h" + +#include "tree-const.h" +#include "variables.h" +#include "gripes.h" +#include "error.h" +#include "utils.h" + +// Global pointer for user defined function required by lsode. +static tree *lsode_fcn; + + +#ifdef WITH_DLD +tree_constant * +builtin_lsode_2 (tree_constant *args, int nargin, int nargout) +{ + return lsode (args, nargin, nargout); +} +#endif + +ColumnVector +lsode_user_function (const ColumnVector& x, double t) +{ + ColumnVector retval; + + int nstates = x.capacity (); + +// tree_constant name (lsode_fcn->name ()); + tree_constant *args = new tree_constant [3]; +// args[0] = name; + args[2] = tree_constant (t); + + if (nstates > 1) + { + Matrix m (nstates, 1); + for (int i = 0; i < nstates; i++) + m (i, 0) = x.elem (i); + tree_constant state (m); + args[1] = state; + } + else + { + double d = x.elem (0); + tree_constant state (d); + args[1] = state; + } + + if (lsode_fcn != NULL_TREE) + { + tree_constant *tmp = lsode_fcn->eval (args, 3, 1, 0); + delete [] args; + if (tmp != NULL_TREE_CONST && tmp[0].is_defined ()) + { + retval = tmp[0].to_vector (); + delete [] tmp; + } + else + { + delete [] tmp; + gripe_user_supplied_eval ("lsode"); + jump_to_top_level (); + } + } + + return retval; +} + +tree_constant * +lsode (tree_constant *args, int nargin, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + tree_constant *retval = NULL_TREE_CONST; + + lsode_fcn = is_valid_function (args[1], "lsode", 1); + if (lsode_fcn == NULL_TREE + || takes_correct_nargs (lsode_fcn, 3, "lsode", 1) != 1) + return retval; + + ColumnVector state = args[2].to_vector (); + ColumnVector out_times = args[3].to_vector (); + ColumnVector crit_times; + int crit_times_set = 0; + if (nargin > 4) + { + crit_times = args[4].to_vector (); + crit_times_set = 1; + } + + double tzero = out_times.elem (0); + int nsteps = out_times.capacity (); + + ODEFunc func (lsode_user_function); + ODE ode (state, tzero, func); + + int nstates = state.capacity (); + Matrix output (nsteps, nstates + 1); + + if (crit_times_set) + output = ode.integrate (out_times, crit_times); + else + output = ode.integrate (out_times); + + retval = new tree_constant [2]; + retval[0] = tree_constant (output); + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/