Minor cleanup.

Co-authored-by: Copilot <copilot@github.com>

Author: Jonas Breuling

solve_dae/integrate/_dae/base.py

@@ -180,7 +180,7 @@ def __init__(self, fun, t0, y0, yp0, t_bound, rtol, atol,
         self.nfev = 0
         self.njev = 0
         self.nlu = 0
-        self._nlusove = 0
+        self._nlusolve = 0
 
         if vectorized:
             def fun_single(t, y, yp):
@@ -224,20 +224,17 @@ def lu(A):
                 return splu(A)
 
             def solve_lu(LU, b):
-                self._nlusove += 1
+                self._nlusolve += 1
                 return LU.solve(b)
-
-            self.I = eye(self.n, format='csc')
         else:
             def lu(A):
                 self.nlu += 1
                 return lu_factor(A, overwrite_a=True)
 
             def solve_lu(LU, b):
-                self._nlusove += 1
+                self._nlusolve += 1
                 return lu_solve(LU, b, overwrite_b=True)
 
-            self.I = np.identity(self.n)
         self.lu = lu
         self.solve_lu = solve_lu
 
@@ -402,12 +399,7 @@ def dense_output(self):
         if self.t_old is None:
             raise RuntimeError("Dense output is available after a successful "
                                "step was made.")
-
-        if self.n == 0 or self.t == self.t_old:
-            # Handle corner cases of empty solver and no integration.
-            return ConstantDAEDenseOutput(self.t_old, self.t, self.y, self.yp)
-        else:
-            return self._dense_output_impl()
+        return self._dense_output_impl()
 
     def _step_impl(self):
         raise NotImplementedError
@@ -455,25 +447,3 @@ def __call__(self, t):
 
     def _call_impl(self, t):
         raise NotImplementedError
-
-
-class ConstantDAEDenseOutput(DAEDenseOutput):
-    """Constant value interpolator.
-
-    This class used for degenerate integration cases: equal integration limits
-    or a system with 0 equations.
-    """
-    def __init__(self, t_old, t, value, derivative):
-        super().__init__(t_old, t)
-        self.value = value
-        self.derivative = derivative
-
-    def _call_impl(self, t):
-        if t.ndim == 0:
-            return self.value, self.derivative
-        else:
-            ret = np.empty((self.value.shape[0], t.shape[0]))
-            ret[:] = self.value[:, None]
-            ret_der = np.empty((self.derivative.shape[0], t.shape[0]))
-            ret_der[:] = self.derivative[:, None]
-            return ret, ret_der

solve_dae/integrate/_dae/bdf.py

@@ -104,7 +104,7 @@ class BDFDAE(DaeSolver):
     max_order : int, optional
         Highest order of the method with 1 <= max_order <= 6, although 
         max_order = 6 should be used with care due to the limited stability 
-        propoerties of the corresponding BDF method.
+        properties of the corresponding BDF method.
     NDF_strategy : string, optional
         The strategy that is applied for obtaining numerical differentiation 
         formulas (NDF):
@@ -235,7 +235,7 @@ def __init__(self, fun, t0, y0, yp0, t_bound, max_step=np.inf,
 
         assert 1 <= max_order <= 6, "Ensure that 1 <= max_order <= 6."
         if max_order == 6:
-            warn("Choosing `max_order = 6` is not recomended due to its poor stability properties.",
+            warn("Choosing `max_order = 6` is not recommended due to its poor stability properties.",
                  stacklevel=3)
         self.max_order = max_order
 
@@ -565,3 +565,58 @@ def _call_impl(self, t):
         #     yp3 = np.squeeze(yp3)
 
         # return y3, yp3
+
+
+# class BdfDenseOutput(DAEDenseOutput):
+#     """Newton backward-difference interpolant for the BDF solution.
+
+#     Both y(t) and yp(t) are evaluated in a single O(order) Horner pass
+#     using dual-number differentiation.
+
+#     The interpolant P(t) satisfies P(t_n - j·h) = y_{n-j} for j = 0..order-1.
+#     P'(t) is the analytical time-derivative of P.
+#     """
+
+#     def __init__(self, t_old, t, h, order, D):
+#         super().__init__(t_old, t)
+#         self.order = order
+#         self.h = h          # signed step size: h = direction * h_abs
+#         self.t_n = t        # current time = right end of the interpolation window
+#         self.D = D          # shape (order+1, n): divided differences
+
+#     def _call_impl(self, t):
+#         """Evaluate (y(t), yp(t)) using Horner's rule with dual-number derivative.
+
+#         The Newton form of P(t) is:
+#             P = D[0] + x_1*(D[1] + x_2*(D[2] + ... + x_q*D[q]))
+#         where x_j = (t - (t_n - (j-1)*h)) / (j*h).
+
+#         Carrying dual variables (v, dv/dt) through the Horner recurrence gives
+#         both P(t) and P'(t) in a single O(order) pass — no second traversal needed.
+#         """
+#         scalar = t.ndim == 0
+#         t = np.atleast_1d(t)           # (n_pts,)
+#         n_pts = len(t)
+#         n     = self.D.shape[1]        # state dimension
+
+#         # Initialise with the innermost term D[order] (constant in t)
+#         # v : (n, n_pts) — running polynomial value
+#         # dv: (n, n_pts) — running derivative d/dt of v
+#         v  = np.tile(self.D[self.order, :, None], (1, n_pts))   # (n, n_pts)
+#         dv = np.zeros_like(v)
+
+#         for j in range(self.order, 0, -1):
+#             # Node for this Horner step
+#             node = self.t_n - (j - 1) * self.h   # scalar
+#             x_j  = (t - node) / (j * self.h)     # (n_pts,)
+#             dx_j = 1.0 / (j * self.h)            # scalar (same for all pts)
+
+#             # Dual-number multiply: new = D[j-1] + x_j * v
+#             # d(new)/dt = dx_j * v + x_j * dv
+#             new_v  = self.D[j - 1, :, None] + x_j[None, :] * v
+#             new_dv = dx_j * v + x_j[None, :] * dv
+#             v, dv = new_v, new_dv
+
+#         if scalar:
+#             return np.squeeze(v), np.squeeze(dv)
+#         return v, dv

solve_dae/integrate/_dae/radau.py

@@ -288,7 +288,7 @@ class RadauDAE(DaeSolver):
         bounded and determined solely by the solver.
     rtol, atol : float and array_like, optional
         Relative and absolute tolerances. The solver keeps the local error
-        estimates less than ``atol + rtol * abs(y)``. HHere `rtol` controls a
+        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
         relative accuracy (number of correct digits), while `atol` controls
         absolute accuracy (number of correct decimal places). To achieve the
         desired `rtol`, set `atol` to be smaller than the smallest value that
@@ -413,9 +413,7 @@ def __init__(self, fun, t0, y0, yp0, t_bound, stages=3,
         ) = radau_constants(stages)
 
         self.h_abs_old = None
-        self.h_abs_oldold = None
         self.error_norm_old = None
-        self.error_norm_oldold = None
 
         # modify tolerances as in radau.f line 824ff and 920ff
         # TODO: Document this rescaling
@@ -460,7 +458,6 @@ def __init__(self, fun, t0, y0, yp0, t_bound, stages=3,
         self.current_jac = True
         self.LU_real = None
         self.LU_complex = None
-        self.Z = None
 
     def _step_impl(self):
         t = self.t
@@ -643,7 +640,6 @@ def _step_impl(self):
         self.y = y_new
         self.yp = yp_new
 
-        self.Z = Z
         self.Y = Y
         self.Yp = Yp
 

solve_dae/integrate/_dae/tests/test_integration_robertson.py

@@ -34,7 +34,7 @@ def F_robertson(t, state, statep):
         ):
             with suppress_warnings() as sup:
                 sup.filter(UserWarning,
-                        "Choosing `max_order = 6` is not recomended due to its "
+                        "Choosing `max_order = 6` is not recommended due to its "
                         "poor stability properties.")
                 res = solve_dae(F_robertson, tspan, y0, yp0, rtol=rtol,
                                 atol=atol, method=method, max_order=max_order,
@@ -93,7 +93,7 @@ def F_robertson(t, y, yp):
         ):
             with suppress_warnings() as sup:
                 sup.filter(UserWarning,
-                        "Choosing `max_order = 6` is not recomended due to its "
+                        "Choosing `max_order = 6` is not recommended due to its "
                         "poor stability properties.")
                 res = solve_dae(F_robertson, tspan, y0, yp0, rtol=rtol,
                                 atol=atol, method=method, max_order=max_order,