octave_packages/ocs-0.1.3/tst/tst_backward_euler.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},@var{niter}] =}  @
  24 ## tst_backward_euler(@var{cirstruct},@var{x},@var{t},@var{tol},@
  25 ## @var{maxit},@var{pltvars},@var{verbosity},@var{dae_fun})
  26 ##
  27 ## Perform a transient simulation of the system described by
  28 ## @var{cirstruct} over the time interval @var{t} using the backward
  29 ## Euler algorithm.
  30 ##
  31 ## The initial value for the state vector is computed by solving a
  32 ## steady state problem at @var{t}(1), with starting guess @var{x}.
  33 ##
  34 ## @var{tol} and @var{maxit} are parameters passed to
  35 ## @code{nls_newton_raphson}.
  36 ##
  37 ## The output @var{out} will contain the value of the state vector at
  38 ## each point of @var{t}.
  39 ##
  40 ## The optional parameter @var{verbosity} controls the amount of
  41 ## output produced:
  42 ##
  43 ## @itemize @minus
  44 ## @item if verbosity(1) != 0, information on the progress
  45 ## of the algorithm are output at runtime
  46 ## @item if verbosity(2) != 0, the plot of the variables whose names
  47 ## are listed in @var{pltvars} is
  48 ## produced after the computation
  49 ## @end itemize
  50 ##
  51 ## For special purposes one may need to pass modified jacobian and
  52 ## residual functions. This can be done via the cell array of function
  53 ## handles @var{dae_fun}.
  54 ##
  55 ## Such functions should have the same input and output
  56 ## parameter list as the default sub-functions
  57 ## TSTBWEFUNJAC0,TSTBWEFUNRES0, TSTBWEFUNJAC,TSTBWEFUNRES.
  58 ##
  59 ## The optional output @var{niter} returns the number of Newton iterations
  60 ## needed to reach convergence.
  61 ##
  62 ## @seealso{tst_daspk,tst_theta_method,tst_odepkg,nls_newton_raphson}
  63 ##
  64 ## @end deftypefn
  65
  66 function [out, varargout] = tst_backward_euler (outstruct, x, t, tol, maxit, pltvars, verbosity, dae_fun)
  67
  68   ## Check input
  69   ## FIXME: add input check!
  70   if ((nargin < 6) || (nargin > 8))
  71     error ("tst_backward_euler: wrong number of input parameters.");
  72   endif
  73
  74   if ~exist ("verbosity")
  75     verbosity = [0,0];
  76   elseif (length (verbosity) < 2)
  77     verbosity(2) = 0;
  78   endif
  79
  80   out      = zeros (rows (x), columns (t));
  81   out(:,1) = x;
  82
  83   if nargout > 1
  84     niter = zeros (length(t),1);
  85   endif
  86
  87   if (verbosity(1))
  88     fprintf (1, "Initial value.\n");
  89   endif
  90
  91
  92   [A0, B, C] = asm_initialize_system (outstruct, x);
  93
  94   if (nargin > 8)
  95     JAC = @(x) dae_fun{1} (outstruct,x,t(1),B);
  96     RES = @(x) dae_fun{2} (outstruct,x,t(1),B,C);
  97     [out(:,1), ii, resnrm] = nls_newton_raphson (x, RES, JAC, tol, maxit, verbosity(1));
  98   else
  99     [out(:,1),ii] = nls_stationary (outstruct, x, tol, maxit);
 100   endif
 101
 102   if (nargout > 1)
 103     niter(1) = ii;
 104   endif
 105
 106   for it = 2:length (t)
 107
 108     if (verbosity(1))
 109       fprintf (1,"Timestep #%d.\n",it);
 110     endif
 111
 112     if nargin > 8
 113       JAC = @(x) dae_fun{3} (outstruct, x, t(it-1), t(it), A0, B);
 114       RES = @(x) dae_fun{4} (outstruct, x, out(:,it-1), t(it-1), t(it), A0, B, C);
 115     else
 116       JAC = @(x,A1,Jac,res) TSTBWEFUNJAC1(outstruct, x, t(it-1),
 117                                           t(it), A0, B, A1, Jac, res);
 118       RES = @(x,A1,Jac,res) TSTBWEFUNRES1(outstruct, x, out(:,it-1),
 119                                           t(it-1), t(it), A0, B, C,
 120                                           A1, Jac, res);
 121       UPDT = @(x) TSTBWEFUNUP1 (outstruct, x, t(it));
 122     endif
 123
 124     [out(:,it),ii,resnrm] = nls_newton_raphson (out(:,it-1), RES, JAC, tol, maxit, verbosity(1), UPDT);
 125
 126     if nargout > 1
 127       niter(it) = ii;
 128     endif
 129
 130     if (verbosity(2))
 131      utl_plot_by_name (t(1:it), out(:,1:it), outstruct, pltvars);
 132      drawnow ();
 133     endif
 134
 135     ## Stop at runtime
 136     ## FIXME: maintain this part?
 137     if exist("~/.stop_ocs","file")
 138       printf("stopping at timestep %d\n",it);
 139       unix("rm ~/.stop_ocs");
 140       break
 141     end
 142
 143   endfor
 144
 145   if nargout > 1
 146     varargout{1} = niter;
 147   endif
 148
 149 endfunction
 150
 151 ## Jacobian for transient problem
 152 function lhs = TSTBWEFUNJAC1 (outstruct, x, t0, t1, A0, B, A1, Jac, res)
 153
 154   DT = t1-t0;
 155   if (nargin < 9)
 156     [A1, Jac, res] = asm_build_system (outstruct, x, t1);
 157   endif
 158   lhs = ((A0+A1)/DT + B + Jac);
 159
 160 endfunction
 161
 162 ## Residual for transient problem
 163 function rhs = TSTBWEFUNRES1 (outstruct, x, xold, t0, t1, A0, B, C, A1, Jac, res)
 164
 165   DT = t1-t0;
 166   if ( nargin < 11 )
 167     [A1, Jac, res] = asm_build_system (outstruct, x, t1);
 168   endif
 169   rhs = (res + C + B*x + (A0+A1)*(x-xold)/DT);
 170
 171 endfunction
 172
 173 ## Update for transient problem
 174 function update = TSTBWEFUNUP1 (outstruct, x, t1)
 175
 176   [A1, Jac, res] = asm_build_system (outstruct, x, t1);
 177   update = {A1, Jac, res};
 178
 179 endfunction