fixes and updated tests, bumping version to 0.4.2

.gitignore

@@ -46,6 +46,10 @@ nosetests.xml
 coverage.xml
 *.cover
 .hypothesis/
+test_reports/
+*.pyc
+__pycache__/
+.pytest_cache/
 
 # Jupyter Notebook
 .ipynb_checkpoints
@@ -60,6 +64,9 @@ env/
 venv/
 ENV/
 
+# Shell
+*.bat
+
 # IDEs
 .vscode/
 .idea/

examples/example_vanderpol.py

@@ -11,7 +11,7 @@
 
 from pathsim import Simulation, Connection
 from pathsim.blocks import Scope, ODE
-from pathsim.solvers import ESDIRK32, ESDIRK54
+from pathsim.solvers import ESDIRK32, ESDIRK43
 
 
 # VAN DER POL OSCILLATOR INITIAL VALUE PROBLEM ==========================================
@@ -42,7 +42,7 @@ def jac(x, u, t):
     ]
 
 #initialize simulation with the blocks, connections, timestep and logging enabled
-Sim = Simulation(blocks, connections, dt=0.5, log=True, Solver=ESDIRK54, tolerance_lte_abs=1e-6, tolerance_lte_rel=1e-4)
+Sim = Simulation(blocks, connections, dt=0.5, log=True, Solver=ESDIRK43, tolerance_lte_abs=1e-6, tolerance_lte_rel=1e-4)
 
 Sim.run(3*mu)
 

pathsim/solvers/_solver.py

@@ -124,18 +124,18 @@ def buffer(self):
         self.x_0 = self.x
 
 
-    def change(self, Solver, tolerance_lte=None):
+    def change(self, Sol, **solver_args):
         """
         Change the integration engine to a new type and initialize 
         with previous solver arguments so it can continue from where 
         the 'old' solver stopped.
         """
 
-        #check if new tolerance is defined
-        tol = self.tolerance_lte if tolerance_lte is None else tolerance_lte
-
         #create new engine from self
-        engine = Solver(self.initial_value, self.func, self.jac, tol)
+        engine = Sol(initial_value=self.initial_value, 
+                     func=self.func, 
+                     jac=self.jac, 
+                     **solver_args)
 
         #set internal state of new engine from self
         engine.set(self.get())

pathsim/solvers/dirk2.py

@@ -71,6 +71,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/dirk3.py

@@ -64,6 +64,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/esdirk32.py

@@ -110,6 +110,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/esdirk4.py

@@ -70,6 +70,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/esdirk43.py

@@ -118,6 +118,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/esdirk54.py

@@ -120,6 +120,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/esdirk85.py

@@ -177,6 +177,9 @@ def solve(self, u, t, dt):
         #use the jacobian
         if self.jac is not None:
 
+            #most recent butcher coefficient
+            b = self.BT[self.stage][self.stage]
+
             #compute jacobian of fixed-point equation
             jac_g = dt * b * self.jac(self.x, u, t)
 

pathsim/solvers/rkdp87.py

@@ -117,7 +117,7 @@ def step(self, u, t, dt):
         self.x = dt * sum(k*b for k, b in zip(self.Ks.values(), self.BT[self.stage])) + self.x_0
 
         #error and step size control
-        if self.stage < 8:
+        if self.stage < 12:
             self.stage += 1
             return True, 0.0, 0.0, 1.0
         else: 

pathsim/solvers/rkv65.py

@@ -60,8 +60,10 @@ def __init__(self,
                    7:[           11/144,      0,              0,          256/693,                0,   125/504, 125/528,        5/72], 
                    8:[           28/477,      0,              0,          212/441,   -312500/366177, 2125/1764,       0, -2105/35532, 2995/17766]}
 
-        #coefficients for truncation error
-        self.TR = [a-b for a, b in zip(self.BT[7], self.BT[8])]
+        #compute coefficients for truncation error
+        _A1 = [11/144, 0, 0, 256/693,              0,   125/504, 125/528,        5/72, 0]
+        _A2 = [28/477, 0, 0, 212/441, -312500/366177, 2125/1764,       0, -2105/35532, 2995/17766]        
+        self.TR = [a-b for a, b in zip(_A1, _A2)]
 
 
     def error_controller(self, dt):

setup.py

@@ -5,7 +5,7 @@
 
 setup(
     name="pathsim",
-    version="0.4.1",
+    version="0.4.2",
     author="Milan Rother",
     author_email="[email protected]",
     description="A block based time domain system simulation framework.",

tests/blocks/test_block.py

@@ -128,7 +128,7 @@ def test_step(self):
         B = Block()
 
         #test default implementation 
-        self.assertEqual(B.step(None, None), (True, 0.0, 1.0))
+        self.assertEqual(B.step(None, None), (True, 0.0, 0.0, 1.0))
 
 
 

tests/blocks/test_integrator.py

@@ -60,16 +60,18 @@ def test_set_solver(self):
         #test that no solver is initialized
         self.assertEqual(I.engine, None)
 
-        I.set_solver(Solver, tolerance_lte=1e-6)
+        I.set_solver(Solver, tolerance_lte_abs=1e-6, tolerance_lte_rel=1e-3)
 
         #test that solver is now available
         self.assertTrue(isinstance(I.engine, Solver))
-        self.assertEqual(I.engine.tolerance_lte, 1e-6)
+        self.assertEqual(I.engine.tolerance_lte_rel, 1e-3)
+        self.assertEqual(I.engine.tolerance_lte_abs, 1e-6)
 
-        I.set_solver(Solver, tolerance_lte=1e-3)
+        I.set_solver(Solver, tolerance_lte_abs=1e-4, tolerance_lte_rel=1e-2)
 
         #test that solver tolerance is changed
-        self.assertEqual(I.engine.tolerance_lte, 1e-3)
+        self.assertEqual(I.engine.tolerance_lte_rel, 1e-2)
+        self.assertEqual(I.engine.tolerance_lte_abs, 1e-4)
 
 
     def test_update(self):

tests/solvers/_referenceproblems.py

@@ -0,0 +1,40 @@
+
+########################################################################################
+##
+##                      REFERENCE PROBLEMS FOR SOLVER TESTING
+##
+##                              Milan Rother 2024
+##
+########################################################################################
+
+# IMPORTS ==============================================================================
+
+import numpy as np
+
+
+# TEST PROBLEMS ========================================================================
+
+class Problem:
+    def __init__(self, name, func, jac, x0, solution):
+        self.name = name
+        self.func = func
+        self.jac = jac
+        self.x0 = x0
+        self.solution = solution
+
+
+#create some reference problems for testing
+problems = [
+    Problem(name="linear_feedback", 
+            func=lambda x, u, t: -x, 
+            jac=lambda x, u, t: -1, 
+            x0=1.0, 
+            solution=lambda t: np.exp(-t)
+            ),
+    Problem(name="logistic", 
+            func=lambda x, u, t: x*(1-x), 
+            jac=lambda x, u, t: 1-2*x, 
+            x0=0.5, 
+            solution=lambda t: 1/(1 + np.exp(-t))
+            )
+]