octave_packages/ocs-0.1.3/tst/tst_odepkg.m

   1 ## Copyright (C) 2006,2007,2008  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}, [@var{tout}]] =}  tst_odepkg @
  24 ## (@var{cirstruct},@var{x},@var{t},@var{tol},@var{maxit},@
  25 ## @var{pltvars},@var{solver},@var{odestruct},@var{verbosity});
  26 ##
  27 ## Perform a transient simulation of the system described by
  28 ## @var{cirstruct} over the time interval @var{t} using the @code{odepkg} DAE
  29 ## solver specified in @var{solver}.
  30 ##
  31 ## Pssible values for @var{solver} are @code{ode2r}, @code{ode5r},
  32 ## @code{oders} or @code{odesx}.
  33 ##
  34 ## The initial value for the state vector is computed by solving a
  35 ## steady state problem at @var{t}(1), with starting guess @var{x}.
  36 ##
  37 ## @var{tol} and @var{maxit} are parameters passed to @code{nls_newton_raphson}.
  38 ##
  39 ## If one output is requested @var{out} will contain the value of the state vector
  40 ## at each point of @var{t}.
  41 ##
  42 ## If two outputs are requested @var{out} will contain the value of the state vector
  43 ## at each point of @var{tout}.
  44 ##
  45 ## Extra options for options for the solver can be passed to the solver
  46 ## via @var{odestruct}.
  47 ##
  48 ## The optional parameter  @var{verbosity} controls the amount of
  49 ## output produced:
  50 ##
  51 ## @itemize @minus
  52 ## @item if verbosity(1) != 0, information on the progress
  53 ## of the algorithm are output at runtime
  54 ## @item if verbosity(2) != 0, the plot of the variables whose names
  55 ## are listed in @var{pltvars} is
  56 ## produced after the computation
  57 ## @end itemize
  58 ##
  59 ## @seealso{tst_backward_euler,tst_theta_method,tst_daspk,nls_newton_raphson,odepkg,odeset,@
  60 ## ode2r,ode5r,oders,odesx}
  61 ##
  62 ## @end deftypefn
  63
  64 function [out, tout] = tst_odepkg (outstruct,x,t,tol,maxit,\
  65                                    pltvars,solver,verbosity,odestruct)
  66
  67   ## Check input
  68   ## FIXME: add input check!
  69   if ((nargin < 7) || (nargin > 9))
  70     error("tst_odepkg:wrong number of input parameters");
  71   endif
  72
  73   if ~exist("verbosity")
  74     verbosity = [0,0];
  75   elseif length(verbosity)<2
  76     verbosity(2) =0;
  77   endif
  78
  79   if(verbosity(1))
  80     fprintf(1,"initial value:\n");
  81   endif
  82
  83   if ~exist("odestruct")
  84     odestruct = odeset();
  85   endif
  86
  87   [A0,B,C] = asm_initialize_system(outstruct,x);
  88
  89   [out(:,1),ii] = nls_stationary(outstruct,x,tol,maxit);
  90
  91   JAC = @(t, x) TSTODEPKGFUNJAC(outstruct, x, A0, B, t);
  92   RES = @(t, x) TSTODEPKGFUNRES(outstruct, x, A0, B, C, t);
  93   MASS= @(t, x) TSTODEPKGFUNMASS(outstruct, x, A0, t);
  94
  95   odestruct = odeset(odestruct,"Jacobian", JAC);
  96   odestruct = odeset(odestruct,"Mass",MASS(0,x));
  97   odestruct = odeset(odestruct,"RelTol", 1e-6,"AbsTol",1e6*eps,
  98                      "MaxStep", max(diff(t)),"InitialStep",(diff(t))(1));
  99
 100   if verbosity(2)
 101     odestruct = odeset(odestruct, "OutputFcn",
 102                        @(t, y, flag)  plotfun(t, y, flag, outstruct, pltvars) );
 103   endif
 104
 105   [tout, out] = feval( solver, RES, t([1 end]), x, odestruct);
 106   if (nargout < 2)
 107     out = interp1(tout, out, t).';
 108   endif
 109
 110 endfunction
 111
 112 function [varargout] = plotfun (vt, vy, vflag, outstruct, pltvars)
 113   ## this function is a modified version of odeplot distributed
 114   ## with odepkg (c) Thomas Treichl
 115
 116   %# No input argument check is done for a higher processing speed
 117   persistent vfigure; persistent vtold;
 118   persistent vyold;   persistent vcounter;
 119
 120   if (strcmp (vflag, "init"))
 121     %# Nothing to return, vt is either the time slot [tstart tstop]
 122     %# or [t0, t1, ..., tn], vy is the inital value vector 'vinit'
 123     vfigure  = figure; vtold = vt(1,1); vyold = vy(:,1);
 124     vcounter = 1;
 125
 126   elseif (isempty (vflag))
 127     %# Return something in varargout{1}, either false for 'not stopping
 128     %# the integration' or true for 'stopping the integration'
 129     vcounter = vcounter + 1; figure (vfigure);
 130     vtold(vcounter,1) = vt(1,1);
 131     vyold(:,vcounter) = vy(:,1);
 132     utl_plot_by_name(vtold, vyold, outstruct, pltvars); drawnow;
 133     varargout{1} = false;
 134
 135   elseif (strcmp (vflag, "done"))
 136     %# Cleanup has to be done, clear the persistent variables because
 137     %# we don't need them anymore
 138     clear ("vfigure","vtold","vyold","vcounter");
 139
 140   endif
 141
 142 endfunction
 143 ## Jacobian for transient problems
 144 function lhs = TSTODEPKGFUNJAC(outstruct,x,A0,B,t)
 145
 146   [A1,Jac,res] = asm_build_system(outstruct,x,t);
 147   lhs = ( B + Jac);
 148
 149 endfunction
 150
 151 function lhs = TSTODEPKGFUNMASS(outstruct,x,A0,t)
 152
 153   [A1,Jac,res] = asm_build_system(outstruct,x,t);
 154   lhs = -(A0+A1);
 155
 156 endfunction
 157
 158 ## Residual for transient problems
 159 function rhs = TSTODEPKGFUNRES(outstruct,x,A0,B,C,t)
 160
 161   [A1,Jac,res] = asm_build_system(outstruct,x,t);
 162   rhs =  B*x + C + res;
 163
 164 endfunction