Qucs-core  0.0.19
transient.cpp
Go to the documentation of this file.
00001 /*
00002  * transient.cpp - transient helper class implementation
00003  *
00004  * Copyright (C) 2004, 2006 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 
00033 #include "object.h"
00034 #include "complex.h"
00035 #include "circuit.h"
00036 #include "net.h"
00037 #include "tvector.h"
00038 #include "tmatrix.h"
00039 #include "eqnsys.h"
00040 #include "transient.h"
00041 
00042 #define COEFFDEBUG 0
00043 #define FIXEDCOEFF 0
00044 
00045 // Defines where the equivalent admittance coefficient is going to be stored.
00046 #define COEFF_G 0
00047 
00048 namespace qucs {
00049 
00050 using namespace transient;
00051 
00052 /* The function calculates the integration coefficient for numerical
00053    integration methods.  Supported methods are: Gear (order 1-6),
00054    Trapezoidal, backward Euler and Adams-Moulton (order 1-6). */
00055 void transient::calcCorrectorCoeff (int Method, int order,
00056                                     nr_double_t * coefficients,
00057                                     nr_double_t * delta) {
00058 
00059   tmatrix<nr_double_t> A (order + 1);
00060   tvector<nr_double_t> x (order + 1);
00061   tvector<nr_double_t> b (order + 1);
00062   eqnsys<nr_double_t> e;
00063   e.setAlgo (ALGO_LU_DECOMPOSITION);
00064 
00065   switch (Method) {
00066   case INTEGRATOR_GEAR: // GEAR order 1 to 6
00067     {
00068 #if FIXEDCOEFF
00069       int i, r, c;
00070       // right hand side vector
00071       for (i = 0; i < order + 1; i++) b.set (i, 1);
00072       for (i = 1; i < order + 1; i++) {
00073         A.set (i, 0, i); // first column
00074         A.set (0, i, 1); // first row
00075       }
00076       for (c = 1; c <= order - 1; c++) {
00077         nr_double_t entry = -c;
00078         for (r = 1; r <= order; r++) {
00079           A.set (r, c + 1, entry);
00080           entry *= -c;
00081         }
00082       }
00083       e.passEquationSys (&A, &x, &b);
00084       e.solve ();
00085 
00086       // vector x consists of b_{-1}, a_{0}, a_{1} ... a_{k-1} right here
00087 #if COEFFDEBUG
00088       logprint (LOG_STATUS, "DEBUG: Gear order %d:", order);
00089       for (i = 0; i < x.getRows (); i++) {
00090         logprint (LOG_STATUS, " %g", x.get (i));
00091       }
00092       logprint (LOG_STATUS, "\n");
00093 #endif
00094       nr_double_t k = x.get (0);
00095       coefficients[COEFF_G] = 1 / delta[0] / k;
00096       for (i = 1; i <= order; i++) {
00097         coefficients[i] = - 1 / delta[0] / k * x.get (i);
00098       }
00099 #else /* !FIXEDCOEFF */
00100       int c, r;
00101       // right hand side vector
00102       b.set (1, -1 / delta[0]);
00103       // first row
00104       for (c = 0; c < order + 1; c++) A.set (0, c, 1);
00105       nr_double_t f, a;
00106       for (f = 0, c = 0; c < order; c++) {
00107         f += delta[c];
00108         for (a = 1, r = 0; r < order; r++) {
00109           a *= f / delta[0];
00110           A.set (r + 1, c + 1, a);
00111         }
00112       }
00113       e.passEquationSys (&A, &x, &b);
00114       e.solve ();
00115       for (r = 0; r <= order; r++) coefficients[r] = x.get (r);
00116 #endif /* !FIXEDCOEFF */
00117     }
00118     break;
00119   case INTEGRATOR_EULER: // BACKWARD EULER
00120     coefficients[COEFF_G] = 1 / delta[0];
00121     coefficients[1] = - 1 / delta[0];
00122     break;
00123   case INTEGRATOR_TRAPEZOIDAL: // TRAPEZOIDAL (bilinear)
00124     coefficients[COEFF_G] = 2 / delta[0];
00125     coefficients[1] = - 2 / delta[0];
00126     break;
00127   case INTEGRATOR_ADAMSMOULTON: // ADAMS-MOULTON order 1 to 6
00128     {
00129       int i, r, c;
00130       // right hand side vector
00131       for (i = 0; i < order + 1; i++) b.set (i, 1);
00132       for (i = 1; i < order + 1; i++) {
00133         A.set (i, 1, i); // second column
00134         A.set (1, i, 1); // second row
00135       }
00136       A.set (0, 0, 1);
00137       for (c = 1; c <= order - 2; c++) {
00138         nr_double_t entry = -c;
00139         for (r = 2; r <= order; r++) {
00140           A.set (r, c + 2, r * entry);
00141           entry *= -c;
00142         }
00143       }
00144       e.passEquationSys (&A, &x, &b);
00145       e.solve ();
00146 
00147       // vector x consists of a_{0}, b_{-1}, b_{0} ... b_{k-2} right here
00148 #if COEFFDEBUG
00149       logprint (LOG_STATUS, "DEBUG: Moulton order %d:", order);
00150       for (i = 0; i < x.getRows (); i++) {
00151         logprint (LOG_STATUS, " %g", x.get (i));
00152       }
00153       logprint (LOG_STATUS, "\n");
00154 #endif
00155       nr_double_t k = x.get (1);
00156       coefficients[COEFF_G] = 1 / delta[0] / k;
00157       coefficients[1] = -x.get (0) / delta[0] / k;
00158       for (i = 2; i <= order; i++) {
00159         coefficients[i] = -x.get (i) / k;
00160       }
00161     }
00162     break;
00163   }
00164 }
00165 
00166 /* The function calculates the integration coefficient for numerical
00167    integration methods.  Supported methods are: Adams-Bashford (order
00168    1-6), forward Euler and explicit Gear (order 1-6). */
00169 void transient::calcPredictorCoeff (int Method, int order,
00170                                     nr_double_t * coefficients,
00171                                     nr_double_t * delta) {
00172 
00173   tmatrix<nr_double_t> A (order + 1);
00174   tvector<nr_double_t> x (order + 1);
00175   tvector<nr_double_t> b (order + 1);
00176   eqnsys<nr_double_t> e;
00177   e.setAlgo (ALGO_LU_DECOMPOSITION);
00178 
00179   switch (Method) {
00180   case INTEGRATOR_GEAR: // explicit GEAR order 1 to 6
00181     {
00182       int c, r;
00183       // right hand side vector
00184       b.set (0, 1);
00185       // first row
00186       for (c = 0; c < order + 1; c++) A.set (0, c, 1);
00187       nr_double_t f, a;
00188       for (f = 0, c = 0; c < order + 1; c++) {
00189         f += delta[c];
00190         for (a = 1, r = 0; r < order; r++) {
00191           a *= f / delta[0];
00192           A.set (r + 1, c, a);
00193         }
00194       }
00195       e.passEquationSys (&A, &x, &b);
00196       e.solve ();
00197       for (r = 0; r <= order; r++) coefficients[r] = x.get (r);
00198     }
00199     break;
00200   case INTEGRATOR_ADAMSBASHFORD: // ADAMS-BASHFORD order 1 to 6
00201     {
00202       int i, r, c;
00203       // right hand side vector
00204       for (i = 0; i < order + 1; i++) b.set (i, 1);
00205       for (i = 1; i < order + 1; i++) A.set (1, i, 1); // second row
00206       A.set (0, 0, 1);
00207       for (c = 1; c <= order - 1; c++) {
00208         nr_double_t entry = -c;
00209         for (r = 2; r <= order; r++) {
00210           A.set (r, c + 1, r * entry);
00211           entry *= -c;
00212         }
00213       }
00214       e.passEquationSys (&A, &x, &b);
00215       e.solve ();
00216 
00217       // vector x consists of a_{0}, b_{0}, b_{1} ... b_{k-1} right here
00218 #if COEFFDEBUG
00219       logprint (LOG_STATUS, "DEBUG: Bashford order %d:", order);
00220       for (i = 0; i < x.getRows (); i++) {
00221         logprint (LOG_STATUS, " %g", x.get (i));
00222       }
00223       logprint (LOG_STATUS, "\n");
00224 #endif
00225       coefficients[COEFF_G] = x.get (0);
00226       for (i = 1; i <= order; i++) {
00227         coefficients[i] = x.get (i) * delta[0];
00228       }
00229 #if !FIXEDCOEFF
00230       if (order == 2) {
00231         nr_double_t f = - delta[0] / (2 * delta[1]);
00232         coefficients[0] = 1;
00233         coefficients[1] = (1 - f) * delta[0];
00234         coefficients[2] = f * delta[0];
00235       }
00236 #endif
00237     }
00238     break;
00239   case INTEGRATOR_EULER: // FORWARD EULER
00240     coefficients[COEFF_G] = 1;
00241     coefficients[1] = delta[0];
00242     break;
00243   }
00244 }
00245 
00246 // Loads the equivalent conductance.
00247 void transient::getConductance (integrator * c, nr_double_t cap,
00248                                 nr_double_t& geq) {
00249   nr_double_t * coeff = c->getCoefficients ();
00250   geq = cap * coeff[COEFF_G];
00251 }
00252 
00253 // This is the implicit Euler integrator.
00254 void transient::integrateEuler (integrator * c, int qstate, nr_double_t cap,
00255                                 nr_double_t& geq, nr_double_t& ceq) {
00256   nr_double_t * coeff = c->getCoefficients ();
00257   int cstate = qstate + 1;
00258   nr_double_t cur;
00259   geq = cap * coeff[COEFF_G];
00260   ceq = c->getState (qstate, 1) * coeff[1];
00261   cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
00262   c->setState (cstate, cur);
00263 }
00264 
00265 // Trapezoidal integrator.
00266 void transient::integrateBilinear (integrator * c, int qstate, nr_double_t cap,
00267                                    nr_double_t& geq, nr_double_t& ceq) {
00268   nr_double_t * coeff = c->getCoefficients ();
00269   int cstate = qstate + 1;
00270   nr_double_t cur;
00271   geq = cap * coeff[COEFF_G];
00272   ceq = c->getState (qstate, 1) * coeff[1] - c->getState (cstate, 1);
00273   cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
00274   c->setState (cstate, cur);
00275 }
00276 
00277 // Integrator using the Gear coefficients.
00278 void transient::integrateGear (integrator * c, int qstate, nr_double_t cap,
00279                                nr_double_t& geq, nr_double_t& ceq) {
00280   nr_double_t * coeff = c->getCoefficients ();
00281   int i, cstate = qstate + 1;
00282   nr_double_t cur;
00283   geq = cap * coeff[COEFF_G];
00284   for (ceq = 0, i = 1; i <= c->getOrder (); i++) {
00285     ceq += c->getState (qstate, i) * coeff[i];
00286   }
00287   cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
00288   c->setState (cstate, cur);
00289 }
00290 
00291 // Integrator using the Adams-Moulton coefficients.
00292 void transient::integrateMoulton (integrator * c, int qstate, nr_double_t cap,
00293                                   nr_double_t& geq, nr_double_t& ceq) {
00294   nr_double_t * coeff = c->getCoefficients ();
00295   int i, cstate = qstate + 1;
00296   nr_double_t cur;
00297   geq = cap * coeff[COEFF_G];
00298   ceq = c->getState (qstate, 1) * coeff[1];
00299   for (i = 2; i <= c->getOrder (); i++) {
00300     ceq += c->getState (cstate, i - 1) * coeff[i];
00301   }
00302   cur = c->getState (qstate) * coeff[COEFF_G] + ceq;
00303   c->setState (cstate, cur);
00304 }
00305 
00306 /* The function applies the appropriate integration function to the
00307    given circuit object. */
00308 void transient::setIntegrationMethod (circuit * c, int Method) {
00309   switch (Method) {
00310   case INTEGRATOR_GEAR:
00311     c->setIntegration (integrateGear);
00312     break;
00313   case INTEGRATOR_TRAPEZOIDAL:
00314     c->setIntegration (integrateBilinear);
00315     break;
00316   case INTEGRATOR_EULER:
00317     c->setIntegration (integrateEuler);
00318     break;
00319   case INTEGRATOR_ADAMSMOULTON:
00320     c->setIntegration (integrateMoulton);
00321     break;
00322   default:
00323     c->setIntegration (NULL);
00324     break;
00325   }
00326   c->setConductance (getConductance);
00327 }
00328 
00329 /* Returns an appropriate integrator type identifier and the maximum
00330    order depending on the given string argument. */
00331 int transient::correctorType (const char * const Method, int& MaxOrder) {
00332   if (!strcmp (Method, "Gear")) {
00333     if (MaxOrder > 6) MaxOrder = 6;
00334     if (MaxOrder < 1) MaxOrder = 1;
00335     return INTEGRATOR_GEAR;
00336   }
00337   else if (!strcmp (Method, "Trapezoidal")) {
00338     MaxOrder = 2;
00339     return INTEGRATOR_TRAPEZOIDAL;
00340   }
00341   else if (!strcmp (Method, "Euler")) {
00342     MaxOrder = 1;
00343     return INTEGRATOR_EULER;
00344   }
00345   else if (!strcmp (Method, "AdamsMoulton")) {
00346     if (MaxOrder > 6) MaxOrder = 6;
00347     if (MaxOrder < 1) MaxOrder = 1;
00348     return INTEGRATOR_ADAMSMOULTON;
00349   }
00350   else if (!strcmp (Method, "AdamsBashford")) {
00351     if (MaxOrder > 6) MaxOrder = 6;
00352     if (MaxOrder < 1) MaxOrder = 1;
00353     return INTEGRATOR_ADAMSBASHFORD;
00354   }
00355   return INTEGRATOR_UNKNOWN;
00356 }
00357 
00358 /* The function returns the appropriate predictor integration method
00359    for the given corrector method and adjusts the order of the
00360    predictor as well based on the given corrector method. */
00361 int transient::predictorType (int corrMethod, int corrOrder, int& predOrder) {
00362   int predMethod = INTEGRATOR_UNKNOWN;
00363   switch (corrMethod) {
00364   case INTEGRATOR_GEAR:
00365     predMethod = INTEGRATOR_GEAR;
00366     break;
00367   case INTEGRATOR_ADAMSMOULTON:
00368     predMethod = INTEGRATOR_ADAMSBASHFORD;
00369     break;
00370   case INTEGRATOR_TRAPEZOIDAL:
00371     predMethod = INTEGRATOR_ADAMSBASHFORD;
00372     break;
00373   case INTEGRATOR_EULER:
00374     predMethod = INTEGRATOR_EULER;
00375     break;
00376   }
00377   predOrder = corrOrder;
00378   return predMethod;
00379 }
00380 
00381 // Structure defining integration algorithm for each possible order.
00382 struct integration_types_t {
00383   int Method;
00384   int integratorType[6];
00385   nr_double_t corrErrorConstant[6];
00386   nr_double_t predErrorConstant[6];
00387 };
00388 
00389 static struct integration_types_t integration_types[] = {
00390   { INTEGRATOR_EULER,
00391     { INTEGRATOR_EULER },
00392     { -1.0/2 },
00393     { +1.0/2 }
00394   },
00395   { INTEGRATOR_TRAPEZOIDAL,
00396     { INTEGRATOR_EULER, INTEGRATOR_TRAPEZOIDAL },
00397     { -1.0/2, -1.0/12 },
00398     { +1.0/2, +5.0/12 }
00399   },
00400   { INTEGRATOR_GEAR,
00401     { INTEGRATOR_GEAR, INTEGRATOR_GEAR, INTEGRATOR_GEAR,
00402       INTEGRATOR_GEAR, INTEGRATOR_GEAR, INTEGRATOR_GEAR },
00403     { -1.0/2, -2.0/9, -3.0/22, -12.0/125, -10.0/137, -20.0/343 },
00404     { +1.0, +1.0, +1.0, +1.0, +1.0, +1.0 }
00405   },
00406   { INTEGRATOR_ADAMSMOULTON,
00407     { INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON,
00408       INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON,
00409       INTEGRATOR_ADAMSMOULTON, INTEGRATOR_ADAMSMOULTON },
00410     { -1.0/2, -1.0/12, -1.0/24, -19.0/720, -3.0/160, -863.0/60480 },
00411     { +1.0/2, +1.0/12, +1.0/24, +19.0/720, +3.0/160, +863.0/60480 }
00412   },
00413   { INTEGRATOR_ADAMSBASHFORD,
00414     { INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD,
00415       INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD,
00416       INTEGRATOR_ADAMSBASHFORD, INTEGRATOR_ADAMSBASHFORD },
00417     { -1.0/2, -5.0/12, -3.0/8, -251.0/720, -95.0/288, -19087.0/60480 },
00418     { +1.0/2, +5.0/12, +3.0/8, +251.0/720, +95.0/288, +19087.0/60480 }
00419   }
00420 };
00421 
00422 /* The function returns the appropriate integration type for the given
00423    corrector integration type and order. */
00424 int transient::correctorType (int Method, int order) {
00425   return integration_types[Method].integratorType[order - 1];
00426 }
00427 
00428 // Returns the error constant for the given corrector.
00429 nr_double_t transient::getCorrectorError (int Method, int order) {
00430   return integration_types[Method].corrErrorConstant[order - 1];
00431 }
00432 
00433 // Returns the error constant for the given predictor.
00434 nr_double_t transient::getPredictorError (int Method, int order) {
00435   return integration_types[Method].predErrorConstant[order - 1];
00436 }
00437 
00438 }