handle input containing division by zero better

.github/workflows/ode-toolbox-build.yml

@@ -43,7 +43,7 @@ jobs:
       fail-fast: false
       matrix:
         with_gsl: ["0", "1"]
-        sympy_version: ["==1.4", ""]  # empty string for "latest"
+        sympy_version: [""]  # empty string for "latest"
 
     steps:
       - name: Checkout ODE-toolbox code

odetoolbox/shapes.py

@@ -29,7 +29,7 @@
 import sympy
 import sympy.parsing.sympy_parser
 
-from sympy.core.expr import Expr as SympyExpr   # works for both sympy 1.4 and 1.8
+from sympy.core.expr import Expr as SympyExpr
 
 from .config import Config
 from .sympy_helpers import _check_numerical_issue, _check_forbidden_name, _custom_simplify_expr, _is_constant_term, _is_sympy_type, _is_zero

odetoolbox/system_of_shapes.py

@@ -195,17 +195,14 @@ def get_sub_system(self, symbols):
 
 
     def _generate_propagator_matrix(self, A):
-        r"""Generate the propagator matrix by matrix exponentiation.
-
-        XXX: the default custom simplification expression does not work well with sympy 1.4 here. Consider replacing sympy.simplify() with _custom_simplify_expr() if sympy 1.4 support is dropped.
-        """
+        r"""Generate the propagator matrix by matrix exponentiation."""
 
         # naive: calculate propagators in one step
-        # P_naive = sympy.simplify(sympy.exp(A * sympy.Symbol(Config().output_timestep_symbol)))
+        # P_naive = _custom_simplify_expr(sympy.exp(A * sympy.Symbol(Config().output_timestep_symbol)))
 
         # optimized: be explicit about block diagonal elements; much faster!
         blocks = get_block_diagonal_blocks(np.array(A))
-        propagators = [sympy.simplify(sympy.exp(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol))) for block in blocks]
+        propagators = [_custom_simplify_expr(sympy.exp(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol))) for block in blocks]
         P = sympy.Matrix(scipy.linalg.block_diag(*propagators))
 
         # check the result