improved error contron for GEAR methods

examples/example_nested_subsystems.py

@@ -86,7 +86,7 @@
     ]
 
 #initialize simulation with the blocks, connections, timestep and logging enabled
-Sim = Simulation(blocks, connections, dt=dt, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-5, tolerance_lte_rel=1e-3)
+Sim = Simulation(blocks, connections, dt=dt, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-5, tolerance_lte_rel=1e-4)
 
 #run simulation for some number of seconds
 Sim.run(3*mu)

examples/example_pendulum.py

@@ -19,7 +19,7 @@
     Scope
     )
 
-from pathsim.solvers import RKCK54, GEAR52A
+from pathsim.solvers import RKCK54
 
 
 # MATHEMATICAL PENDULUM =================================================================

examples/example_stickslip.py

@@ -46,7 +46,7 @@ def _f(x, u, t):
     ]
 
 #initialize simulation with the blocks, connections, timestep and logging enabled
-Sim = Simulation(blocks, connections, dt=0.01, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-6, tolerance_lte_rel=1e-3)
+Sim = Simulation(blocks, connections, dt=0.01, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-5, tolerance_lte_rel=1e-4)
 
 #run the simulation for some time
 Sim.run(2*T)

examples/example_vanderpol.py

@@ -42,7 +42,7 @@ def jac(x, u, t):
     ]
 
 #initialize simulation with the blocks, connections, timestep and logging enabled
-Sim = Simulation(blocks, connections, dt=5, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-5, tolerance_lte_rel=1e-3)
+Sim = Simulation(blocks, connections, dt=1, log=True, Solver=GEAR52A, tolerance_lte_abs=1e-6, tolerance_lte_rel=1e-4)
 
 Sim.run(3*mu)
 

pathsim/solvers/gear.py

@@ -19,8 +19,7 @@
 def compute_bdf_coefficients(order, timesteps):
     """
     Computes the coefficients for backward differentiation formulas for a given order.
-    The timesteps can be specified if the coefficients are for variable timestep BDF 
-    methods. 
+    The timesteps can be specified for variable timestep BDF methods. 
 
     For m-th order BDF we have for the n-th timestep:
         sum(alpha_i * x_i; i=n-m,...,n) = h_n * f_n(x_n, t_n)
@@ -190,8 +189,8 @@ def error_controller(self, x_m):
         #compute timestep scale factor using accuracy order of truncation error
         timestep_rescale = self.beta / error_norm ** (1/self.n)
 
-        #clip the rescale factor to a reasonable range
-        timestep_rescale = np.clip(timestep_rescale, 0.1, 10.0)
+        # #clip the rescale factor to a reasonable range
+        # timestep_rescale = np.clip(timestep_rescale, 0.1, 10.0)
 
         return success, error_norm, timestep_rescale
 
@@ -391,36 +390,23 @@ def error_controller(self, x_m, x_p):
 
         #compute the error norm and clip it
         error_norm_m = np.clip(float(np.max(scaled_error_m)), 1e-18, None)
-        error_norm_p = np.clip(float(np.max(scaled_error_p)), 1e-18, None)        
+        error_norm_p = np.clip(float(np.max(scaled_error_p)), 1e-18, None)      
 
-        #decrease the order if smaller order is more accurate
-        if error_norm_m < error_norm_p:
-
-            #success metric
-            success = error_norm_m <= 1.0
+        #success metric (use lower order estimate)
+        success = error_norm_m <= 1.0
 
-            #compute timestep scale factor using accuracy order of truncation error
-            timestep_rescale = self.beta / error_norm_m ** (1/self.n)
-            timestep_rescale = np.clip(timestep_rescale, 0.1, 10.0)
+        #compute timestep scale factor using accuracy order of truncation error
+        timestep_rescale = self.beta / error_norm_m ** (1/self.n)  
 
-            #decrease method order by one
+        #decrease the order if smaller order is more accurate (stability)
+        if error_norm_m < error_norm_p:
             self.n = max(self.n-1, self.n_min)
-
-            return success, error_norm_m, timestep_rescale
-
+
+        #increase the order if larger order is more accurate (accuracy -> larger steps)
         else:
-
-            #success metric
-            success = error_norm_p <= 1.0
-
-            #compute timestep scale factor using accuracy order of truncation error
-            timestep_rescale = self.beta / error_norm_p ** (1/(self.n + 1))
-            timestep_rescale = np.clip(timestep_rescale, 0.1, 10.0)
-
-            #increase method order by one
             self.n = min(self.n+1, self.n_max)
 
-            return success, error_norm_p, timestep_rescale
+        return success, error_norm_p, timestep_rescale
 
 
     # methods for timestepping ---------------------------------------------------------