Implementation of events with signature event(t, y, yp). (#38)

Author: JonasBreuling

.github/workflows/main.yml

@@ -24,7 +24,7 @@ jobs:
     - name: Setup python
       uses: actions/setup-python@v5
       with:
-        python-version: "3.8" 
+        python-version: "3.10" 
         cache: "pip" # caching pip dependencies
 
     - name: Install dependencies and scipy_dae

README.md

@@ -1,4 +1,4 @@
-# scipy_dae - solving differential algebraic equations (DAE's) and implicit differential equations (IDE's) in Python
+# scipy_dae - solving differential algebraic equations (DAEs) and implicit differential equations (IDEs) in Python
 
 <p align="center">
 <a href="https://github.com/JonasBreuling/scipy_dae/actions/workflows/main.yml/badge.svg"><img alt="Actions Status" src="https://github.com/JonasBreuling/scipy_dae/actions/workflows/main.yml/badge.svg"></a>
@@ -9,7 +9,7 @@
 <a href="https://pypi.org/project/scipy_dae/"><img alt="PyPI" src="https://img.shields.io/pypi/v/scipy_dae"></a>
 </p>
 
-Python implementation of solvers for differential algebraic equations (DAE's) and implicit differential equations (IDE's) that should be added to scipy one day.
+Python implementation of solvers for differential algebraic equations (DAEs) and implicit differential equations (IDEs) that should be added to scipy one day.
 
 Currently, two different methods are implemented.
 
@@ -111,11 +111,11 @@ plt.show()
 
 More examples are given in the [examples](examples/) directory, which includes
 
-* ordinary differential equations (ODE's)
+* ordinary differential equations (ODEs)
     * [Van der Pol oscillator](examples/odes/van_der_pol.py)
     * [Sparse brusselator](examples/odes/sparse_brusselator.py)
     * [Stiff SE(3) Cosserat rod](examples/odes/se3_cosserat_rod.py)
-* differential algebraic equations (DAE's)
+* differential algebraic equations (DAEs)
     * [Robertson problem (index 1)](examples/daes/robertson.py)
     * [Akzo Nobel problem (index 1)](examples/daes/akzo_nobel.py)
     * [Stiff transistor amplifier (index 1)](examples/daes/stiff_transistor_amplifier.py)
@@ -126,7 +126,7 @@ More examples are given in the [examples](examples/) directory, which includes
     * [Cartesian pendulum (index 3)](examples/daes/pendulum.py)
     * [Particle on circular track (index 3)](examples/daes/arevalo.py)
     * [Andrews' squeezer mechanism (index 3)](examples/daes/andrews.py)
-* implicit differential equations (IDE's)
+* implicit differential equations (IDEs)
     * [Weissinger's implicit equation](examples/ides/weissinger.py)
     * [Problem I.542 of E. Kamke](examples/ides/kamke.py)
     * [Jackiewicz' implicit equation](examples/ides/jackiewicz.py)
@@ -177,7 +177,7 @@ $$
 \end{aligned}
 $$
 
-Since the implemented solvers are designed for index 1 DAE's we have to perform some sort of index reduction. Therefore, we transform the semi-explicit form into a general form as proposed by [Gear](https://doi.org/10.1137/0909004). The resulting index 1 system is given as
+Since the implemented solvers are designed for index 1 DAEs we have to perform some sort of index reduction. Therefore, we transform the semi-explicit form into a general form as proposed by [Gear](https://doi.org/10.1137/0909004). The resulting index 1 system is given as
 
 $$
 \begin{aligned}
@@ -226,7 +226,7 @@ $$
 \end{aligned}
 $$
 
-Since the implemented solvers are designed for index 1 DAE's we have to perform some sort of index reduction. Therefore, we use the [stabilized index 1 formulation of Hiller and Anantharaman](https://doi.org/10.1002/nme.1620320803). The resulting index 1 system is given as
+Since the implemented solvers are designed for index 1 DAEs we have to perform some sort of index reduction. Therefore, we use the [stabilized index 1 formulation of Hiller and Anantharaman](https://doi.org/10.1002/nme.1620320803). The resulting index 1 system is given as
 
 $$
 \begin{aligned}

pyproject.toml

@@ -4,16 +4,16 @@ requires = ["meson-python>=0.16.0"]
 
 [project]
 name = "scipy_dae"
-version = "0.0.6"
+version = "0.0.8"
 authors = [
   { name="Jonas Breuling", email="jonas.breuling@inm.uni-stuttgart.de" },
 ]
-description = "Python implementation of solvers for differential algebraic equations (DAE's) and implicit differential equations (IDE's) that should be added to scipy one day."
+description = "Python implementation of solvers for differential algebraic equations (DAEs) and implicit differential equations (IDEs) that should be added to scipy one day."
 readme = "README.md"
 requires-python = ">=3.8"
 dependencies = [
-  "numpy>=1.24.4",
-  "scipy>=1.10.1",
+  "numpy>=2.0.1",
+  "scipy>=1.14.0",
   "matplotlib>=3.4.3",
 ]
 classifiers = [

scipy_dae/integrate/_dae/base.py

@@ -46,9 +46,7 @@ def fun_wrapped(t, y, yp):
 
     return fun_wrapped, y0, yp0
 
-# TODO:
-# - check documentation
-# - add consistent initial conditions somehow
+
 class DaeSolver:
     """Base class for DAE solvers.
 
@@ -102,17 +100,11 @@ class DaeSolver:
     y0 : array_like, shape (n,)
         Initial state.
     yp0 : array_like, shape (n,), optional
-        Initial derivative. If the derivative is not given by the user, it is 
-        esimated using??? => TODO: See Petzold/ Shampine and coworkers how 
-        this is done.
+        Initial derivative. It it is assumed that the value is consistent, 
+        i.e. f(t0, y0, yp0) = 0.
     t_bound : float
         Boundary time --- the integration won't continue beyond it. It also
         determines the direction of the integration.
-    var_index : array_like, shape (n,), optional
-        Differentiation index of the respective row of f(t, y, y') = 0. 
-        Depending on this index, the error estimates are scaled by the 
-        stepsize h**(index - 1) in order to ensure convergence.
-        Default is None, which means all equations are differential equations.
     vectorized : bool
         Whether `fun` can be called in a vectorized fashion. Default is False.
 
@@ -166,7 +158,7 @@ class DaeSolver:
     TOO_SMALL_STEP = "Required step size is less than spacing between numbers."
 
     def __init__(self, fun, t0, y0, yp0, t_bound, rtol, atol, 
-                 first_step=None, max_step=np.inf, vectorized=False,
+                 first_step=None, max_step=np.inf, vectorized=False, 
                  jac=None, jac_sparsity=None, support_complex=False):
         self.t_old = None
         self.t = t0
@@ -281,7 +273,7 @@ def jac_wrapped(t, y, yp):
                 #     lambda t, yp: self.fun_vectorized(t, y, yp), 
                 #     t, yp, f, self.atol, self.jac_factor_yp, sparsity_yp)
 
-                # TODO: This choice is better but not optimal!
+                # note: this choice is better but not optimal
                 threshold = self.atol / self.rtol
                 Jy, self.jac_factor_y = num_jac(
                     lambda t, y: self.fun_vectorized(t, y, yp), 
@@ -423,6 +415,7 @@ def _step_impl(self):
     def _dense_output_impl(self):
         raise NotImplementedError
 
+
 class DAEDenseOutput:
     """Base class for local interpolant over step made by an ODE solver.
 

scipy_dae/integrate/_dae/bdf.py

@@ -152,17 +152,13 @@ class BDFDAE(DaeSolver):
 
         It is generally recommended to provide the Jacobians rather than
         relying on a finite-difference approximation.
-    # TODO: Adapt and test this.
     jac_sparsity : {None, array_like, sparse matrix}, optional
         Defines a sparsity structure of the Jacobian matrix for a
         finite-difference approximation. Its shape must be (n, n). This argument
         is ignored if `jac` is not `None`. If the Jacobian has only few non-zero
         elements in *each* row, providing the sparsity structure will greatly
         speed up the computations [4]_. A zero entry means that a corresponding
         element in the Jacobian is always zero. If None (default), the Jacobian
-    # t_eval = np.concatenate((t, t + t1)) / 2
-    t_eval = np.array([t0, *(t0 + 1e-3 + np.cumsum(np.diff(t)))])
-    # t_eval = np.concatenate((t, t + t1, t + 2 * t1)) / 3
         is assumed to be dense.
     vectorized : bool, optional
         Whether `fun` can be called in a vectorized fashion. Default is False.
@@ -231,8 +227,9 @@ def __init__(self, fun, t0, y0, yp0, t_bound, max_step=np.inf,
                  NDF_strategy="stability", **extraneous):
         warn_extraneous(extraneous)
         super().__init__(fun, t0, y0, yp0, t_bound, rtol, atol, 
-                         first_step, max_step, vectorized, jac, 
-                         jac_sparsity, support_complex=True)
+                         first_step=first_step, max_step=max_step, 
+                         vectorized=vectorized, jac=jac, 
+                         jac_sparsity=jac_sparsity, support_complex=True)
 
         self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))
 

scipy_dae/integrate/_dae/common.py

@@ -1,6 +1,6 @@
 from itertools import groupby
 import numpy as np
-# TODO: use sparse QR when available in scipy
+# note: use sparse QR when available in scipy
 from scipy.linalg import qr, solve_triangular
 from scipy.integrate._ivp.common import norm, EPS
 from scipy.optimize._numdiff import approx_derivative