1 ## Copyright (C) 2006,2007,2008 Carlo de Falco
3 ## This file is part of:
4 ## OCS - A Circuit Simulator for Octave
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.
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.
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/>.
19 ## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
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});
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}.
31 ## Pssible values for @var{solver} are @code{ode2r}, @code{ode5r},
32 ## @code{oders} or @code{odesx}.
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}.
37 ## @var{tol} and @var{maxit} are parameters passed to @code{nls_newton_raphson}.
39 ## If one output is requested @var{out} will contain the value of the state vector
40 ## at each point of @var{t}.
42 ## If two outputs are requested @var{out} will contain the value of the state vector
43 ## at each point of @var{tout}.
45 ## Extra options for options for the solver can be passed to the solver
46 ## via @var{odestruct}.
48 ## The optional parameter @var{verbosity} controls the amount of
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
59 ## @seealso{tst_backward_euler,tst_theta_method,tst_daspk,nls_newton_raphson,odepkg,odeset,@
60 ## ode2r,ode5r,oders,odesx}
64 function [out, tout] = tst_odepkg (outstruct,x,t,tol,maxit,\
65 pltvars,solver,verbosity,odestruct)
68 ## FIXME: add input check!
69 if ((nargin < 7) || (nargin > 9))
70 error("tst_odepkg:wrong number of input parameters");
73 if ~exist("verbosity")
75 elseif length(verbosity)<2
80 fprintf(1,"initial value:\n");
83 if ~exist("odestruct")
87 [A0,B,C] = asm_initialize_system(outstruct,x);
89 [out(:,1),ii] = nls_stationary(outstruct,x,tol,maxit);
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);
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));
101 odestruct = odeset(odestruct, "OutputFcn",
102 @(t, y, flag) plotfun(t, y, flag, outstruct, pltvars) );
105 [tout, out] = feval( solver, RES, t([1 end]), x, odestruct);
107 out = interp1(tout, out, t).';
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
116 %# No input argument check is done for a higher processing speed
117 persistent vfigure; persistent vtold;
118 persistent vyold; persistent vcounter;
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);
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;
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");
143 ## Jacobian for transient problems
144 function lhs = TSTODEPKGFUNJAC(outstruct,x,A0,B,t)
146 [A1,Jac,res] = asm_build_system(outstruct,x,t);
151 function lhs = TSTODEPKGFUNMASS(outstruct,x,A0,t)
153 [A1,Jac,res] = asm_build_system(outstruct,x,t);
158 ## Residual for transient problems
159 function rhs = TSTODEPKGFUNRES(outstruct,x,A0,B,C,t)
161 [A1,Jac,res] = asm_build_system(outstruct,x,t);