X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Focs-0.1.3%2Ftst%2Ftst_daspk.m;fp=octave_packages%2Focs-0.1.3%2Ftst%2Ftst_daspk.m;h=92cf69ee1ce0ffaf7cfd885bc9a26c7cdbd3bbac;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/ocs-0.1.3/tst/tst_daspk.m b/octave_packages/ocs-0.1.3/tst/tst_daspk.m new file mode 100644 index 0000000..92cf69e --- /dev/null +++ b/octave_packages/ocs-0.1.3/tst/tst_daspk.m @@ -0,0 +1,135 @@ +## Copyright (C) 2006,2007,2008,2011 Carlo de Falco +## +## This file is part of: +## OCS - A Circuit Simulator for Octave +## +## OCS is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program (see the file LICENSE); if not, +## see . +## +## author: Carlo de Falco + +## -*- texinfo -*- +## +## @deftypefn{Function File} {[@var{out}] =} tst_daspk @ +## (@var{cirstruct},@var{x},@var{t},@var{tol},@var{maxit},@ +## @var{pltvars},@var{verbosity},@var{daskopts},@var{dae_fun}); +## +## Perform a transient simulation of the system described by +## @var{cirstruct} over the time interval @var{t} using @code{daspk}. +## +## The initial value for the state vector is computed by solving a +## steady state problem at @var{t}(1), with starting guess @var{x}. +## +## @var{tol} and @var{maxit} are parameters passed to +## @code{nls_newton_raphson}. +## +## The output @var{out} will contain the value of the state vector +## at each point of @var{t}. +## +## Extra options for @code{daspk} can be passed as name/value pairs in +## the cellarray @var{daskopts}. +## +## The optional parameter @var{verbosity} controls the amount of +## output produced: +## +## @itemize @minus +## @item if verbosity(1) != 0, information on the progress +## of the algorithm are output at runtime +## @item if verbosity(2) != 0, the plot of the variables whose names +## are listed in @var{pltvars} is +## produced after the computation +## @end itemize +## +## For special purposes one may need to pass modified jacobian and +## residual functions. This can be done via the cell array of function +## handles @var{dae_fun}. +## +## Such functions should have the same input and output +## parameter list as the default sub-functions +## TSTBWEFUNJAC0,TSTBWEFUNRES0, TSTBWEFUNJAC,TSTBWEFUNRES. +## +## @seealso{tst_backward_euler,tst_odepkg,tst_theta_method,nls_newton_raphson,daspk} +## +## @end deftypefn + +function [out] = tst_daspk (outstruct, x, t, tol, maxit, pltvars, verbosity, daspkopts, dae_fun) + + ## FIXME: add input check! + if ((nargin < 6) || (nargin > 9)) + error ("tst_daspk: wrong number of input parameters."); + endif + + if (! exist ("verbosity")) + verbosity = [0, 0]; + elseif (length (verbosity) < 2) + verbosity(2) = 0; + endif + + if (verbosity(1)) + fprintf (1, "initial value:\n"); + endif + + daspk_options ("print initial condition info",1); + daspk_options ("maximum order",2); + daspk_options ("initial step size",t(2)-t(1)); + daspk_options ("relative tolerance",1e-3); + + if (nargin > 8) + for ii = 1:2:length (daspkopts) + daspk_options (daspkopts{ii}, daspkopts{ii+1}); + endfor + endif + + + [A0, B, C] = asm_initialize_system (outstruct, x); + + if (nargin > 9) + JAC = @(x) dae_fun{1} (outstruct, x, t(1), B); + RES = @(x) dae_fun{2} (outstruct, x, t(1), B, C); + [x, ii, resnrm] = nls_newton_raphson (x, RES, JAC, tol, maxit, verbosity); + else + [out(:,1), ii] = nls_stationary (outstruct, x, tol, maxit); + endif + + if (nargin > 9) + JAC = @(x) dae_fun{3} (outstruct, x, t(1), B); + RES = @(x) dae_fun{4} (outstruct, x, t(1), B, C); + else + JAC = @(x,xdot,t,c) TSTDASPKFUNJAC (outstruct, x, xdot, A0, B, t, c); + RES = @(x,xdot,t) TSTDASPKFUNRES (outstruct, x, xdot, A0, B, C, t); + endif + + [out, xdot, istate, msg] = daspk ({RES,JAC}, out(:,1), zeros (size (x)), t); + + out = out.'; + if (verbosity(2)) + utl_plot_by_name (t, out, outstruct, pltvars) + endif + +endfunction + +## Jacobian for transient problem +function lhs = TSTDASPKFUNJAC (outstruct, x, xdot, A0, B, t, c) + + [A1, Jac, res] = asm_build_system (outstruct, x, t); + lhs = (c*(A0+A1) + B + Jac); + +endfunction + +## Residual for transient problem +function rhs = TSTDASPKFUNRES (outstruct, x, xdot, A0, B, C, t) + + [A1, Jac, res] = asm_build_system (outstruct, x, t); + rhs = (A0+A1)*xdot + B*x + C + res; + +endfunction \ No newline at end of file