Branch data Line data Source code
1 : : /*
2 : : * transient.cpp - transient helper class implementation
3 : : *
4 : : * Copyright (C) 2004, 2006 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This 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; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : :
33 : : #include "object.h"
34 : : #include "complex.h"
35 : : #include "circuit.h"
36 : : #include "net.h"
37 : : #include "tvector.h"
38 : : #include "tmatrix.h"
39 : : #include "eqnsys.h"
40 : : #include "transient.h"
41 : :
42 : : #define COEFFDEBUG 0
43 : : #define FIXEDCOEFF 0
44 : :
45 : : // Defines where the equivalent admittance coefficient is going to be stored.
46 : : #define COEFF_G 0
47 : :
48 : : namespace qucs {
49 : :
50 : : using namespace transient;
51 : :
52 : : /* The function calculates the integration coefficient for numerical
53 : : integration methods. Supported methods are: Gear (order 1-6),
54 : : Trapezoidal, backward Euler and Adams-Moulton (order 1-6). */
55 : 214550 : void transient::calcCorrectorCoeff (int Method, int order,
56 : : nr_double_t * coefficients,
57 : : nr_double_t * delta) {
58 : :
59 [ + - ]: 214550 : tmatrix<nr_double_t> A (order + 1);
60 [ + - ]: 214550 : tvector<nr_double_t> x (order + 1);
61 [ + - ]: 214550 : tvector<nr_double_t> b (order + 1);
62 : 214550 : eqnsys<nr_double_t> e;
63 : 214550 : e.setAlgo (ALGO_LU_DECOMPOSITION);
64 : :
65 [ + + + - : 214550 : switch (Method) {
- ]
66 : : case INTEGRATOR_GEAR: // GEAR order 1 to 6
67 : : {
68 : : #if FIXEDCOEFF
69 : : int i, r, c;
70 : : // right hand side vector
71 : : for (i = 0; i < order + 1; i++) b.set (i, 1);
72 : : for (i = 1; i < order + 1; i++) {
73 : : A.set (i, 0, i); // first column
74 : : A.set (0, i, 1); // first row
75 : : }
76 : : for (c = 1; c <= order - 1; c++) {
77 : : nr_double_t entry = -c;
78 : : for (r = 1; r <= order; r++) {
79 : : A.set (r, c + 1, entry);
80 : : entry *= -c;
81 : : }
82 : : }
83 : : e.passEquationSys (&A, &x, &b);
84 : : e.solve ();
85 : :
86 : : // vector x consists of b_{-1}, a_{0}, a_{1} ... a_{k-1} right here
87 : : #if COEFFDEBUG
88 : : logprint (LOG_STATUS, "DEBUG: Gear order %d:", order);
89 : : for (i = 0; i < x.getRows (); i++) {
90 : : logprint (LOG_STATUS, " %g", x.get (i));
91 : : }
92 : : logprint (LOG_STATUS, "\n");
93 : : #endif
94 : : nr_double_t k = x.get (0);
95 : : coefficients[COEFF_G] = 1 / delta[0] / k;
96 : : for (i = 1; i <= order; i++) {
97 : : coefficients[i] = - 1 / delta[0] / k * x.get (i);
98 : : }
99 : : #else /* !FIXEDCOEFF */
100 : : int c, r;
101 : : // right hand side vector
102 [ + - ]: 1672 : b.set (1, -1 / delta[0]);
103 : : // first row
104 [ + - ][ + + ]: 5016 : for (c = 0; c < order + 1; c++) A.set (0, c, 1);
105 : : nr_double_t f, a;
106 [ + + ]: 3344 : for (f = 0, c = 0; c < order; c++) {
107 : 1672 : f += delta[c];
108 [ + + ]: 3344 : for (a = 1, r = 0; r < order; r++) {
109 : 1672 : a *= f / delta[0];
110 [ + - ]: 1672 : A.set (r + 1, c + 1, a);
111 : : }
112 : : }
113 [ + - ]: 1672 : e.passEquationSys (&A, &x, &b);
114 [ + - ]: 1672 : e.solve ();
115 [ + - ][ + + ]: 5016 : for (r = 0; r <= order; r++) coefficients[r] = x.get (r);
116 : : #endif /* !FIXEDCOEFF */
117 : : }
118 : 1672 : break;
119 : : case INTEGRATOR_EULER: // BACKWARD EULER
120 : 198 : coefficients[COEFF_G] = 1 / delta[0];
121 : 198 : coefficients[1] = - 1 / delta[0];
122 : 198 : break;
123 : : case INTEGRATOR_TRAPEZOIDAL: // TRAPEZOIDAL (bilinear)
124 : 212680 : coefficients[COEFF_G] = 2 / delta[0];
125 : 212680 : coefficients[1] = - 2 / delta[0];
126 : 212680 : break;
127 : : case INTEGRATOR_ADAMSMOULTON: // ADAMS-MOULTON order 1 to 6
128 : : {
129 : : int i, r, c;
130 : : // right hand side vector
131 [ # # ][ # # ]: 0 : for (i = 0; i < order + 1; i++) b.set (i, 1);
132 [ # # ]: 0 : for (i = 1; i < order + 1; i++) {
133 [ # # ]: 0 : A.set (i, 1, i); // second column
134 [ # # ]: 0 : A.set (1, i, 1); // second row
135 : : }
136 [ # # ]: 0 : A.set (0, 0, 1);
137 [ # # ]: 0 : for (c = 1; c <= order - 2; c++) {
138 : 0 : nr_double_t entry = -c;
139 [ # # ]: 0 : for (r = 2; r <= order; r++) {
140 [ # # ]: 0 : A.set (r, c + 2, r * entry);
141 : 0 : entry *= -c;
142 : : }
143 : : }
144 [ # # ]: 0 : e.passEquationSys (&A, &x, &b);
145 [ # # ]: 0 : e.solve ();
146 : :
147 : : // vector x consists of a_{0}, b_{-1}, b_{0} ... b_{k-2} right here
148 : : #if COEFFDEBUG
149 : : logprint (LOG_STATUS, "DEBUG: Moulton order %d:", order);
150 : : for (i = 0; i < x.getRows (); i++) {
151 : : logprint (LOG_STATUS, " %g", x.get (i));
152 : : }
153 : : logprint (LOG_STATUS, "\n");
154 : : #endif
155 [ # # ]: 0 : nr_double_t k = x.get (1);
156 : 0 : coefficients[COEFF_G] = 1 / delta[0] / k;
157 [ # # ]: 0 : coefficients[1] = -x.get (0) / delta[0] / k;
158 [ # # ]: 0 : for (i = 2; i <= order; i++) {
159 [ # # ]: 0 : coefficients[i] = -x.get (i) / k;
160 : : }
161 : : }
162 : 0 : break;
163 : 214550 : }
164 : 214550 : }
165 : :
166 : : /* The function calculates the integration coefficient for numerical
167 : : integration methods. Supported methods are: Adams-Bashford (order
168 : : 1-6), forward Euler and explicit Gear (order 1-6). */
169 : 214550 : void transient::calcPredictorCoeff (int Method, int order,
170 : : nr_double_t * coefficients,
171 : : nr_double_t * delta) {
172 : :
173 [ + - ]: 214550 : tmatrix<nr_double_t> A (order + 1);
174 [ + - ]: 214550 : tvector<nr_double_t> x (order + 1);
175 [ + - ]: 214550 : tvector<nr_double_t> b (order + 1);
176 : 214550 : eqnsys<nr_double_t> e;
177 : 214550 : e.setAlgo (ALGO_LU_DECOMPOSITION);
178 : :
179 [ + + + - ]: 214550 : switch (Method) {
180 : : case INTEGRATOR_GEAR: // explicit GEAR order 1 to 6
181 : : {
182 : : int c, r;
183 : : // right hand side vector
184 [ + - ]: 1672 : b.set (0, 1);
185 : : // first row
186 [ + - ][ + + ]: 5016 : for (c = 0; c < order + 1; c++) A.set (0, c, 1);
187 : : nr_double_t f, a;
188 [ + + ]: 5016 : for (f = 0, c = 0; c < order + 1; c++) {
189 : 3344 : f += delta[c];
190 [ + + ]: 6688 : for (a = 1, r = 0; r < order; r++) {
191 : 3344 : a *= f / delta[0];
192 [ + - ]: 3344 : A.set (r + 1, c, a);
193 : : }
194 : : }
195 [ + - ]: 1672 : e.passEquationSys (&A, &x, &b);
196 [ + - ]: 1672 : e.solve ();
197 [ + - ][ + + ]: 5016 : for (r = 0; r <= order; r++) coefficients[r] = x.get (r);
198 : : }
199 : 1672 : break;
200 : : case INTEGRATOR_ADAMSBASHFORD: // ADAMS-BASHFORD order 1 to 6
201 : : {
202 : : int i, r, c;
203 : : // right hand side vector
204 [ + - ][ + + ]: 850720 : for (i = 0; i < order + 1; i++) b.set (i, 1);
205 [ + - ][ + + ]: 638040 : for (i = 1; i < order + 1; i++) A.set (1, i, 1); // second row
206 [ + - ]: 212680 : A.set (0, 0, 1);
207 [ + + ]: 425360 : for (c = 1; c <= order - 1; c++) {
208 : 212680 : nr_double_t entry = -c;
209 [ + + ]: 425360 : for (r = 2; r <= order; r++) {
210 [ + - ]: 212680 : A.set (r, c + 1, r * entry);
211 : 212680 : entry *= -c;
212 : : }
213 : : }
214 [ + - ]: 212680 : e.passEquationSys (&A, &x, &b);
215 [ + - ]: 212680 : e.solve ();
216 : :
217 : : // vector x consists of a_{0}, b_{0}, b_{1} ... b_{k-1} right here
218 : : #if COEFFDEBUG
219 : : logprint (LOG_STATUS, "DEBUG: Bashford order %d:", order);
220 : : for (i = 0; i < x.getRows (); i++) {
221 : : logprint (LOG_STATUS, " %g", x.get (i));
222 : : }
223 : : logprint (LOG_STATUS, "\n");
224 : : #endif
225 [ + - ]: 212680 : coefficients[COEFF_G] = x.get (0);
226 [ + + ]: 638040 : for (i = 1; i <= order; i++) {
227 [ + - ]: 425360 : coefficients[i] = x.get (i) * delta[0];
228 : : }
229 : : #if !FIXEDCOEFF
230 [ + - ]: 212680 : if (order == 2) {
231 : 212680 : nr_double_t f = - delta[0] / (2 * delta[1]);
232 : 212680 : coefficients[0] = 1;
233 : 212680 : coefficients[1] = (1 - f) * delta[0];
234 : 212680 : coefficients[2] = f * delta[0];
235 : : }
236 : : #endif
237 : : }
238 : 212680 : break;
239 : : case INTEGRATOR_EULER: // FORWARD EULER
240 : 198 : coefficients[COEFF_G] = 1;
241 : 198 : coefficients[1] = delta[0];
242 : 198 : break;
243 : 214550 : }
244 : 214550 : }
245 : :
246 : : // Loads the equivalent conductance.
247 : 72850 : void transient::getConductance (integrator * c, nr_double_t cap,
248 : : nr_double_t& geq) {
249 : 72850 : nr_double_t * coeff = c->getCoefficients ();
250 : 72850 : geq = cap * coeff[COEFF_G];
251 : 72850 : }
252 : :
253 : : // This is the implicit Euler integrator.
254 : 10073 : void transient::integrateEuler (integrator * c, int qstate, nr_double_t cap,
255 : : nr_double_t& geq, nr_double_t& ceq) {
256 : 10073 : nr_double_t * coeff = c->getCoefficients ();
257 : 10073 : int cstate = qstate + 1;
258 : : nr_double_t cur;
259 : 10073 : geq = cap * coeff[COEFF_G];
260 : 10073 : ceq = c->getState (qstate, 1) * coeff[1];
261 : 10073 : cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
262 : 10073 : c->setState (cstate, cur);
263 : 10073 : }
264 : :
265 : : // Trapezoidal integrator.
266 : 1426826 : void transient::integrateBilinear (integrator * c, int qstate, nr_double_t cap,
267 : : nr_double_t& geq, nr_double_t& ceq) {
268 : 1426826 : nr_double_t * coeff = c->getCoefficients ();
269 : 1426826 : int cstate = qstate + 1;
270 : : nr_double_t cur;
271 : 1426826 : geq = cap * coeff[COEFF_G];
272 : 1426826 : ceq = c->getState (qstate, 1) * coeff[1] - c->getState (cstate, 1);
273 : 1426826 : cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
274 : 1426826 : c->setState (cstate, cur);
275 : 1426826 : }
276 : :
277 : : // Integrator using the Gear coefficients.
278 : 25417 : void transient::integrateGear (integrator * c, int qstate, nr_double_t cap,
279 : : nr_double_t& geq, nr_double_t& ceq) {
280 : 25417 : nr_double_t * coeff = c->getCoefficients ();
281 : 25417 : int i, cstate = qstate + 1;
282 : : nr_double_t cur;
283 : 25417 : geq = cap * coeff[COEFF_G];
284 [ + + ]: 50834 : for (ceq = 0, i = 1; i <= c->getOrder (); i++) {
285 : 25417 : ceq += c->getState (qstate, i) * coeff[i];
286 : : }
287 : 25417 : cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
288 : 25417 : c->setState (cstate, cur);
289 : 25417 : }
290 : :
291 : : // Integrator using the Adams-Moulton coefficients.
292 : 0 : void transient::integrateMoulton (integrator * c, int qstate, nr_double_t cap,
293 : : nr_double_t& geq, nr_double_t& ceq) {
294 : 0 : nr_double_t * coeff = c->getCoefficients ();
295 : 0 : int i, cstate = qstate + 1;
296 : : nr_double_t cur;
297 : 0 : geq = cap * coeff[COEFF_G];
298 : 0 : ceq = c->getState (qstate, 1) * coeff[1];
299 [ # # ]: 0 : for (i = 2; i <= c->getOrder (); i++) {
300 : 0 : ceq += c->getState (cstate, i - 1) * coeff[i];
301 : : }
302 : 0 : cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
303 : 0 : c->setState (cstate, cur);
304 : 0 : }
305 : :
306 : : /* The function applies the appropriate integration function to the
307 : : given circuit object. */
308 : 2467 : void transient::setIntegrationMethod (circuit * c, int Method) {
309 [ + + + - : 2467 : switch (Method) {
- ]
310 : : case INTEGRATOR_GEAR:
311 : 225 : c->setIntegration (integrateGear);
312 : 225 : break;
313 : : case INTEGRATOR_TRAPEZOIDAL:
314 : 1516 : c->setIntegration (integrateBilinear);
315 : 1516 : break;
316 : : case INTEGRATOR_EULER:
317 : 726 : c->setIntegration (integrateEuler);
318 : 726 : break;
319 : : case INTEGRATOR_ADAMSMOULTON:
320 : 0 : c->setIntegration (integrateMoulton);
321 : 0 : break;
322 : : default:
323 : 0 : c->setIntegration (NULL);
324 : 0 : break;
325 : : }
326 : 2467 : c->setConductance (getConductance);
327 : 2467 : }
328 : :
329 : : /* Returns an appropriate integrator type identifier and the maximum
330 : : order depending on the given string argument. */
331 : 65 : int transient::correctorType (const char * const Method, int& MaxOrder) {
332 [ + + ]: 65 : if (!strcmp (Method, "Gear")) {
333 [ - + ]: 1 : if (MaxOrder > 6) MaxOrder = 6;
334 [ - + ]: 1 : if (MaxOrder < 1) MaxOrder = 1;
335 : 1 : return INTEGRATOR_GEAR;
336 : : }
337 [ + - ]: 64 : else if (!strcmp (Method, "Trapezoidal")) {
338 : 64 : MaxOrder = 2;
339 : 64 : return INTEGRATOR_TRAPEZOIDAL;
340 : : }
341 [ # # ]: 0 : else if (!strcmp (Method, "Euler")) {
342 : 0 : MaxOrder = 1;
343 : 0 : return INTEGRATOR_EULER;
344 : : }
345 [ # # ]: 0 : else if (!strcmp (Method, "AdamsMoulton")) {
346 [ # # ]: 0 : if (MaxOrder > 6) MaxOrder = 6;
347 [ # # ]: 0 : if (MaxOrder < 1) MaxOrder = 1;
348 : 0 : return INTEGRATOR_ADAMSMOULTON;
349 : : }
350 [ # # ]: 0 : else if (!strcmp (Method, "AdamsBashford")) {
351 [ # # ]: 0 : if (MaxOrder > 6) MaxOrder = 6;
352 [ # # ]: 0 : if (MaxOrder < 1) MaxOrder = 1;
353 : 0 : return INTEGRATOR_ADAMSBASHFORD;
354 : : }
355 : 65 : return INTEGRATOR_UNKNOWN;
356 : : }
357 : :
358 : : /* The function returns the appropriate predictor integration method
359 : : for the given corrector method and adjusts the order of the
360 : : predictor as well based on the given corrector method. */
361 : 208 : int transient::predictorType (int corrMethod, int corrOrder, int& predOrder) {
362 : 208 : int predMethod = INTEGRATOR_UNKNOWN;
363 [ + - + + : 208 : switch (corrMethod) {
- ]
364 : : case INTEGRATOR_GEAR:
365 : 16 : predMethod = INTEGRATOR_GEAR;
366 : 16 : break;
367 : : case INTEGRATOR_ADAMSMOULTON:
368 : 0 : predMethod = INTEGRATOR_ADAMSBASHFORD;
369 : 0 : break;
370 : : case INTEGRATOR_TRAPEZOIDAL:
371 : 128 : predMethod = INTEGRATOR_ADAMSBASHFORD;
372 : 128 : break;
373 : : case INTEGRATOR_EULER:
374 : 64 : predMethod = INTEGRATOR_EULER;
375 : 64 : break;
376 : : }
377 : 208 : predOrder = corrOrder;
378 : 208 : return predMethod;
379 : : }
380 : :
381 : : // Structure defining integration algorithm for each possible order.
382 : : struct integration_types_t {
383 : : int Method;
384 : : int integratorType[6];
385 : : nr_double_t corrErrorConstant[6];
386 : : nr_double_t predErrorConstant[6];
387 : : };
388 : :
389 : : static struct integration_types_t integration_types[] = {
390 : : { INTEGRATOR_EULER,
391 : : { INTEGRATOR_EULER },
392 : : { -1.0/2 },
393 : : { +1.0/2 }
394 : : },
395 : : { INTEGRATOR_TRAPEZOIDAL,
396 : : { INTEGRATOR_EULER, INTEGRATOR_TRAPEZOIDAL },
397 : : { -1.0/2, -1.0/12 },
398 : : { +1.0/2, +5.0/12 }
399 : : },
400 : : { INTEGRATOR_GEAR,
401 : : { INTEGRATOR_GEAR, INTEGRATOR_GEAR, INTEGRATOR_GEAR,
402 : : INTEGRATOR_GEAR, INTEGRATOR_GEAR, INTEGRATOR_GEAR },
403 : : { -1.0/2, -2.0/9, -3.0/22, -12.0/125, -10.0/137, -20.0/343 },
404 : : { +1.0, +1.0, +1.0, +1.0, +1.0, +1.0 }
405 : : },
406 : : { INTEGRATOR_ADAMSMOULTON,
407 : : { INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON,
408 : : INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON,
409 : : INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON },
410 : : { -1.0/2, -1.0/12, -1.0/24, -19.0/720, -3.0/160, -863.0/60480 },
411 : : { +1.0/2, +1.0/12, +1.0/24, +19.0/720, +3.0/160, +863.0/60480 }
412 : : },
413 : : { INTEGRATOR_ADAMSBASHFORD,
414 : : { INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD,
415 : : INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD,
416 : : INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD },
417 : : { -1.0/2, -5.0/12, -3.0/8, -251.0/720, -95.0/288, -19087.0/60480 },
418 : : { +1.0/2, +5.0/12, +3.0/8, +251.0/720, +95.0/288, +19087.0/60480 }
419 : : }
420 : : };
421 : :
422 : : /* The function returns the appropriate integration type for the given
423 : : corrector integration type and order. */
424 : 143 : int transient::correctorType (int Method, int order) {
425 : 143 : return integration_types[Method].integratorType[order - 1];
426 : : }
427 : :
428 : : // Returns the error constant for the given corrector.
429 : 214290 : nr_double_t transient::getCorrectorError (int Method, int order) {
430 : 214290 : return integration_types[Method].corrErrorConstant[order - 1];
431 : : }
432 : :
433 : : // Returns the error constant for the given predictor.
434 : 214290 : nr_double_t transient::getPredictorError (int Method, int order) {
435 : 214290 : return integration_types[Method].predErrorConstant[order - 1];
436 : : }
437 : :
438 : : }
|
