Mercurial > octave
changeset 22901:4c56f3ffec04
Add functions ode15i and ode15s
* libinterp/dldfcn/__ode15__.cc: Add oct-file backend for
m-files ode15i.m and ode15s.m.
* scripts/ode/ode15i.m: Add wrapper function ode15i.m.
* scripts/ode/ode15s.m: Add wrapper function ode15s.m.
* scripts/ode/private/check_default_input.m: Add private
function to check default input of ode15i and ode15s.
* libinterp/dldfcn/module-files: Add __ode15__.cc and its flags.
| author | Francesco Faccio <francesco.faccio@mail.polimi.it> |
|---|---|
| date | Tue, 23 Aug 2016 03:19:11 +0200 |
| parents | ee1c77705fcd |
| children | 284bbb0328f2 |
| files | libinterp/dldfcn/__ode15__.cc libinterp/dldfcn/module-files scripts/ode/decic.m scripts/ode/ode15i.m scripts/ode/ode15s.m scripts/ode/private/check_default_input.m |
| diffstat | 6 files changed, 2617 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/dldfcn/__ode15__.cc Tue Aug 23 03:19:11 2016 +0200 @@ -0,0 +1,1321 @@ +/* + +Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> + +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 3 of the License, 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, see +<http://www.gnu.org/licenses/>. + +*/ + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#if defined (HAVE_IDA_IDA_H) +#include <ida/ida.h> +#include <ida/ida_dense.h> +#include <ida/ida_klu.h> +#include <sundials/sundials_sparse.h> +#endif + +#if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H) +#include <nvector/nvector_serial.h> +#endif + + +#include <oct.h> +#include <parse.h> +#include <time.h> + +#include "defun-dld.h" +#include "iostream" +#include "error.h" +#include "ov.h" +#include "ov-struct.h" +#include "errwarn.h" + +#include "algorithm" +#include "iterator" + + +#if defined (HAVE_SUNDIALS) + +namespace octave +{ + class IDA + { + public: + + typedef + ColumnVector (*DAERHSFuncIDA) (const ColumnVector& x, + const ColumnVector& xdot, + realtype t, octave_function *idaf); + + typedef + Matrix (*DAEJacFuncDense) (const ColumnVector& x, + const ColumnVector& xdot, realtype t, + realtype cj, octave_function *idaj); + + typedef + SparseMatrix (*DAEJacFuncSparse) (const ColumnVector& x, + const ColumnVector& xdot, + realtype t, realtype cj, + octave_function *idaj); + + typedef + Matrix (*DAEJacCellDense) (Matrix *dfdy, Matrix *dfdyp, + realtype cj); + + typedef + SparseMatrix (*DAEJacCellSparse) (SparseMatrix *dfdy, + SparseMatrix *dfdyp, realtype cj); + + //Default + IDA (void) + : t0 (0.0), y0 (), yp0 (), havejac (false), havejacfun (false), + havejacsparse (false), mem (NULL), num (), ida_fun (NULL), + ida_jac (NULL), dfdy (NULL), dfdyp (NULL), spdfdy (NULL), + spdfdyp (NULL), fun (0), jacfun (NULL), jacspfun (0), + jacdcell (0), jacspcell (0) + { }; + + + IDA (realtype t, ColumnVector y, ColumnVector yp, + octave_function *ida_fcn, DAERHSFuncIDA daefun) + : t0 (t), y0 (y), yp0 (yp), havejac (false), havejacfun (false), + havejacsparse (false), mem (NULL), num (), ida_fun (ida_fcn), + ida_jac (NULL), dfdy (NULL), dfdyp (NULL), spdfdy (NULL), + spdfdyp (NULL), fun (daefun), jacfun (NULL), jacspfun (0), + jacdcell (0), jacspcell (0) + { }; + + + ~IDA (void) { IDAFree (&mem); } + + IDA& + set_jacobian (octave_function *jac, DAEJacFuncDense j) + { + jacfun = j; + ida_jac = jac; + havejac = true; + havejacfun = true; + havejacsparse = false; + return *this; + } + + IDA& + set_jacobian (octave_function *jac, DAEJacFuncSparse j) + { + jacspfun = j; + ida_jac = jac; + havejac = true; + havejacfun = true; + havejacsparse = true; + return *this; + } + + IDA& + set_jacobian (Matrix *dy, Matrix *dyp, DAEJacCellDense j) + { + jacdcell = j; + dfdy = dy; + dfdyp = dyp; + havejac = true; + havejacfun = false; + havejacsparse = false; + return *this; + } + + IDA& + set_jacobian (SparseMatrix *dy, SparseMatrix *dyp, + DAEJacCellSparse j) + { + jacspcell = j; + spdfdy = dy; + spdfdyp = dyp; + havejac = true; + havejacfun = false; + havejacsparse = true; + return *this; + } + + void + set_userdata (void); + + void + initialize (void); + + static ColumnVector + NVecToCol (N_Vector& v, long int n); + + static N_Vector + ColToNVec (const ColumnVector& data, long int n); + + void + set_up (void); + + void + set_tolerance (ColumnVector& abstol, realtype reltol); + + void + set_tolerance (realtype abstol, realtype reltol); + + static int + resfun (realtype t, N_Vector yy, N_Vector yyp, + N_Vector rr, void *user_data); + + void + resfun_impl (realtype t, N_Vector& yy, + N_Vector& yyp, N_Vector& rr); + + static int + jacdense (long int Neq, realtype t, realtype cj, + N_Vector yy, N_Vector yyp, N_Vector resvec, + DlsMat JJ, void *user_data, N_Vector tempv1, + N_Vector tempv2, N_Vector tempv3) + { + IDA *self = + static_cast <IDA *> (user_data); + + self -> jacdense_impl (Neq, t, cj, yy, yyp, JJ); + return 0; + } + + void + jacdense_impl (long int Neq, realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, DlsMat& JJ); + + static int + jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp, + N_Vector r, SlsMat Jac, void *user_data, N_Vector tmp1, + N_Vector tmp2, N_Vector tmp3) + { + IDA *self = + static_cast <IDA *> (user_data); + + self -> jacsparse_impl (t, cj, yy, yyp, Jac); + return 0; + } + + void + jacsparse_impl (realtype t, realtype cj, N_Vector& yy, + N_Vector& yyp, SlsMat& Jac); + + void + set_maxstep (realtype maxstep); + + void + set_initialstep (realtype initialstep); + + bool + interpolate (int& cont, Matrix& output, ColumnVector& tout, + int refine, realtype tend, bool haveoutputfcn, + bool haveoutputsel, octave_function *output_fcn, + ColumnVector& outputsel, bool haveeventfunction, + octave_function * event_fcn, ColumnVector& te, + Matrix& ye, ColumnVector& ie, ColumnVector& oldval, + ColumnVector& oldisterminal, ColumnVector& olddir, + int& temp, ColumnVector& yold); + + bool + outputfun (octave_function *output_fcn, bool haveoutputsel, + const ColumnVector& output, realtype tout, realtype tend, + ColumnVector& outputsel, std::string flag); + + + bool + event (octave_function *event_fcn, + ColumnVector& te, Matrix& ye, ColumnVector& ie, + realtype tsol, const ColumnVector& y, std::string flag, + const ColumnVector& yp, ColumnVector& oldval, + ColumnVector& oldisterminal, ColumnVector& olddir, + int cont, int& temp, realtype told, ColumnVector& yold); + + void + set_maxorder (int maxorder); + + octave_value_list + integrate (const int numt, const ColumnVector& tt, + const ColumnVector& y0, const ColumnVector& yp0, + const int refine, bool haverefine, bool haveoutputfcn, + octave_function *output_fcn, bool haveoutputsel, + ColumnVector& outputsel, bool haveeventfunction, + octave_function *event_fcn); + + void + print_stat (void); + + private: + + realtype t0; + ColumnVector y0; + ColumnVector yp0; + bool havejac; + bool havejacfun; + bool havejacsparse; + void *mem; + int num; + octave_function *ida_fun; + octave_function *ida_jac; + Matrix *dfdy; + Matrix *dfdyp; + SparseMatrix *spdfdy; + SparseMatrix *spdfdyp; + DAERHSFuncIDA fun; + DAEJacFuncDense jacfun; + DAEJacFuncSparse jacspfun; + DAEJacCellDense jacdcell; + DAEJacCellSparse jacspcell; + + }; + + int + IDA::resfun (realtype t, N_Vector yy, N_Vector yyp, N_Vector rr, + void *user_data) + { + IDA *self = + static_cast <IDA *> (user_data); + + self -> resfun_impl (t, yy, yyp, rr); + return 0; + } + + void + IDA::resfun_impl (realtype t, N_Vector& yy, + N_Vector& yyp, N_Vector& rr) + { + BEGIN_INTERRUPT_WITH_EXCEPTIONS; + + ColumnVector y = + IDA::NVecToCol (yy, num); + + ColumnVector yp = + IDA::NVecToCol (yyp, num); + + ColumnVector res = + (*fun) (y, yp, t, ida_fun); + + realtype *puntrr = + NV_DATA_S (rr); + + for (octave_idx_type i = 0; i < num; i++) + puntrr [i] = res (i); + + END_INTERRUPT_WITH_EXCEPTIONS; + } + + void + IDA::set_up (void) + { + if (havejacsparse) + { + int flag = + IDAKLU (mem, num, num*num); + if (flag != 0) + { + error ("IDAKLU solver not initialized"); + } + flag = IDASlsSetSparseJacFn (mem, IDA::jacsparse); + } + else + { + int flag = + IDADense (mem, num); + if (flag != 0) + { + error ("IDADense solver not initialized"); + } + if (havejac) + { + flag = IDADlsSetDenseJacFn (mem, IDA::jacdense); + if (flag != 0) + { + error ("Dense Jacobian not set"); + } + } + } + } + + void + IDA::jacdense_impl (long int Neq, realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, DlsMat& JJ) + + { + BEGIN_INTERRUPT_WITH_EXCEPTIONS; + + ColumnVector y = + NVecToCol (yy, Neq); + + ColumnVector yp = + NVecToCol (yyp, Neq); + + Matrix jac; + + if (havejacfun) + jac = (*jacfun) (y, yp, t, cj, ida_jac); + else + jac = (*jacdcell) (dfdy, dfdyp, cj); + + std::copy (jac.fortran_vec (), + jac.fortran_vec () + jac.numel (), + JJ -> data); + + END_INTERRUPT_WITH_EXCEPTIONS; + } + + void + IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, + SlsMat& Jac) + + { + BEGIN_INTERRUPT_WITH_EXCEPTIONS; + + ColumnVector y = + NVecToCol (yy, num); + + ColumnVector yp = + NVecToCol (yyp, num); + + SparseMatrix jac; + + if (havejacfun) + jac = (*jacspfun) (y, yp, t, cj, ida_jac); + else + jac = (*jacspcell) (spdfdy, spdfdyp, cj); + + SlsSetToZero (Jac); + + for (int i = 0; i < num + 1; i++) + Jac -> colptrs [i] = jac.cidx (i); + + for (int i = 0; i < jac.nnz (); i++) + { + Jac -> rowvals [i] = jac.ridx (i); + Jac -> data [i] = jac.data (i); + } + + END_INTERRUPT_WITH_EXCEPTIONS; + } + + ColumnVector + IDA::NVecToCol (N_Vector& v, long int n) + { + ColumnVector data (n); + realtype *punt; + punt = NV_DATA_S (v); + + for (octave_idx_type i = 0; i < n; i++) + data (i) = punt [i]; + + return data; + } + + N_Vector + IDA::ColToNVec (const ColumnVector &data, long int n) + { + N_Vector v = + N_VNew_Serial (n); + + realtype * punt; + punt = NV_DATA_S (v); + + for (octave_idx_type i = 0; i < n; i++) + punt [i] = data (i); + + return v; + } + + void + IDA::set_userdata (void) + { + void * userdata = this; + + int flag = + IDASetUserData(mem, userdata); + if (flag != 0) + { + error ("User data not set"); + } + } + + void + IDA::initialize (void) + { + num = y0.numel(); + mem = IDACreate (); + + N_Vector yy = + ColToNVec(y0, num); + + N_Vector yyp = + ColToNVec(yp0, num); + + IDA::set_userdata (); + + int flag = + IDAInit (mem, IDA::resfun, t0, yy, yyp); + if (flag != 0) + { + error ("IDA not initialized"); + } + } + + void + IDA::set_tolerance (ColumnVector& abstol, realtype reltol) + { + N_Vector abs_tol = + ColToNVec (abstol, num); + + int flag = + IDASVtolerances (mem, reltol, abs_tol); + if (flag != 0) + { + error ("IDA: Tolerance not set"); + } + N_VDestroy_Serial (abs_tol); + } + + void + IDA::set_tolerance (realtype abstol, realtype reltol) + { + int flag = + IDASStolerances (mem, reltol, abstol); + if (flag != 0) + { + error ("IDA: Tolerance not set"); + } + } + + octave_value_list + IDA::integrate (const int numt, const ColumnVector& tspan, + const ColumnVector& y, const ColumnVector& yp, + const int refine, bool haverefine, bool haveoutputfcn, + octave_function *output_fcn, bool haveoutputsel, + ColumnVector& outputsel, bool haveeventfunction, + octave_function *event_fcn) + { + // Set up output + ColumnVector tout, yout (num), ypout (num), ysel (outputsel.numel ()); + ColumnVector ie, te, oldval, oldisterminal, olddir; + Matrix output, ye; + int cont = 0, flag = 0, temp = 0; + bool status = 0; + std::string string = ""; + ColumnVector yold = y; + octave_value_list retval; + + realtype tsol = + tspan (0); + realtype tend = + tspan (numt - 1); + + N_Vector yyp = + ColToNVec (yp, num); + + N_Vector yy = + ColToNVec (y, num); + + // Initialize OutputFcn + if (haveoutputfcn) + status = IDA::outputfun (output_fcn, haveoutputsel, y, + tsol, tend, outputsel, "init"); + + // Initialize Events + if (haveeventfunction) + status = IDA::event (event_fcn, te, ye, ie, tsol, y, + "init", yp, oldval, oldisterminal, + olddir, cont, temp, tsol, yold); + + if (numt > 2) + { + // First output value + tout.resize (numt); + tout (0) = tsol; + output.resize (numt, num); + + for (octave_idx_type i = 0; i < num; i++) + output.elem (0, i) = y.elem (i); + + //Main loop + for (octave_idx_type j = 1; j < numt && status == 0; j++) + { + // IDANORMAL already interpolates tspan(j) + flag = IDASolve (mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL); + if (flag != 0) + { + error ("IDASolve failed"); + } + + yout = NVecToCol (yy, num); + ypout = NVecToCol (yyp, num); + tout (j) = tsol; + + for (octave_idx_type i = 0; i < num; i++) + output.elem (j, i) = yout.elem (i); + + if (haveoutputfcn) + status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, + tend, outputsel, string); + + if (haveeventfunction) + status = IDA::event (event_fcn, te, ye, ie, tout (j), yout, + string, ypout, oldval, oldisterminal, + olddir, j, temp, tout (j - 1), yold); + + // If integration is stopped, return only the reached steps + if (status == 1) + { + output.resize (j + 1, num); + tout.resize (j + 1); + } + + } + } + else // numel (tspan) == 2 + { + // First output value + tout.resize (1); + tout (0) = tsol; + output.resize (1, num); + + for (octave_idx_type i = 0; i < num; i++) + output.elem (0, i) = y.elem (i); + + bool posdirection = + (tend > tsol); + + //main loop + while (((posdirection == 1 && tsol < tend) + || (posdirection == 0 && tsol > tend)) + && status == 0) + { + flag = IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP); + if (flag != 0) + { + error ("IDASolve failed"); + } + + if (haverefine) + status = IDA::interpolate (cont, output, tout, refine, tend, + haveoutputfcn, haveoutputsel, + output_fcn, outputsel, + haveeventfunction, event_fcn, te, + ye, ie, oldval, oldisterminal, + olddir, temp, yold); + + ypout = NVecToCol (yyp, num); + cont += 1; + output.resize (cont + 1, num); // This may be not efficient + tout.resize (cont + 1); + tout (cont) = tsol; + yout = NVecToCol (yy, num); + + for (octave_idx_type i = 0; i < num; i++) + output.elem (cont, i) = yout.elem (i); + + if (haveoutputfcn && ! haverefine && tout (cont) < tend) + status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, + tend, outputsel, string); + + if (haveeventfunction && ! haverefine && tout (cont) < tend) + status = IDA::event (event_fcn, te, ye, ie, tout (cont), yout, + string, ypout, oldval, oldisterminal, + olddir, cont, temp, tout (cont - 1), yold); + } + if (status == 0) + { + // Interpolate in tend + N_Vector dky = + N_VNew_Serial (num); + + flag = IDAGetDky (mem, tend, 0, dky); + if (flag != 0) + { + error ("IDA failed to interpolate y"); + } + + tout (cont) = tend; + yout = NVecToCol (dky, num); + + for (octave_idx_type i = 0; i < num; i++) + output.elem (cont, i) = yout.elem (i); + + // Plot final value + if (haveoutputfcn) + { + status = IDA::outputfun (output_fcn, haveoutputsel, yout, + tend, tend, outputsel, string); + + // Events during last step + if (haveeventfunction) + status = IDA::event (event_fcn, te, ye, ie, tend, yout, + string, ypout, oldval, oldisterminal, + olddir, cont, temp, tout (cont - 1), + yold); + } + + N_VDestroy_Serial (dky); + } + // Cleanup plotter + status = IDA::outputfun (output_fcn, haveoutputsel, yout, tend, tend, + outputsel, "done"); + + } + + retval (0) = tout; + retval (1) = output; + retval (2) = te; + retval (3) = ye; + retval (4) = ie; + + return retval; + } + + bool + IDA::event (octave_function *event_fcn, + ColumnVector& te, Matrix& ye, ColumnVector& ie, + realtype tsol, const ColumnVector& y, std::string flag, + const ColumnVector& yp, ColumnVector& oldval, + ColumnVector& oldisterminal, ColumnVector& olddir, int cont, + int& temp, realtype told, ColumnVector& yold) + { + bool status = 0; + octave_value_list args (3); + octave_value_list output (3); + ColumnVector val, isterminal, dir; + args (0) = tsol; + args (1) = y; + args (2) = yp; + + // cont is the number of steps reached by the solver + // temp is the number of events registered + + if (flag == "init") + { + output = feval (event_fcn, args, 3); + oldval = output(0).vector_value (); + oldisterminal = output(1).vector_value (); + olddir = output(2).vector_value (); + } + else if (flag == "") + { + ColumnVector index (0); + output = feval (event_fcn, args, 3); + val = output(0).vector_value (); + isterminal = output(1).vector_value (); + dir = output(2).vector_value (); + + // Get the index of the changed values + for (octave_idx_type i = 0; i < val.numel (); i++) + { + if ((val(i) > 0 && oldval(i) < 0 && dir(i) != -1) // increasing + || (val(i) < 0 && oldval(i) > 0 && dir(i) != 1)) // decreasing + { + index.resize (index.numel () + 1); + index (index.numel () - 1) = i; + } + } + + if (cont == 1 && index.numel () > 0) // Events in first step + { + temp = 1; // register only the first event + te.resize (1); + ye.resize (1, num); + ie.resize (1); + + // Linear interpolation + ie (0) = index (0); + te (0) = tsol - val (index (0)) * (tsol - told) + / (val (index (0)) - oldval (index (0))); + + ColumnVector ytemp = + y - ((tsol - te (0)) * (y - yold ) / (tsol - told)); + + for (octave_idx_type i = 0; i < num; i++) + ye.elem (0, i) = ytemp.elem (i); + + } + else if (index.numel () > 0) + // Not first step: register all events and test if stop integration or not + { + te.resize (temp + index.numel ()); + ye.resize (temp + index.numel (), num); + ie.resize (temp + index.numel ()); + + for (octave_idx_type i = 0; i < index.numel (); i++) + { + + if (isterminal (index (i)) == 1) + status = 1; // Stop integration + + // Linear interpolation + ie (temp + i) = index (i); + te (temp + i) = tsol - + val (index (i)) * (tsol - told) / + (val (index (i)) - oldval (index (i))); + + ColumnVector ytemp = + y - (tsol - te (temp + i)) * (y - yold) / (tsol - told); + + for (octave_idx_type j = 0; j < num; j++) + ye.elem (temp + i, j) = ytemp.elem (j); + + } + + temp += index.numel (); + } + + // Update variables + yold = y; + told = tsol; + olddir = dir; + oldval = val; + oldisterminal = isterminal; + } + + return status; + } + + bool + IDA::interpolate (int& cont, Matrix& output, ColumnVector& tout, + int refine, realtype tend, bool haveoutputfcn, + bool haveoutputsel, octave_function *output_fcn, + ColumnVector& outputsel, bool haveeventfunction, + octave_function * event_fcn, ColumnVector& te, + Matrix& ye, ColumnVector& ie, ColumnVector& oldval, + ColumnVector& oldisterminal, ColumnVector& olddir, + int& temp, ColumnVector& yold) + { + realtype h = 0, tcur = 0; + bool status = 0; + + N_Vector dky = + N_VNew_Serial (num); + + N_Vector dkyp = + N_VNew_Serial (num); + + ColumnVector yout (num); + ColumnVector ypout (num); + std::string string = ""; + + int flag = + IDAGetLastStep (mem, &h); + + if (flag != 0) + { + error ("IDA failed to return last step"); + } + flag = IDAGetCurrentTime (mem, &tcur); + if (flag != 0) + { + error ("IDA failed to return the current time"); + } + + realtype tin = + tcur - h; + + realtype step = + h / refine; + + for (octave_idx_type i = 1; i < refine && tin + step * i < tend + && status == 0; i++) + { + flag = IDAGetDky (mem, tin + step*i, 0, dky); + if (flag != 0) + { + error ("IDA failed to interpolate y"); + } + flag = IDAGetDky (mem, tin + step*i, 1, dkyp); + if (flag != 0) + { + error ("IDA failed to interpolate yp"); + } + cont += 1; + output.resize (cont + 1, num); + tout.resize (cont + 1); + + tout (cont) = tin + step * i; + yout = NVecToCol (dky, num); + ypout = NVecToCol (dkyp, num); + + for (octave_idx_type j = 0; j < num; j++) + output.elem (cont, j) = yout.elem (j); + + if (haveoutputfcn) + status = IDA::outputfun (output_fcn, haveoutputsel, yout, + tout (cont), tend, outputsel, ""); + + if (haveeventfunction) + status = IDA::event (event_fcn, te, ye, ie, tout (cont), + yout, string, ypout, oldval, + oldisterminal, olddir, cont, temp, + tout (cont - 1), yold); + } + N_VDestroy_Serial (dky); + return status; + } + + bool + IDA::outputfun (octave_function *output_fcn, bool haveoutputsel, + const ColumnVector& yout, realtype tsol, + realtype tend, ColumnVector& outputsel, + std::string flag) + { + bool status = 0; + octave_value_list output, val; + ColumnVector ysel (outputsel.numel ()); + + if (haveoutputsel) + { + for (octave_idx_type i = 0; i < outputsel.numel (); i++) + ysel (i) = yout (outputsel (i)); + + output (1) = ysel; + } + else + output (1) = yout; + + output (2) = flag; + + if (flag == "init") + { + ColumnVector toutput (2); + toutput (0) = tsol; + toutput (1) = tend; + output (0) = toutput; + + feval (output_fcn, output, 0); + } + else if (flag == "") + { + output (0) = tsol; + val = feval (output_fcn, output, 1); + status = val(0).bool_value (); + } + else + { // Cleanup plotter + output (0) = tend; + feval (output_fcn, output, 0); + } + + return status; + } + + void + IDA::set_maxstep (realtype maxstep) + { + int flag = + IDASetMaxStep (mem, maxstep); + + if (flag != 0) + { + error ("IDA: Max Step not set"); + } + } + + void + IDA::set_initialstep (realtype initialstep) + { + int flag = + IDASetInitStep (mem, initialstep); + + if (flag != 0) + { + error ("IDA: Initial Step not set"); + } + } + + void + IDA::set_maxorder (int maxorder) + { + int flag = + IDASetMaxOrd (mem, maxorder); + + if (flag != 0) + { + error ("IDA: Max Order not set"); + } + } + + void + IDA::print_stat (void) + { + long int nsteps = 0, netfails = 0, nrevals = 0; + int flag = + IDAGetNumSteps(mem, &nsteps); + + if (flag != 0) + { + error ("IDA failed to return the number of internal steps"); + } + flag = IDAGetNumErrTestFails(mem, &netfails); + if (flag != 0) + { + error ("IDA failed to return the number of internal errors"); + } + flag = IDAGetNumResEvals(mem, &nrevals); + if (flag != 0) + { + error ("IDA failed to return the number of residual evaluations"); + } + + std::cout<<nsteps<<" successful steps\n"; + std::cout<<netfails<<" failed attempts\n"; + std::cout<<nrevals<<" function evaluations\n"; + //std::cout<<<<" partial derivatives\n"; + //std::cout<<<<" LU decompositions\n"; + //std::cout<<<<" solutions of linear systems\n"; + } + + ColumnVector + ida_user_function (const ColumnVector& x, const ColumnVector& xdot, + double t, octave_function *ida_fc) + { + ColumnVector retval; + octave_value_list args; + + args(2) = xdot; + args(1) = x; + args(0) = t; + + octave_value_list tmp; + + try + { + tmp = ida_fc -> do_multi_index_op (1, args); + } + catch (octave::execution_exception& e) + { + err_user_supplied_eval (e, "__ode15__"); + } + + retval = tmp(0).vector_value (); + + return retval; + } + + Matrix + ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot, + double t, double cj, octave_function *ida_jc) + { + Matrix retval; + octave_value_list newargs (3); + + newargs (0) = t; + newargs (1) = x; + newargs (2) = xdot; + + octave_value_list tmp; + + try + { + tmp = ida_jc -> do_multi_index_op (2, newargs); + } + catch (octave::execution_exception& e) + { + err_user_supplied_eval (e, "__ode15__"); + } + + retval = tmp(0).matrix_value () + cj * tmp(1).matrix_value (); + + return retval; + } + + SparseMatrix + ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot, + double t, double cj, octave_function *ida_jc) + { + SparseMatrix retval; + octave_value_list newargs (3); + + newargs (0) = t; + newargs (1) = x; + newargs (2) = xdot; + + octave_value_list tmp; + + try + { + tmp = ida_jc -> do_multi_index_op (2, newargs); + } + catch (octave::execution_exception& e) + { + err_user_supplied_eval (e, "__ode15__"); + } + + retval = tmp(0).sparse_matrix_value () + + cj * tmp(1).sparse_matrix_value (); + + return retval; + } + + Matrix + ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj) + { + Matrix retval; + retval = (*dfdy) + cj * (*dfdyp); + return retval; + } + + SparseMatrix + ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp, + double cj) + { + SparseMatrix retval; + retval = (*spdfdy) + cj * (*spdfdyp); + return retval; + } + + octave_value_list + do_ode15 (octave_function *ida_fcn, + const ColumnVector &tspan, + const int numt, + const realtype t0, + const ColumnVector &y0, + const ColumnVector &yp0, + const octave_scalar_map &options) + { + octave_value_list retval; + + // Create object + IDA dae (t0, y0, yp0, ida_fcn, ida_user_function); + + // Set Jacobian + bool havejac = + options.getfield ("havejac").bool_value (); + + bool havejacsparse = + options.getfield ("havejacsparse").bool_value (); + + bool havejacfun = + options.getfield ("havejacfun").bool_value (); + + Matrix ida_dfdy, ida_dfdyp, *dfdy, *dfdyp; + SparseMatrix ida_spdfdy, ida_spdfdyp, *spdfdy, *spdfdyp; + octave_function *ida_jac; + Cell jaccell; + + if (havejac) + { + if (havejacfun) + { + ida_jac = options.getfield ("Jacobian").function_value (); + if (havejacsparse) + dae.set_jacobian (ida_jac, ida_sparse_jac); + else + dae.set_jacobian (ida_jac, ida_dense_jac); + } + else + { + jaccell = options.getfield ("Jacobian").cell_value (); + + if (havejacsparse) + { + ida_spdfdy = jaccell(0).sparse_matrix_value (); + ida_spdfdyp = jaccell(1).sparse_matrix_value (); + spdfdy = &ida_spdfdy; + spdfdyp = &ida_spdfdyp; + dae.set_jacobian (spdfdy, spdfdyp, ida_sparse_cell_jac); + } + else + { + ida_dfdy = jaccell(0).matrix_value (); + ida_dfdyp = jaccell(1).matrix_value (); + dfdy = &ida_dfdy; + dfdyp = &ida_dfdyp; + dae.set_jacobian (dfdy, dfdyp, ida_dense_cell_jac); + } + } + } + + // Initialize IDA + dae.initialize (); + + // Set tolerances + realtype rel_tol = + options.getfield("RelTol").double_value (); + + bool haveabstolvec = + options.getfield ("haveabstolvec").bool_value (); + + if (haveabstolvec) + { + ColumnVector abs_tol = + options.getfield("AbsTol").vector_value (); + + dae.set_tolerance (abs_tol, rel_tol); + } + else + { + realtype abs_tol = + options.getfield("AbsTol").double_value (); + + dae.set_tolerance (abs_tol, rel_tol); + } + + //Set max step + realtype maxstep = + options.getfield("MaxStep").double_value (); + + dae.set_maxstep (maxstep); + + //Set initial step + if (!(options.getfield("InitialStep").is_empty ())) + { + realtype initialstep = + options.getfield("InitialStep").double_value (); + + dae.set_initialstep (initialstep); + } + + //Set max order FIXME: it doesn't work + int maxorder = + options.getfield("MaxOrder").int_value (); + + dae.set_maxorder (maxorder); + + //Set Refine + const int refine = + options.getfield("Refine").int_value (); + + bool haverefine = + (refine > 1); + + octave_function *output_fcn = NULL; + ColumnVector outputsel; + + // OutputFcn + bool haveoutputfunction = + options.getfield("haveoutputfunction").bool_value (); + + if (haveoutputfunction) + output_fcn = options.getfield("OutputFcn").function_value (); + + // OutputSel + bool haveoutputsel = + options.getfield("haveoutputselection").bool_value (); + + if (haveoutputsel) + outputsel = options.getfield("OutputSel").vector_value (); + + octave_function *event_fcn = NULL; + + // Events + bool haveeventfunction = + options.getfield("haveeventfunction").bool_value (); + + if (haveeventfunction) + event_fcn = options.getfield("Events").function_value (); + + // Set up linear solver + dae.set_up (); + + // Integrate + retval = dae.integrate (numt, tspan, y0, yp0, refine, + haverefine, haveoutputfunction, + output_fcn, haveoutputsel, outputsel, + haveeventfunction, event_fcn); + + // Statistics + bool havestats = + options.getfield("havestats").bool_value (); + + if (havestats) + dae.print_stat (); + + return retval; + } +} +#endif + + +DEFUN_DLD (__ode15__, args, nargout, doc: /* -*- texinfo -*- +@deftypefn {} {@var{t}, @var{y} =} __ode15__ (@var{fun}, @ +@var{tspan}, @var{y0}, @var{yp0}, @var{options}) +Undocumented internal function. +@end deftypefn */) +{ + +#if defined (HAVE_SUNDIALS) + + // Check number of parameters + int nargin = args.length (); + + if (nargin != 5 || nargout != 5) + print_usage (); + + // Check odefun + octave_function *ida_fcn = NULL; + + octave_value f_arg = args(0); + + if (f_arg.is_function_handle ()) + ida_fcn = f_arg.function_value (); + else + error ("__ode15__: odefun must be a function handle"); + + // Check input tspan + ColumnVector tspan = + args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers"); + + int numt = + tspan.numel (); + + realtype t0 = + tspan (0); + + if (numt < 2) + error ("__ode15__: TRANGE must contain at least 2 elements"); + else if (!(tspan.is_sorted ()) || (tspan(0) == tspan(numt - 1))) + error ("__ode15__: TRANGE must be strictly monotonic"); + + // input y0 and yp0 + ColumnVector y0 = + args(2).xvector_value ("__ode15__: initial state y0 must be a vector"); + + ColumnVector yp0 = + args(3).xvector_value ("__ode15__: initial state yp0 must be a vector"); + + + if (y0.numel () != yp0.numel ()) + error ("__ode15__: initial state y0 and yp0 must have the same length"); + else if (y0.numel () < 1) + error ("__ode15__: initial state yp0 must be a vector or a scalar"); + + + if (! args(4).is_map ()) + error ("__ode15__: OPTS argument must be a structure"); + + octave_scalar_map options = + args(4).xscalar_map_value ("__ode15__:", + "OPTS argument must be a scalar structure"); + + + return octave::do_ode15 (ida_fcn, tspan, numt, t0, + y0, yp0, options); + + +#else + + octave_unused_parameter (args); + err_disabled_feature ("__ode15__", "sundials_ida, sundials_nvecserial"); + +#endif +} + + +
--- a/libinterp/dldfcn/module-files Tue Aug 23 02:31:51 2016 +0200 +++ b/libinterp/dldfcn/module-files Tue Aug 23 03:19:11 2016 +0200 @@ -5,6 +5,7 @@ __glpk__.cc|$(GLPK_CPPFLAGS)|$(GLPK_LDFLAGS)|$(GLPK_LIBS) __init_fltk__.cc|$(FLTK_CPPFLAGS) $(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)|$(FLTK_LDFLAGS) $(FT2_LDFLAGS)|$(FLTK_LIBS) $(FT2_LIBS) $(OPENGL_LIBS) __init_gnuplot__.cc|$(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)|| +__ode15__.cc|$(SUNDIALS_XCPPFLAGS)|$(SUNDIALS_XLDFLAGS)|$(SUNDIALS_XLIBS) __osmesa_print__.cc|$(OSMESA_CPPFLAGS) $(FT2_CPPFLAGS)|$(OSMESA_LDFLAGS)|$(OSMESA_LIBS) __voronoi__.cc|$(QHULL_CPPFLAGS)|$(QHULL_LDFLAGS)|$(QHULL_LIBS) amd.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
--- a/scripts/ode/decic.m Tue Aug 23 02:31:51 2016 +0200 +++ b/scripts/ode/decic.m Tue Aug 23 03:19:11 2016 +0200 @@ -91,7 +91,7 @@ if (nargin < 6 || nargin > 7 || nargout > 3) print_usage (); endif - + #Check input if (! isa (odefun, "function_handle")) error ("Octave:invalid-input-arg",
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/ode15i.m Tue Aug 23 03:19:11 2016 +0200 @@ -0,0 +1,543 @@ +## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}) +## @deftypefnx {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}, @var{ode_opt}) +## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15i (@dots{}) +## @deftypefnx {} {@var{solution} =} ode15i (@dots{}) +## +## Solve a set of full-implicit Ordinary Differential Equations and +## Differential Algebraic Equations (DAEs) of index 1, with the variable-step, +## variable order BDF (Backward Differentiation Formula) method, which +## ranges from order 1 to 5. +## +## @var{fun} is a function handle, inline function, or string containing the +## name of the function that defines the ODE: @code{f(@var{t},@var{y},@var{yp})}. +## The function must accept three inputs where the first is time @var{t}, the +## second is a column vector of unknowns @var{y} and the third is a column +## vector of unknowns @var{yp}. +## +## @var{trange} specifies the time interval over which the ODE will be +## evaluated. Typically, it is a two-element vector specifying the initial and +## final times (@code{[tinit, tfinal]}). If there are more than two elements +## then the solution will also be evaluated at these intermediate time. +## +## @var{y0} and @var{yp0} contain the initial values for the unknowns @var{y} +## and @var{yp}. If they are row vectors then the solution @var{y} will be a +## matrix in which each column is the solution for the corresponding initial +## value in @var{y0} and @var{yp0}. +## +## @var{y0} and @var{yp0} must be consistent initial conditions, meaning that +## @code{f(@var{t},@var{y0},@var{yp0})=0} is satisfied. You can use function +## decic to compute consistent initial conditions, given initial guesses. +## +## The optional fifth argument @var{ode_opt} specifies non-default options to +## the ODE solver. It is a structure generated by @code{odeset}. +## +## The function typically returns two outputs. Variable @var{t} is a +## column vector and contains the times where the solution was found. The +## output @var{y} is a matrix in which each column refers to a different +## unknown of the problem and each row corresponds to a time in @var{t}. +## +## The output can also be returned as a structure @var{solution} which +## has field @var{x} containing the time where the solution was evaluated and +## field @var{y} containing the solution matrix for the times in @var{x}. +## Use @code{fieldnames (@var{solution})} to see the other fields and +## additional information returned. +## +## If using the @qcode{"Events"} option then three additional outputs may +## be returned. @var{te} holds the time when an Event function returned a +## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} +## contains an index indicating which Event function was triggered in the case +## of multiple Event functions. +## +## This function can be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, instead +## @var{y} is a matrix in which each column refers to a different unknown of +## the problem and the rows number is the same of @var{t} rows number so +## that each row of @var{y} contains the values of all unknowns at the time +## value contained in the corresponding row in @var{t}. +## +## Example: Solve the @nospell{Robetson}'s equations: +## +## @example +## @group +## function res = robertsidae(@var{t}, @var{y}, @var{yp}) +## res = [-(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3)); +## -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + +## 3e7*@var{y}(2)^2); +## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1]; +## endfunction +## [@var{t},@var{y}] = ode15i (@@robertsidae, [0 1e3], [1; 0; 0],[-1e-4; 1e-4; 0]); +## @end group +## @end example +## @seealso{decic, odeset, odeget} +## @end deftypefn + +function varargout = ode15i (fun, trange, y0, yp0, varargin) + + solver = 'ode15i'; + + if (nargin < 4) + print_usage (); + endif + + n = numel (y0); + + if (nargin > 4) + options = varargin{1}; + else + options = odeset (); + endif + + ## Check fun, trange, y0, yp0 + fun = check_default_input (fun, trange, solver, y0, yp0); + + if (! isempty (options.Jacobian)) + if (ischar (options.Jacobian)) + try + options.Jacobian = str2func (options.Jacobian); + catch + warning (lasterr); + end_try_catch + if (! isa (options.Jacobian, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + endif + endif + + if (! isempty (options.OutputFcn)) + if (ischar (options.OutputFcn)) + try + options.OutputFcn = str2func (options.OutputFcn); + catch + warning (lasterr); + end_try_catch + if (! isa (options.OutputFcn, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "OutputFcn"); + endif + endif + endif + + if (! isempty (options.Events)) + if (ischar (options.Events)) + try + options.Events = str2func (options.Events); + catch + warning (lasterr); + end_try_catch + if (! isa (options.Events, "function_handle") && ! ismatrix (options.Events)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Events"); + endif + endif + endif + + persistent defaults = []; + persistent classes = []; + persistent attributes = []; + + [defaults, classes, attributes] = odedefaults (n, trange(1), + trange(end)); + + defaults = rmfield (defaults, {"NonNegative", "Mass", ... + "MStateDependence", "MvPattern", ... + "MassSingular", "InitialSlope", "BDF"}); + classes = rmfield (classes, {"NonNegative", "Mass", ... + "MStateDependence", "MvPattern", ... + "MassSingular", "InitialSlope", "BDF"}); + attributes = rmfield (attributes, {"NonNegative", "Mass", ... + "MStateDependence", "MvPattern", ... + "MassSingular", "InitialSlope", "BDF"}); + + classes = odeset (classes, 'Vectorized', {}); + attributes = odeset (attributes, 'Jacobian', {}, 'Vectorized', {}); + + options = odemergeopts (options, defaults, classes, attributes, solver); + + ## Jacobian + options.havejac = false; + options.havejacsparse = false; + options.havejacfun = false; + + if (! isempty (options.Jacobian)) + options.havejac = true; + if (iscell (options.Jacobian)) + if (numel (options.Jacobian) == 2) + if (issparse (options.Jacobian{1}) && issparse (options.Jacobian{2})) ## Jac is sparse cell + options.havejacsparse = true; + endif + + if (any (size (options.Jacobian{1}) != [n n]) + || any (size (options.Jacobian{2}) != [n n]) + || ! isnumeric (options.Jacobian{1}) + || ! isnumeric (options.Jacobian{2}) + || ! isreal (options.Jacobian{1}) + || ! isreal (options.Jacobian{2})) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + + elseif (isa (options.Jacobian, "function_handle")) + options.havejacfun = true; + if (nargin (options.Jacobian) == 3) + [A, B] = options.Jacobian (trange(1), y0, yp0); + if (issparse (A) && issparse (B)) + options.havejacsparse = true; ## Jac is sparse fun + endif + + if (any (size (A) != [n n]) || any (size (B) != [n n]) + || ! isnumeric (A) || ! isnumeric (B) || ! isreal (A) + || ! isreal (B)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + endif + + ## Abstol and Reltol + + options.haveabstolvec = false; + + if (numel (options.AbsTol) != 1 && numel (options.AbsTol) != n) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "AbsTol"); + + elseif (numel (options.AbsTol) == n) + options.haveabstolvec = true; + endif + + ## Stats + options.havestats = false; + if (strcmp (options.Stats, "on")) + options.havestats = true; + endif + + ## Don't use Refine when the output is a structure + if (nargout == 1) + options.Refine = 1; + endif + + ## OutputFcn and OutputSel + if (isempty (options.OutputFcn) && nargout == 0) + options.OutputFcn = @odeplot; + options.haveoutputfunction = true; + else + options.haveoutputfunction = ! isempty (options.OutputFcn); + endif + + options.haveoutputselection = ! isempty (options.OutputSel); + if (options.haveoutputselection) + options.OutputSel = options.OutputSel - 1; + endif + + ## Events + options.haveeventfunction = ! isempty (options.Events); + + + [t, y, te, ye, ie] = __ode15__ (fun, trange, y0, yp0, options); + + + if (nargout == 2) + varargout{1} = t; + varargout{2} = y; + elseif (nargout == 1) + varargout{1}.x = t; # Time stamps are saved in field x + varargout{1}.y = y; # Results are saved in field y + varargout{1}.solver = solver; + if (options.haveeventfunction) + varargout{1}.xe = te; # Time info when an event occurred + varargout{1}.ye = ye; # Results when an event occurred + varargout{1}.ie = ie; # Index info which event occurred + endif + elseif (nargout == 5) + varargout = cell (1,5); + varargout{1} = t; + varargout{2} = y; + if (options.haveeventfunction) + varargout{3} = te; # Time info when an event occurred + varargout{4} = ye; # Results when an event occurred + varargout{5} = ie; # Index info which event occurred + endif + endif + +endfunction + +%!demo +%! +%! ##Solve Robertson's equations with ode15i +%! fun = @ (t, y, yp) [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); +%! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2); +%! y(1) + y(2) + y(3) - 1]; +%! +%! opt = odeset ('RelTol',1e-4, 'AbsTol', [1e-8, 1e-14, 1e-6]); +%! y0 = [1; 0; 0]; +%! yp0 = [-1e-4; 1e-4; 0]; +%! tspan = [0 4*logspace(-6, 6)]; +%! +%! [t, y] = ode15i (fun, tspan, y0, yp0, opt); +%! +%! y (:,2) = 1e4 * y (:, 2); +%! figure (2); +%! semilogx (t, y, 'o') +%! xlabel ('time'); +%! ylabel ('species concentration'); +%! title ('Robertson DAE problem with a Conservation Law'); +%! legend ('y1', 'y2', 'y3'); + +%!function res = rob (t, y, yp) +%! res =[-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); +%! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2); +%! y(1) + y(2) + y(3) - 1]; +%!endfunction +%! +%!function ref = fref() +%! ref = [100, 0.617234887614937, 0.000006153591397, 0.382758958793666]; +%!endfunction +%! +%!function ref2 = fref2() +%! ref2 = [4e6 0 0 1]; +%!endfunction +%! +%!function [DFDY, DFDYP] = jacfundense(t, y, yp) +%! DFDY = [-0.04, 1e4*y(3), 1e4*y(2); +%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); +%! 1, 1, 1]; +%! DFDYP = [-1, 0, 0; +%! 0, -1, 0; +%! 0, 0, 0]; +%!endfunction +%! +%!function [DFDY, DFDYP] = jacfunsparse(t, y, yp) +%! DFDY = sparse ([-0.04, 1e4*y(3), 1e4*y(2); +%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); +%! 1, 1, 1]); +%! DFDYP = sparse ([-1, 0, 0; +%! 0, -1, 0; +%! 0, 0, 0]); +%!endfunction +%!function [DFDY, DFDYP] = jacwrong(t, y, yp) +%! DFDY = [-0.04, 1e4*y(3); +%! 0.04, -1e4*y(3)-6e7*y(2)]; +%! DFDYP = [-1, 0; +%! 0, -1]; +%!endfunction +%!function [DFDY, DFDYP, A] = jacwrong2(t, y, yp) +%! DFDY = [-0.04, 1e4*y(3), 1e4*y(2); +%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); +%! 1, 1, 1]; +%! DFDYP = [-1, 0, 0; +%! 0, -1, 0; +%! 0, 0, 0]; +%! A = DFDY; +%!endfunction +%!function [val, isterminal, direction] = ff (t, y, yp) +%! isterminal = [0 1]; +%! if (t < 1e1) +%! val = [-1, -2]; +%! else +%! val = [1 3]; +%! endif +%! +%! direction = [1 0]; +%!endfunction + +%!test # anonymous function instead of real function +%! ref = [0.049787079136413]; +%! ff = @(t, u, udot) udot + 3 * u; +%! [t, y] = ode15i (ff, 0:1, 1, -3); +%! assert ([t(end), y(end)], [1, ref], 1e-3); +%!test # function passed as string +%! [t, y] = ode15i ('rob',[0 100 200], [1;0;0], [-1e-4;1e-4;0]); +%! assert ([t(2), y(2,:)], fref, 1e-3); +%!test # solve in intermidiate step +%! [t, y] = ode15i (@rob,[0 100 200], [1;0;0], [-1e-4;1e-4;0]); +%! assert ([t(2), y(2,:)], fref, 1e-3); +%!test # numel(trange) = 2 final value +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! assert ([t(end), y(end,:)], fref, 1e-5); +%!test # With empty options +%! opt = odeset(); +%! [t, y] = ode15i (@rob,[0 1e6 2e6 3e6 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref2, 1e-3); +%! opt = odeset(); +%!test # Without options +%! [t, y] = ode15i (@rob,[0 1e6 2e6 3e6 4e6], [1;0;0], [-1e-4;1e-4;0]); +%! assert ([t(end), y(end,:)], fref2, 1e-3); +%!test # InitialStep option +%! opt = odeset ("InitialStep", 1e-8); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(2)-t(1)], [1e-8], 1e-9); +%!test # MaxStep option +%! opt = odeset ("MaxStep", 1e-3); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! assert ([t(5)-t(4)], [1e-3], 1e-3); +%!test # AbsTol scalar option +%! opt = odeset ("AbsTol", 1e-8); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # AbsTol scalar and RelTol option +%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-6); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # AbsTol vector option +%! opt = odeset ("AbsTol", [1e-8, 1e-14,1e-6]); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # AbsTol vector and RelTol option +%! opt = odeset ("AbsTol", [1e-8, 1e-14,1e-6], "RelTol", 1e-6); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # RelTol option +%! opt = odeset ("RelTol", 1e-6); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # Jacobian fun dense +%! opt = odeset ("Jacobian", @jacfundense); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # Jacobian fun dense as string +%! opt = odeset ("Jacobian", 'jacfundense'); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # Jacobian fun sparse +%! opt = odeset ("Jacobian", @jacfunsparse, "AbsTol", 1e-7, "RelTol", 1e-7); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), y(end,:)], fref, 1e-3); +%!test # Solve in backward direction starting at t=100 +%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; +%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! [t2, y2] = ode15i (@rob,[100 0], Yref', YPref); +%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); +#%!test # Solve in backward direction with MaxStep option +%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; +%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; +%! opt = odeset ('MaxStep', 1e-2); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! [t2, y2] = ode15i (@rob,[100 0], Yref', YPref, opt); +%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); +%! assert ([t2(9)-t2(10)], [1e-2], 1e-2); +#%!test # Solve in backward direction starting with intermidiate step +%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; +%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! [t2, y2] = ode15i (@rob,[100 5 0], Yref', YPref); +%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); +%!test # Refine +%! opt = odeset ("Refine", 3); +%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); +%! [t2, y2] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([numel(t2)], numel(t)*3, 3); +%!test # Refine ignored if numel (trange) > 2 +%! opt = odeset ("Refine", 3); +%! [t, y] = ode15i (@rob,[0 10 100], [1;0;0], [-1e-4;1e-4;0]); +%! [t2, y2] = ode15i (@rob,[0 10 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([numel(t2)], numel(t)); +%!test # Events option add further elements in sol +%! opt = odeset ("Events", @ff); +%! sol = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert (isfield (sol, "ie")); +%! assert (sol.ie, [0;1]); +%! assert (isfield (sol, "xe")); +%! assert (isfield (sol, "ye")); +%! assert (sol.x(end), 10, 1); +%!test # Events option, five output arguments +%! opt = odeset ("Events", @ff); +%! [t, y, te, ye, ie] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); +%! assert ([t(end), te', ie'], [10, 10, 10, 0, 1], [1, 0.2, 0.2, 0, 0]); + +%!error # Jacobian fun wrong dimension +%! opt = odeset ("Jacobian", @jacwrong); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian cell dense wrong dimension +%! DFDY = [-0.04, 1; +%! 0.04, 1]; +%! DFDYP = [-1, 0, 0; +%! 0, -1, 0; +%! 0, 0, 0]; +%! opt = odeset ("Jacobian", {DFDY, DFDYP}); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian cell sparse wrong dimension +%! DFDY = sparse ([-0.04, 1; +%! 0.04, 1]); +%! DFDYP = sparse ([-1, 0, 0; +%! 0, -1, 0; +%! 0, 0, 0]); +%! opt = odeset ("Jacobian", {DFDY, DFDYP}); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian cell wrong number of matrices +%! A = [1 2 3; 4 5 6; 7 8 9]; +%! opt = odeset ("Jacobian", {A,A,A}); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian single matrix +%! opt = odeset ("Jacobian", [1 2 3; 4 5 6; 7 8 9]); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian single matrix wrong dimension +%! opt = odeset ("Jacobian", [1 2 3; 4 5 6]); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!error # Jacobian strange field +%! opt = odeset ("Jacobian", "foo"); +%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); +%!function ydot = fun (t, y, yp) +%! ydot = [y - yp]; +%!endfunction +%!error ode15i (); +%!error ode15i (1); +%!error ode15i (1, 1, 1); +%!error ode15i (1, 1, 1); +%!error ode15i (1, 1, 1, 1); +%!error ode15i (1, 1, 1, 1, 1); +%!error ode15i (1, 1, 1, 1, 1, 1); +%!error ode15i (@fun, 1, 1, 1); +%!error ode15i (@fun, [1, 1], 1, 1); +%!error ode15i (@fun, [1, 2], [1], [1, 2]); +%!error +%! opt = odeset ('RelTol', "foo"); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); +%!error +%! opt = odeset ('RelTol', [1, 2]); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); +%!error +%! opt = odeset ('RelTol', -2); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); +%!error +%! opt = odeset ('AbsTol', "foo"); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); +%!error +%! opt = odeset ('AbsTol', -1); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); +%!error +%! opt = odeset ('AbsTol', [1, 1, 1]); +%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/ode15s.m Tue Aug 23 03:19:11 2016 +0200 @@ -0,0 +1,686 @@ +## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}) +## @deftypefnx {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}, @var{ode_opt}) +## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15s (@dots{}) +## @deftypefnx {} {@var{solution} =} ode15s (@dots{}) +## +## Solve a set of stiff Ordinary Differential Equations and stiff semi-explicit +## Differential Algebraic Equations (DAEs) of index 1, with the variable-step, +## variable order BDF (Backward Differentiation Formula) method, which +## ranges from order 1 to 5. +## +## @var{fun} is a function handle, inline function, or string containing the +## name of the function that defines the ODE: @code{f(@var{t},@var{y})}. +## The function must accept two inputs where the first is time @var{t} and the +## second is a column vector of unknowns @var{y}. +## +## @var{trange} specifies the time interval over which the ODE will be +## evaluated. Typically, it is a two-element vector specifying the initial and +## final times (@code{[tinit, tfinal]}). If there are more than two elements +## then the solution will also be evaluated at these intermediate time. +## +## @var{init} contains the initial value for the unknowns. If it is a row +## vector then the solution @var{y} will be a matrix in which each column is +## the solution for the corresponding initial value in @var{init}. +## +## The optional fourth argument @var{ode_opt} specifies non-default options to +## the ODE solver. It is a structure generated by @code{odeset}. +## +## The function typically returns two outputs. Variable @var{t} is a +## column vector and contains the times where the solution was found. The +## output @var{y} is a matrix in which each column refers to a different +## unknown of the problem and each row corresponds to a time in @var{t}. +## +## The output can also be returned as a structure @var{solution} which +## has field @var{x} containing the time where the solution was evaluated and +## field @var{y} containing the solution matrix for the times in @var{x}. +## Use @code{fieldnames (@var{solution})} to see the other fields and +## additional information returned. +## +## If using the @qcode{"Events"} option then three additional outputs may +## be returned. @var{te} holds the time when an Event function returned a +## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} +## contains an index indicating which Event function was triggered in the case +## of multiple Event functions. +## +## This function can be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, instead +## @var{y} is a matrix in which each column refers to a different unknown of +## the problem and the rows number is the same of @var{t} rows number so +## that each row of @var{y} contains the values of all unknowns at the time +## value contained in the corresponding row in @var{t}. +## +## Example: Solve the @nospell{Robetson}'s equations: +## @example +## @group +## function res = robertsidae(@var{t}, @var{y}) +## res = [-0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3); +## 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3) - 3e7*@var{y}(2)^2; +## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1]; +## endfunction +## opt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none"); +## [@var{t},@var{y}] = ode15s (@@robertsidae, [0 1e3], [1; 0; 0], opt); +## @end group +## @end example +## @seealso{decic, odeset, odeget} +## @end deftypefn + +function varargout = ode15s (fun, trange, y0, varargin) + + solver = 'ode15s'; + + if (nargin < 3) + print_usage (); + endif + + ## Check fun, trange, y0, yp0 + fun = check_default_input (fun, trange, solver, y0); + + n = numel (y0); + + if (nargin > 3) + options = varargin{1}; + else + options = odeset (); + endif + + if (! isempty (options.Mass)) + if (ischar (options.Mass)) + try + options.Mass = str2func (options.Mass); + catch + warning (lasterr); + end_try_catch + if (! isa (options.Mass, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + endif + endif + + if (! isempty (options.Jacobian)) + if (ischar (options.Jacobian)) + try + options.Jacobian = str2func (options.Jacobian); + catch + warning (lasterr); + end_try_catch + if (! isa (options.Jacobian, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + endif + endif + + if (! isempty (options.OutputFcn)) + if (ischar (options.OutputFcn)) + try + options.OutputFcn = str2func (options.OutputFcn); + catch + warning (lasterr); + end_try_catch + if (! isa (options.OutputFcn, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "OutputFcn"); + endif + endif + endif + + if (! isempty (options.Events)) + if (ischar (options.Events)) + try + options.Events = str2func (options.Events); + catch + warning (lasterr); + end_try_catch + if (! isa (options.Events, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Events"); + endif + endif + endif + + persistent defaults = []; + persistent classes = []; + persistent attributes = []; + + [defaults, classes, attributes] = odedefaults (n, trange(1), + trange(end)); + + classes = odeset (classes, 'Vectorized', {}); + attributes = odeset (attributes, 'Jacobian', {}, 'Vectorized', {}); + + options = odemergeopts (options, defaults, classes, attributes, solver); + + ## Mass + + options.havemassfun = false; + options.havestatedep = false; + options.havetimedep = false; + options.havemasssparse = false; + + if (! isempty (options.Mass)) + if (isa (options.Mass, "function_handle")) + options.havemassfun = true; + if (nargin (options.Mass) == 2) + options.havestatedep = true; + M = options.Mass (trange(1), y0); + options.havemasssparse = issparse (M); + if (any (size (M) != [n n]) || ! isnumeric (M) || ! isreal (M)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + elseif (nargin (options.Mass) == 1) + options.havetimedep = true; + M = options.Mass (trange(1)); + options.havemasssparse = issparse (M); + if (any (size (M) != [n n]) || ! isnumeric (M) || ! isreal (M)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + elseif (ismatrix (options.Mass)) + options.havemasssparse = issparse (options.Mass); + if (any (size (options.Mass) != [n n]) || + ! isnumeric (options.Mass) || ! isreal (options.Mass)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Mass"); + endif + endif + + + ## Jacobian + options.havejac = false; + options.havejacsparse = false; + options.havejacfun = false; + + if (! isempty (options.Jacobian)) + options.havejac = true; + if (isa (options.Jacobian, "function_handle")) + options.havejacfun = true; + if (nargin (options.Jacobian) == 2) + [A] = options.Jacobian (trange(1), y0); + if (issparse (A)) + options.havejacsparse = true; ## Jac is sparse fun + endif + if (any (size (A) != [n n]) || ! isnumeric (A) || ! isreal (A)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + elseif (ismatrix (options.Jacobian)) + if (issparse (options.Jacobian)) + options.havejacsparse = true; ## Jac is sparse matrix + endif + if (! issquare (options.Jacobian)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + else + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "Jacobian"); + endif + endif + + + ## Derivative of M(t,y) for implicit problem not implemented yet + if (! isempty (options.Mass) && ! isempty (options.Jacobian)) + if (options.MStateDependence != "none" || options.havestatedep == true) + options.havejac = false; + options.Jacobian = []; + warning ("ode15s:mass_state_dependent_provided", + ["with MStateDependence != 'none' an internal", ... + " approximation of Jacobian Matrix will be used.", ... + " Set MStateDependence equal to 'none' if you want ", ... + " to provide a constant or time-dependent Jacobian"]); + endif + endif + + ## Use sparse methods only if all matrices are sparse + if (! options.havemasssparse) + options.havejacsparse = false; + endif + + ## If Mass or Jacobian is fun, then new Jacobian is fun + if (options.havejac) + if (options.havejacfun || options.havetimedep) + options.Jacobian = @ (t, y, yp) wrapjacfun (t, y, yp, + options.Jacobian, + options.Mass, + options.havetimedep, + options.havejacfun); + options.havejacfun = true; + else ## All matrices are constant + options.Jacobian = {[- options.Jacobian], [options.Mass]}; + + endif + endif + + + + ## Abstol and Reltol + + options.haveabstolvec = false; + + if (numel (options.AbsTol) != 1 && numel (options.AbsTol) != n) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "AbsTol"); + + elseif (numel (options.AbsTol) == n) + options.haveabstolvec = true; + endif + + ## Stats + options.havestats = false; + if (strcmp (options.Stats, "on")) + options.havestats = true; + endif + + ## Don't use Refine when the output is a structure + if (nargout == 1) + options.Refine = 1; + endif + + ## OutputFcn and OutputSel + if (isempty (options.OutputFcn) && nargout == 0) + options.OutputFcn = @odeplot; + options.haveoutputfunction = true; + else + options.haveoutputfunction = ! isempty (options.OutputFcn); + endif + + options.haveoutputselection = ! isempty (options.OutputSel); + if (options.haveoutputselection) + options.OutputSel = options.OutputSel - 1; + endif + + ## Events + options.haveeventfunction = ! isempty (options.Events); + + yp0 = options.InitialSlope; + + + [t, y, te, ye, ie] = __ode15__ (@ (t, y, yp) wrap (t, y, yp, options.Mass, + options.havetimedep, + options.havestatedep, + fun), + trange, y0, yp0, options); + + if (nargout == 2) + varargout{1} = t; + varargout{2} = y; + elseif (nargout == 1) + varargout{1}.x = t; # Time stamps are saved in field x + varargout{1}.y = y; # Results are saved in field y + varargout{1}.solver = solver; + if (options.haveeventfunction) + varargout{1}.xe = te; # Time info when an event occurred + varargout{1}.ye = ye; # Results when an event occurred + varargout{1}.ie = ie; # Index info which event occurred + endif + elseif (nargout == 5) + varargout = cell (1,5); + varargout{1} = t; + varargout{2} = y; + if (options.haveeventfunction) + varargout{3} = te; # Time info when an event occurred + varargout{4} = ye; # Results when an event occurred + varargout{5} = ie; # Index info which event occurred + endif + endif + +endfunction + +function res = wrap (t, y, yp, Mass, havetimedep, havestatedep, fun) + + if (! isempty (Mass) && havestatedep) + res = Mass (t, y) * yp - fun (t, y); + elseif (! isempty (Mass) && havetimedep) + res = Mass (t) * yp - fun (t, y); + elseif (! isempty (Mass)) + res = Mass * yp - fun (t, y); + else + res = yp - fun (t, y); + endif + +endfunction + +function [jac, jact] = wrapjacfun (t, y, yp, Jac, Mass, + havetimedep, havejacfun) + if (havejacfun) + jac = - Jac (t, y); + else + jac = - Jac; + endif + + if (! isempty (Mass) && havetimedep) + jact = Mass (t); + elseif (! isempty (Mass)) + jact = Mass; + else + jact = speye (numel (y)); + endif + +endfunction + + +%!demo +%! +%! ##Solve Robertson's equations with ode15s +%! fun = @ (t, y) [-0.04*y(1) + 1e4*y(2).*y(3); +%! 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2; +%! y(1) + y(2) + y(3) - 1 ]; +%! +%! y0 = [1; 0; 0]; +%! tspan = [0 4*logspace(-6, 6)]; +%! M = [1 0 0; 0 1 0; 0 0 0]; +%! +%! options = odeset ('RelTol', 1e-4, 'AbsTol', [1e-6 1e-10 1e-6], +%! 'MStateDependence', 'none', 'Mass', M); +%! +%! [t, y] = ode15s (fun, tspan, y0, options); +%! +%! y (:,2) = 1e4 * y (:,2); +%! figure (2); +%! semilogx (t, y, 'o') +%! xlabel ('time'); +%! ylabel ('species concentration'); +%! title ('Robertson DAE problem with a Conservation Law'); +%! legend ('y1', 'y2', 'y3'); + +%!function ydot = fpol (t, y) # The Van der Pol +%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%!endfunction +%!function ref = fref () # The computed reference sol +%! ref = [0.32331666704577, -1.83297456798624]; +%!endfunction +%!function jac = fjac (t, y) # its Jacobian +%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; +%!endfunction +%!function jac = fjcc (t, y) # sparse type +%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); +%!endfunction +%!function mas = fmas (t, y) +%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests +%!endfunction +%!function mas = fmsa (t, y) +%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix +%!endfunction +%!function res = rob (t, y) +%! res = [-0.04*y(1) + 1e4*y(2).*y(3); +%! 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2; +%! y(1) + y(2) + y(3) - 1 ]; +%!endfunction +%! +%!function refrob = frefrob() +%! refrob = [100, 0.617234887614937, 0.000006153591397, 0.382758958793666]; +%!endfunction +%!function [val, isterminal, direction] = feve (t, y) +%! isterminal = [0 1]; +%! if (t < 1e1) +%! val = [-1, -2]; +%! else +%! val = [1 3]; +%! endif +%! +%! direction = [1 0]; +%!endfunction +%!function masrob = massdensefunstate (t, y) +%! masrob = [1 0 0; 0 1 0; 0 0 0]; +%!endfunction +%!function masrob = masssparsefunstate (t, y) +%! masrob = sparse([1 0 0; 0 1 0; 0 0 0]); +%!endfunction +%!function masrob = massdensefuntime (t) +%! masrob = [1 0 0; 0 1 0; 0 0 0]; +%!endfunction +%!function masrob = masssparsefuntime (t) +%! masrob = sparse([1 0 0; 0 1 0; 0 0 0]); +%!endfunction +%!function jac = jacfundense (t, y) +%! jac = [-0.04, 1e4*y(3), 1e4*y(2); +%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); +%! 1, 1, 1]; +%!endfunction +%!function jac = jacfunsparse (t, y) +%! jac = sparse([-0.04, 1e4*y(3), 1e4*y(2); +%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); +%! 1, 1, 1]); +%!endfunction + + + +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", [1 0 0; 0 1 0; 0 0 0]); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0])); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @massdensefunstate); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @masssparsefunstate); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", 'massdensefuntime'); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", [1 0 0; 0 1 0; 0 0 0], +%! "Jacobian", 'jacfundense'); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0]), +%! "Jacobian", @jacfundense); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @massdensefunstate, +%! "Jacobian", @jacfundense); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @masssparsefunstate, +%! "Jacobian", @jacfundense); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @massdensefuntime, +%! "Jacobian", @jacfundense); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%! opt = odeset ("MStateDependence", "none", +%! "Mass", 'masssparsefuntime', +%! "Jacobian", 'jacfundense'); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", [1 0 0; 0 1 0; 0 0 0], +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0]), +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @massdensefunstate, +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @masssparsefunstate, +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @massdensefuntime, +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test +%! opt = odeset ("MStateDependence", "none", +%! "Mass", @masssparsefuntime, +%! "Jacobian", @jacfunsparse); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), y(end,:)], frefrob, 1e-3); +%!test # two output arguments +%! [t, y] = ode15s (@fpol, [0 2], [2 0]); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); +%!test # anonymous function instead of real function +%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! [t, y] = ode15s (fvdb, [0 2], [2 0]); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); +%!test # Solve another anonymous function below zero +%! ref = [0, 14.77810590694212]; +%! [t, y] = ode15s (@(t,y) y, [-2 0], 2); +%! assert ([t(end), y(end,:)], ref, 5e-2); +%!test # InitialStep option +%! opt = odeset ("InitialStep", 1e-8); +%! [t, y] = ode15s (@fpol, [0 0.2], [2 0], opt); +%! assert ([t(2)-t(1)], [1e-8], 1e-9); +%!test # MaxStep option +%! opt = odeset ("MaxStep", 1e-3); +%! sol = ode15s (@fpol, [0 0.2], [2 0], opt); +%! assert ([sol.x(5)-sol.x(4)], [1e-3], 1e-3); +%!test # Solve in backward direction starting at t=0 +%! ref = [-1.205364552835178, 0.951542399860817]; +%! sol = ode15s (@fpol, [0 -2], [2 0]); +%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 5e-3); +%!test # Solve in backward direction starting at t=2 +%! ref = [-1.205364552835178, 0.951542399860817]; +%! sol = ode15s (@fpol, [2 0 -2], fref); +%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 3e-2); +%!test # Solve another anonymous function in backward direction +%! ref = [-1, 0.367879437558975]; +%! sol = ode15s (@(t,y) y, [0 -1], 1); +%! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); +%!test # Solve another anonymous function below zero +%! ref = [0, 14.77810590694212]; +%! sol = ode15s (@(t,y) y, [-2 0], 2); +%! assert ([sol.x(end), sol.y(end,:)], ref, 5e-2); +%!test # Solve in backward direction starting at t=0 with MaxStep option +%! ref = [-1.205364552835178, 0.951542399860817]; +%! opt = odeset ("MaxStep", 1e-3); +%! sol = ode15s (@fpol, [0 -2], [2 0], opt); +%! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); +%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3); +%!test # AbsTol option +%! opt = odeset ("AbsTol", 1e-5); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 4e-3); +%!test # AbsTol and RelTol option +%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%!test # RelTol option -- higher accuracy +%! opt = odeset ("RelTol", 1e-8); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-4); +%!test # Mass option as function +%! opt = odeset ("Mass", @fmas, "MStateDependence", "none"); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); +%!test # Mass option as matrix +%! opt = odeset ("Mass", eye (2,2), "MStateDependence", "none"); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); +%!test # Mass option as sparse matrix +%! opt = odeset ("Mass", sparse (eye (2,2)), "MStateDependence", "none"); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); +%!test # Mass option as function and sparse matrix +%! opt = odeset ("Mass", 'fmsa', "MStateDependence", "none"); +%! sol = ode15s (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); +%!test # Refine +%! opt2 = odeset ("Refine", 3, "Mass", @massdensefunstate, +%! "MStateDependence", "none"); +%! opt1 = odeset ("Mass", @massdensefunstate, "MStateDependence", "none"); +%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt1); +%! [t2, y2] = ode15s (@rob,[0 100], [1;0;0], opt2); +%! assert ([numel(t2)], numel(t)*3, 3); +%!test # Refine ignored if numel (trange) > 2 +%! opt2 = odeset ("Refine", 3, "Mass", 'massdensefunstate', +%! "MStateDependence", "none"); +%! opt1 = odeset ("Mass", @massdensefunstate, "MStateDependence", "none"); +%! [t, y] = ode15s ('rob',[0 10 100], [1;0;0], opt1); +%! [t2, y2] = ode15s ('rob',[0 10 100], [1;0;0], opt2); +%! assert ([numel(t2)], numel(t)); +%!test # Events option add further elements in sol +%! opt = odeset ("Events", @feve, "Mass", @massdensefunstate, +%! "MStateDependence", "none"); +%! sol = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert (isfield (sol, "ie")); +%! assert (sol.ie, [0;1]); +%! assert (isfield (sol, "xe")); +%! assert (isfield (sol, "ye")); +%! assert (sol.x(end), 10, 1); +%!test # Events option, five output arguments +%! opt = odeset ("Events", @feve, "Mass", @massdensefunstate, +%! "MStateDependence", "none"); +%! [t, y, te, ye, ie] = ode15s (@rob,[0 100], [1;0;0], opt); +%! assert ([t(end), te', ie'], [10, 10, 10, 0, 1], [1, 0.5, 0.5, 0, 0]); + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/check_default_input.m Tue Aug 23 03:19:11 2016 +0200 @@ -0,0 +1,65 @@ +## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +function [fun] = check_default_input (fun, trange, solver, varargin); + + ## Check fun + validateattributes (fun, {"function_handle", "char"}, {}, solver, "fun"); + + if (! (nargin (fun) == nargin - 2)) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "fun"); + endif + + if (ischar (fun)) + try + fun = str2func (fun); + catch + warning (lasterr); + end_try_catch + endif + if (! isa (fun, "function_handle")) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "fun"); + endif + + ## Check trange + validateattributes (trange, {"float"}, {"vector", "real"}, solver, "trange"); + + if (numel (trange) < 2) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "trange"); + elseif (! ((all (diff (trange) > 0)) || all (diff (-trange) > 0))) + error ("Octave:invalid-input-arg", + [solver ": invalid value assigned to field '%s'"], "trange"); + endif + + ## Check y0 and yp0 + y0 = varargin{1}; + + if (nargin == 5) + yp0 = varargin{2}; + n = numel (feval (fun, trange(1), y0, yp0)); + validateattributes (yp0, {"float"}, {"numel", n}, solver, "yp0"); + else + n = numel (feval (fun, trange (1), y0)); + endif + + validateattributes (y0, {"float"}, {"numel", n}, solver, "y0"); + +endfunction