Merge pull request #203 from in3otd/issue_34

Fixed errors due to numerical precision.

qucs-core/src/components/tswitch.cpp

@@ -22,6 +22,15 @@
  *
  */
 
+/*!
+ * \file tswitch.cpp
+ * \brief time controlled switch class implementation
+ *
+ * \todo add a property to allow choosing the resistance profile
+ * (cubic spline, linear, abrupt)
+ */
+
+
 #if HAVE_CONFIG_H
 # include <config.h>
 #endif
@@ -86,42 +95,43 @@ void tswitch::initTR (void) {
   qucs::vector * values = getPropertyVector ("time");
   // Find the total time of the switching periods
   T = real (sum (*values));
-  // if the user enters an even number of switchng times
+  // if the user enters an even number of switching times
   // the pattern is repeated continuously
   repeat = (values->getSize () % 2) == 0 ? true : false;
   // make the time taken to go from fully on to fully off
   // the smallest switching time / 100, or the smallest possible
-  // number, but no bogger than the max specified duration
+  // number, but no bigger than the max specified duration
   nr_double_t maxduration = getPropertyDouble("MaxDuration");
   duration = std::min ( std::max (10*NR_TINY, values->minimum() / 100),
                         maxduration );
-
   initDC ();
 }
 
 void tswitch::calcTR (nr_double_t t) {
   const char * const init = getPropertyString ("init");
   nr_double_t ron = getPropertyDouble ("Ron");
   nr_double_t roff = getPropertyDouble ("Roff");
+  const char * const trans_type = getPropertyString ("Transition");
   nr_double_t r = 0;
   nr_double_t rdiff = 0;
-  nr_double_t s_i = 0;
+  //nr_double_t s_i = 0;
   nr_double_t r_0 = 0;
   qucs::vector * values = getPropertyVector ("time");
   bool on = !strcmp (init, "on");
   nr_double_t ti = 0;
 
   if (repeat) {
-    // if the user enters an even number of switchng times
-    // the pattern is repeated continuously. This is acieved by
+    // if the user enters an even number of switching times
+    // the pattern is repeated continuously. This is achieved by
     // subtracting an integer number of total switching periods
     // from the real time
     t = t - T * qucs::floor (t / T);
   }
 
-  // Initialise the last switching time to be a full
-  // switching duration
-  nr_double_t ts = t - duration;
+  // Initialise the last switching time to be well in the past
+  // to avoid having the switch even partially in a transition
+  // state (due to inaccurately computed time differences)
+  nr_double_t ts = -2*duration;
 
   // here we determine whether a switching event should occur
   // by looping through the list of switching times and comparing
@@ -145,39 +155,46 @@ void tswitch::calcTR (nr_double_t t) {
       break;
     }
   }
-  // calculate the time since the last switch occured
-  nr_double_t tdiff = std::max(NR_TINY, t - ts);
-
-  // set the time difference to be no more than the max switch
-  // duration so when we interpolate below we only get the max
-  // or min function value if we are past a switching time
-  if (tdiff > duration) {
-        tdiff = duration;
-  }
-
-  // Set the appropriate resistance. The resistance is interpolated
-  // along a cubic spline with zero derivative at the start and end
-  // points to ensure a smooth derivative
-  if (on)
-  {
-    r_0 = roff;
 
-    rdiff = ron - roff;
-
-    s_i = (rdiff) / (duration);
-  }
-  else
-  {
-    r_0 = ron;
-
-    rdiff = roff - ron;
-
-    s_i = (rdiff) / (duration);
+  if (!strcmp(trans_type, "abrupt")) {
+    r = (on ? ron : roff);
+  } else {
+    // calculate the time since the last switch occurred
+    nr_double_t tdiff = std::max(NR_TINY, t - ts);
+
+    // set the time difference to be no more than the max switch
+    // duration so when we interpolate below we only get the max
+    // or min function value if we are past a switching time
+    if (tdiff > duration) {
+      tdiff = duration;
+    }
+    // Set the appropriate resistance. 
+    if (on) {
+      r_0 = roff;
+      rdiff = ron - roff;
+      //    s_i = (rdiff) / (duration);
+    } else {
+      r_0 = ron;
+      rdiff = roff - ron;
+      //  s_i = (rdiff) / (duration);
+    }
+    if (!strcmp(trans_type, "linear")) {
+      // simple linear transition over the transition time
+      r = r_0 + rdiff * tdiff / duration;
+    } else { // assume trans_type is "spline"
+	// the resistance is interpolated along a constrained cubic spline
+	// with zero derivative at the start and end points to ensure a 
+	// smooth derivative
+	//r = r_0 + ((3. * s_i * qucs::pow (tdiff,2.0)) / (duration))
+	//        + ((-2. * s_i * qucs::pow (tdiff,3.0)) / qucs::pow (duration, 2.0));
+	// use Horner's rule to reduce the numerical errors
+	r = r_0 + (((-2. * rdiff * tdiff / duration) + 3. * rdiff) * qucs::pow(tdiff/duration, 2.0));
+      }
   }
 
-  // perform the interpolation of the constrained cubic spline
-  r = r_0 + ((3. * s_i * qucs::pow (tdiff,2.0)) / (duration))
-          + ((-2. * s_i * qucs::pow (tdiff,3.0)) / qucs::pow (duration, 2.0));
+  // check for (numerical) errors
+  assert(r >= ron);
+  assert(r <= roff);
 
   setD (VSRC_1, VSRC_1, -r);
 }
@@ -192,6 +209,7 @@ PROP_OPT [] = {
   { "Roff", PROP_REAL, { 1e12, PROP_NO_STR }, PROP_POS_RANGE },
   { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
   { "MaxDuration", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_MIN_VAL (10*NR_TINY) },
+  { "Transition", PROP_STR, { PROP_NO_VAL, "spline" }, PROP_RNG_STR3 ("abrupt", "linear", "spline") },
   PROP_NO_PROP };
 struct define_t tswitch::cirdef =
   { "Switch", 2, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };

qucs-core/src/components/vrect.cpp

@@ -52,7 +52,11 @@ void vrect::initDC (void) {
   nr_double_t tf = getPropertyDouble ("Tf");
   if (tr > th) tr = th;
   if (tf > tl) tf = tl;
-  nr_double_t a  = (th + (tf - tr) / 2) / (th + tl);
+  // DC value defined as 0.0 instead of
+  // (th + (tf - tr) / 2) / (th + tl) previously used
+  // so that the transient starting value will also be 0, 
+  // otherwise a discontinuity occurs
+  nr_double_t a  = 0.0;
   nr_double_t u  = getPropertyDouble ("U") * a;
   allocMatrixMNA ();
   voltageSource (VSRC_1, NODE_1, NODE_2, u);

qucs/qucs/components/switch.cpp

@@ -35,6 +35,8 @@ Switch::Switch()
 		QObject::tr("simulation temperature in degree Celsius")));
   Props.append(new Property("MaxDuration", "1e-6", false,
 		QObject::tr("Max possible switch transition time (transition time 1/100 smallest value in 'time', or this number)")));
+  Props.append(new Property("Transition", "spline", false,
+		QObject::tr("Resistance transition shape")+" [abrupt, linear, spline]"));
 
   createSymbol();
   tx = x1+4;