Qucs-core
0.0.19
|
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 }