octave_packages/ocs-0.1.3/tst/tst_theta_method.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{niter}] =}  tst_theta_method @
  24 ## (@var{cirstruct},@var{x},@var{t},@var{tol},@
  25 ## @var{maxit},@var{theta},@var{pltvars},@
  26 ## @var{verbosity});
  27 ##
  28 ## Perform a transient simulation of the system described by
  29 ## @var{cirstruct} over the time interval @var{t} using the
  30 ## theta-method with parameter @var{theta}.
  31 ##
  32 ## The initial value for the state vector is computed by solving a
  33 ## steady state problem at @var{t}(1), with starting guess @var{x}.
  34 ##
  35 ## @var{tol} and @var{maxit} are parameters passed to
  36 ## @code{nls_newton_raphson}.
  37 ##
  38 ## The output @var{out} will contain the value of the state vector at
  39 ## each point of @var{t}.
  40 ##
  41 ## The optional parameter @var{verbosity} controls the amount of
  42 ## output produced:
  43 ##
  44 ## @itemize @minus
  45 ## @item if verbosity(1) != 0, information on the progress
  46 ## of the algorithm are output at runtime
  47 ## @item if verbosity(2) != 0, the plot of the variables whose names
  48 ## are listed in @var{pltvars} is
  49 ## produced after the computation
  50 ## @end itemize
  51 ##
  52 ## The optional output @var{niter} returns the number of Newton iterations
  53 ## needed to reach convergence.
  54 ##
  55 ## @seealso{tst_backward_euler,tst_daspk,tst_odepkg,nls_newton_raphson}
  56 ##
  57 ## @end deftypefn
  58
  59 function [out, varargout] = tst_theta_method(outstruct,x,t,tol,maxit,\
  60                                              theta,pltvars,verbosity)
  61
  62   ## Check input
  63   ## FIXME: add input check!
  64   if ((nargin < 7) || (nargin > 8))
  65     error("tst_theta_method: wrong number of input parameters.");
  66   endif
  67
  68   if ~exist("verbosity")
  69     verbosity = [0,0];
  70   elseif length(verbosity)<2
  71     verbosity(2) = 0;
  72   endif
  73
  74   out = zeros(rows(x),columns(t));
  75   out(:,1) = x;
  76
  77   if nargout > 1
  78     niter = zeros(length(t),1);
  79   endif
  80
  81   if (verbosity(1))
  82     fprintf(1,"Initial value.\n");
  83   endif
  84
  85
  86   [A0,B,C] = asm_initialize_system(outstruct,x);
  87
  88   [out(:,1),ii] = nls_stationary(outstruct,x,tol,maxit);
  89
  90   if nargout > 1
  91     niter(1) = ii;
  92   endif
  93
  94   for it=2:length(t)
  95
  96     if(verbosity)
  97       fprintf(1,"Timestep #%d.\n",it);
  98     endif
  99
 100
 101     [A1old,Jacold,resold] = asm_build_system(outstruct, out(:,it-1), t(it-1));
 102
 103     JAC = @(x,A1,Jac,res) TSTTHETAFUNJAC1(outstruct, x, t(it-1),
 104                                           t(it), A0, B, theta,
 105                                           A1, Jac, res);
 106     RES = @(x,A1,Jac,res) TSTTHETAFUNRES1(outstruct, x, out(:,it-1),
 107                                           t(it-1), t(it), A0, B, C,
 108                                           resold, theta, A1, Jac, res);
 109     UPDT = @(x) TSTTHETAFUNUP1 (outstruct, x, t(it));
 110
 111     [out(:,it),ii,resnrm] = nls_newton_raphson(out(:,it-1),RES,JAC,\
 112                                                tol, maxit,verbosity(1),\
 113                                                UPDT);
 114
 115     if nargout > 1
 116       niter(it) = ii;
 117     endif
 118
 119     if (verbosity(2))
 120       utl_plot_by_name(t(1:it),out(:,1:it),outstruct,pltvars);
 121       pause(.1);
 122     endif
 123
 124     if exist("~/.stop_ocs","file")
 125       break
 126     end
 127   endfor
 128
 129   if nargout > 1
 130     varargout{1} = niter;
 131   endif
 132
 133 endfunction
 134
 135 ## Jacobian for transient problem
 136 function lhs = TSTTHETAFUNJAC1(outstruct, x, t0, t1, A0, B, theta, A1, Jac, res)
 137
 138   DT = t1-t0;
 139   if ( nargin < 10 )
 140     [A1,Jac,res] = asm_build_system(outstruct,x,t1);
 141   endif
 142   lhs = ( (A0+A1)/DT + theta*(B + Jac) );
 143
 144 endfunction
 145 ## Residual for transient problem
 146 function rhs = TSTTHETAFUNRES1(outstruct, x, xold, t0, t1, A0, B, C,
 147                                resold, theta, A1, Jac, res)
 148   DT = t1-t0;
 149   if ( nargin < 13 )
 150     [A1,Jac,res] = asm_build_system(outstruct,x,t1);
 151   endif
 152   rhs = ( (A1+A0)*(x-xold)/DT  + theta * (res + C + B*x) +
 153          (1-theta) * (resold + C + B*xold) );
 154
 155 endfunction
 156 ## Update for transient problem
 157 function update = TSTTHETAFUNUP1(outstruct,x,t1)
 158
 159   [A1,Jac,res] = asm_build_system(outstruct,x,t1);
 160   update = {A1,Jac,res};
 161
 162 endfunction