1 ## Copyright (C) 2006,2007,2008,2011 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{niter}] =} @
24 ## tst_backward_euler(@var{cirstruct},@var{x},@var{t},@var{tol},@
25 ## @var{maxit},@var{pltvars},@var{verbosity},@var{dae_fun})
27 ## Perform a transient simulation of the system described by
28 ## @var{cirstruct} over the time interval @var{t} using the backward
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}.
34 ## @var{tol} and @var{maxit} are parameters passed to
35 ## @code{nls_newton_raphson}.
37 ## The output @var{out} will contain the value of the state vector at
38 ## each point of @var{t}.
40 ## The optional parameter @var{verbosity} controls the amount of
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
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}.
55 ## Such functions should have the same input and output
56 ## parameter list as the default sub-functions
57 ## TSTBWEFUNJAC0,TSTBWEFUNRES0, TSTBWEFUNJAC,TSTBWEFUNRES.
59 ## The optional output @var{niter} returns the number of Newton iterations
60 ## needed to reach convergence.
62 ## @seealso{tst_daspk,tst_theta_method,tst_odepkg,nls_newton_raphson}
66 function [out, varargout] = tst_backward_euler (outstruct, x, t, tol, maxit, pltvars, verbosity, dae_fun)
69 ## FIXME: add input check!
70 if ((nargin < 6) || (nargin > 8))
71 error ("tst_backward_euler: wrong number of input parameters.");
74 if ~exist ("verbosity")
76 elseif (length (verbosity) < 2)
80 out = zeros (rows (x), columns (t));
84 niter = zeros (length(t),1);
88 fprintf (1, "Initial value.\n");
92 [A0, B, C] = asm_initialize_system (outstruct, x);
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));
99 [out(:,1),ii] = nls_stationary (outstruct, x, tol, maxit);
106 for it = 2:length (t)
109 fprintf (1,"Timestep #%d.\n",it);
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);
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,
121 UPDT = @(x) TSTBWEFUNUP1 (outstruct, x, t(it));
124 [out(:,it),ii,resnrm] = nls_newton_raphson (out(:,it-1), RES, JAC, tol, maxit, verbosity(1), UPDT);
131 utl_plot_by_name (t(1:it), out(:,1:it), outstruct, pltvars);
136 ## FIXME: maintain this part?
137 if exist("~/.stop_ocs","file")
138 printf("stopping at timestep %d\n",it);
139 unix("rm ~/.stop_ocs");
146 varargout{1} = niter;
151 ## Jacobian for transient problem
152 function lhs = TSTBWEFUNJAC1 (outstruct, x, t0, t1, A0, B, A1, Jac, res)
156 [A1, Jac, res] = asm_build_system (outstruct, x, t1);
158 lhs = ((A0+A1)/DT + B + Jac);
162 ## Residual for transient problem
163 function rhs = TSTBWEFUNRES1 (outstruct, x, xold, t0, t1, A0, B, C, A1, Jac, res)
167 [A1, Jac, res] = asm_build_system (outstruct, x, t1);
169 rhs = (res + C + B*x + (A0+A1)*(x-xold)/DT);
173 ## Update for transient problem
174 function update = TSTBWEFUNUP1 (outstruct, x, t1)
176 [A1, Jac, res] = asm_build_system (outstruct, x, t1);
177 update = {A1, Jac, res};