octave_packages/ocs-0.1.3/tst/tst_daspk.m

   1 ## Copyright (C) 2006,2007,2008,2011  Carlo de Falco
   2 ##
   3 ## This file is part of:
   4 ## OCS - A Circuit Simulator for Octave
   5 ##
   6 ## OCS is free software; you can redistribute it and/or modify
   7 ## it under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation.
   9 ##
  10 ## This program is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program (see the file LICENSE); if not,
  17 ## see <http://www.gnu.org/licenses/>.
  18 ##
  19 ## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
  20
  21 ## -*- texinfo -*-
  22 ##
  23 ## @deftypefn{Function File} {[@var{out}] =}  tst_daspk @
  24 ## (@var{cirstruct},@var{x},@var{t},@var{tol},@var{maxit},@
  25 ## @var{pltvars},@var{verbosity},@var{daskopts},@var{dae_fun});
  26 ##
  27 ## Perform a transient simulation of the system described by
  28 ## @var{cirstruct} over the time interval @var{t} using @code{daspk}.
  29 ##
  30 ## The initial value for the state vector is computed by solving a
  31 ## steady state problem at @var{t}(1), with starting guess @var{x}.
  32 ##
  33 ## @var{tol} and @var{maxit} are parameters passed to
  34 ## @code{nls_newton_raphson}.
  35 ##
  36 ## The output @var{out} will contain the value of the state vector
  37 ## at each point of @var{t}.
  38 ##
  39 ## Extra options for @code{daspk} can be passed as name/value pairs in
  40 ## the cellarray @var{daskopts}.
  41 ##
  42 ## The optional parameter  @var{verbosity} controls the amount of
  43 ## output produced:
  44 ##
  45 ## @itemize @minus
  46 ## @item if verbosity(1) != 0, information on the progress
  47 ## of the algorithm are output at runtime
  48 ## @item if verbosity(2) != 0, the plot of the variables whose names
  49 ## are listed in @var{pltvars} is
  50 ## produced after the computation
  51 ## @end itemize
  52 ##
  53 ## For special purposes one may need to pass modified jacobian and
  54 ## residual functions. This can be done via the cell array of function
  55 ## handles @var{dae_fun}.
  56 ##
  57 ## Such functions should have the same input and output
  58 ## parameter list as the default sub-functions
  59 ## TSTBWEFUNJAC0,TSTBWEFUNRES0, TSTBWEFUNJAC,TSTBWEFUNRES.
  60 ##
  61 ## @seealso{tst_backward_euler,tst_odepkg,tst_theta_method,nls_newton_raphson,daspk}
  62 ##
  63 ## @end deftypefn
  64
  65 function [out] = tst_daspk (outstruct, x, t, tol, maxit, pltvars, verbosity, daspkopts, dae_fun)
  66
  67   ## FIXME: add input check!
  68   if ((nargin < 6) || (nargin > 9))
  69     error ("tst_daspk: wrong number of input parameters.");
  70   endif
  71
  72   if (! exist ("verbosity"))
  73     verbosity = [0, 0];
  74   elseif (length (verbosity) < 2)
  75     verbosity(2) = 0;
  76   endif
  77
  78   if (verbosity(1))
  79     fprintf (1, "initial value:\n");
  80   endif
  81
  82   daspk_options ("print initial condition info",1);
  83   daspk_options ("maximum order",2);
  84   daspk_options ("initial step size",t(2)-t(1));
  85   daspk_options ("relative tolerance",1e-3);
  86
  87   if (nargin > 8)
  88     for ii = 1:2:length (daspkopts)
  89       daspk_options (daspkopts{ii}, daspkopts{ii+1});
  90     endfor
  91   endif
  92
  93
  94   [A0, B, C] = asm_initialize_system (outstruct, x);
  95
  96   if (nargin > 9)
  97     JAC = @(x) dae_fun{1} (outstruct, x, t(1), B);
  98     RES = @(x) dae_fun{2} (outstruct, x, t(1), B, C);
  99     [x, ii, resnrm] = nls_newton_raphson (x, RES, JAC, tol, maxit, verbosity);
 100   else
 101     [out(:,1), ii] = nls_stationary (outstruct, x, tol, maxit);
 102   endif
 103
 104   if (nargin > 9)
 105     JAC = @(x) dae_fun{3} (outstruct, x, t(1), B);
 106     RES = @(x) dae_fun{4} (outstruct, x, t(1), B, C);
 107   else
 108     JAC = @(x,xdot,t,c) TSTDASPKFUNJAC (outstruct, x, xdot, A0, B, t, c);
 109     RES = @(x,xdot,t)   TSTDASPKFUNRES (outstruct, x, xdot, A0, B, C, t);
 110   endif
 111
 112   [out, xdot, istate, msg] = daspk ({RES,JAC}, out(:,1), zeros (size (x)), t);
 113
 114   out = out.';
 115   if (verbosity(2))
 116     utl_plot_by_name (t, out, outstruct, pltvars)
 117   endif
 118
 119 endfunction
 120
 121 ## Jacobian for transient problem
 122 function lhs = TSTDASPKFUNJAC (outstruct, x, xdot, A0, B, t, c)
 123
 124   [A1, Jac, res] = asm_build_system (outstruct, x, t);
 125   lhs = (c*(A0+A1) + B + Jac);
 126
 127 endfunction
 128
 129 ## Residual for transient problem
 130 function rhs = TSTDASPKFUNRES (outstruct, x, xdot, A0, B, C, t)
 131
 132   [A1, Jac, res] = asm_build_system (outstruct, x, t);
 133   rhs = (A0+A1)*xdot + B*x + C + res;
 134
 135 endfunction