Old dense output again.

examples/van_der_pol.py

@@ -52,8 +52,8 @@ def fun_composite(t, z):
     t1 = 3e3
     t_span = (t0, t1)
 
-    # method = "BDF"
-    method = "Radau"
+    method = "BDF"
+    # method = "Radau"
 
     # initial conditions
     y0 = np.array([2, 0], dtype=float)
@@ -105,7 +105,7 @@ def fun_composite(t, z):
     print(f"elapsed time: {end - start}")
     t = sol.t
     y = sol.y
-    yp = sol.yp
+    # yp = sol.yp
     success = sol.success
     status = sol.status
     message = sol.message
@@ -131,12 +131,12 @@ def fun_composite(t, z):
 
     yp_scipy = np.array([rhs(ti, yi) for ti, yi in zip(t_scipy, y_scipy.T)]).T
 
-    ax[2].plot(t, yp[0], "-ok", label=f"yp0 ({method})", mfc="none")
+    # ax[2].plot(t, yp[0], "-ok", label=f"yp0 ({method})", mfc="none")
     ax[2].plot(t_scipy, yp_scipy[0], "-xr", label="yp0 scipy")
     ax[2].legend()
     ax[2].grid()
 
-    ax[3].plot(t, yp[1], "-ok", label=f"yp1 ({method})", mfc="none")
+    # ax[3].plot(t, yp[1], "-ok", label=f"yp1 ({method})", mfc="none")
     ax[3].plot(t_scipy, yp_scipy[1], "-xr", label="yp1 scipy")
     ax[3].legend()
     ax[3].grid()

scipy_dae/integrate/_dae/base.py

@@ -278,22 +278,22 @@ def _validate_jac(self, jac, sparsity):
 
             def jac_wrapped(t, y, yp, f):
                 self.njev += 1
-                # Jy, self.jac_factor_y = num_jac(
-                #     lambda t, y: self.fun_vectorized(t, y, yp), 
-                #     t, y, f, self.atol, self.jac_factor_y, sparsity_y)
-                # Jyp, self.jac_factor_yp = num_jac(
-                #     lambda t, yp: self.fun_vectorized(t, y, yp), 
-                #     t, yp, f, self.atol, self.jac_factor_yp, sparsity_yp)
+                Jy, self.jac_factor_y = num_jac(
+                    lambda t, y: self.fun_vectorized(t, y, yp), 
+                    t, y, f, self.atol, self.jac_factor_y, sparsity_y)
+                Jyp, self.jac_factor_yp = num_jac(
+                    lambda t, yp: self.fun_vectorized(t, y, yp), 
+                    t, yp, f, self.atol, self.jac_factor_yp, sparsity_yp)
 
-                # test better Jacobian approximation
-                method = "2-point"
-                # method = "3-point"
-                Jy = approx_derivative(
-                    lambda y: self.fun_single(t, y, yp), 
-                    y, f0=f, sparsity=sparsity_y, method=method)
-                Jyp = approx_derivative(
-                    lambda yp: self.fun_single(t, y, yp), 
-                    yp, f0=f, sparsity=sparsity_yp, method=method)
+                # # test better Jacobian approximation
+                # method = "2-point"
+                # # method = "3-point"
+                # Jy = approx_derivative(
+                #     lambda y: self.fun_single(t, y, yp), 
+                #     y, f0=f, sparsity=sparsity_y, method=method)
+                # Jyp = approx_derivative(
+                #     lambda yp: self.fun_single(t, y, yp), 
+                #     yp, f0=f, sparsity=sparsity_yp, method=method)
 
                 return Jy, Jyp
 

scipy_dae/integrate/_dae/bdf.py

@@ -1,9 +1,9 @@
 import numpy as np
 from warnings import warn
 from scipy.integrate._ivp.common import norm, EPS, warn_extraneous
-# from scipy.integrate._ivp.base import DenseOutput
+from scipy.integrate._ivp.base import DenseOutput
+# from .base import DAEDenseOutput as DenseOutput
 from .dae import DaeSolver
-from .base import DAEDenseOutput as DenseOutput
 
 
 NEWTON_MAXITER = 4
@@ -440,4 +440,5 @@ def _call_impl(self, t):
         else:
             y += self.D[0, :, None]
 
-        return y, np.zeros_like(y)
+        return y
+        # return y, np.zeros_like(y)

scipy_dae/integrate/_dae/dae.py

@@ -412,9 +412,9 @@ def fun(t, x, x_dot, fun=fun):
                 if sol is None:
                     sol = solver.dense_output()
                 ts.append(t_eval_step)
-                # ys.append(sol(t_eval_step))
-                y, yp = sol(t_eval_step)
-                ys.append(y)
+                ys.append(sol(t_eval_step))
+                # y, yp = sol(t_eval_step)
+                # ys.append(y)
                 yps.append(yp)
                 t_eval_i = t_eval_i_new
 
@@ -424,7 +424,7 @@ def fun(t, x, x_dot, fun=fun):
     message = MESSAGES.get(status, message)
 
     if t_events is not None:
-        raise NotImplementedError("Events are not ready yet")
+        # raise NotImplementedError("Events are not ready yet")
         t_events = [np.asarray(te) for te in t_events]
         y_events = [np.asarray(ye) for ye in y_events]
 
@@ -440,17 +440,17 @@ def fun(t, x, x_dot, fun=fun):
     if dense_output:
         if t_eval is None:
             sol = OdeSolution(
-                ts, interpolants, #alt_segment=False
-                # ts, interpolants, alt_segment=True if method in [BDFDAE] else False
+                # ts, interpolants, #alt_segment=False
+                ts, interpolants, alt_segment=True if method in [BDFDAE] else False
             )
         else:
             sol = OdeSolution(
-                ti, interpolants, #alt_segment=False
-                # ti, interpolants, alt_segment=True if method in [BDFDAE] else False
+                # ti, interpolants, #alt_segment=False
+                ti, interpolants, alt_segment=True if method in [BDFDAE] else False
             )
     else:
         sol = None
 
     return OdeResult(t=ts, y=ys, yp=yps, sol=sol, t_events=t_events, y_events=y_events,
                      nfev=solver.nfev, njev=solver.njev, nlu=solver.nlu,
-                     status=status, message=message, success=status >= 0)
+                     status=status, message=message, success=status>=0)

scipy_dae/integrate/_dae/radau.py

@@ -3,8 +3,8 @@
 from scipy.linalg import eig, cdf2rdf
 from scipy.integrate._ivp.common import norm, EPS, warn_extraneous
 from scipy.integrate._ivp.base import DenseOutput
+# from .base import DAEDenseOutput as DenseOutput
 from .dae import DaeSolver
-from .base import DAEDenseOutput as DenseOutput
 
 
 def butcher_tableau(s):
@@ -559,8 +559,8 @@ def _step_impl(self):
             if self.sol is None:
                 Z0 = np.zeros((s, y.shape[0]))
             else:
-                # Z0 = self.sol(t + h * C).T - y
-                Z0 = self.sol(t + h * C)[0].T - y
+                Z0 = self.sol(t + h * C).T - y
+                # Z0 = self.sol(t + h * C)[0].T - y
 
             scale = atol + np.abs(y) * rtol
 
@@ -658,15 +658,15 @@ def eval_collocation_polynomial(xi, C, Y):
             p0_2 = coeffs[0] # TODO: That is cheap if V_inv is already available :)
             p0_2 = V_inv[0] @ Y # and that is even cheaper ;)
 
-            # p0_2 = eval_collocation_polynomial2(0.0, C, Y)
-            # p0_2 = compute_interp(C, Y, 0.0)
-            # # p0_2 = compute_interp(C, Y, 0.0, df=Yp)
-            # p0_ = eval_collocation_polynomial2(0.0, C, Y)
-            p0_ = eval_collocation_polynomial(0.0, C, Y)
-            # p1_ = eval_collocation_polynomial(1, C, Y)
-            # print(f"p0_: {p0_}")
-            # print(f"p0_2: {p0_2}")
-            # error_collocation = y - p0_
+            # # p0_2 = eval_collocation_polynomial2(0.0, C, Y)
+            # # p0_2 = compute_interp(C, Y, 0.0)
+            # # # p0_2 = compute_interp(C, Y, 0.0, df=Yp)
+            # # p0_ = eval_collocation_polynomial2(0.0, C, Y)
+            # p0_ = eval_collocation_polynomial(0.0, C, Y)
+            # # p1_ = eval_collocation_polynomial(1, C, Y)
+            # # print(f"p0_: {p0_}")
+            # # print(f"p0_2: {p0_2}")
+            # # error_collocation = y - p0_
             error_collocation = y - p0_2
 
             # c1, c2, c3 = C
@@ -691,7 +691,7 @@ def eval_collocation_polynomial(xi, C, Y):
             yp_hat_new = MU_REAL * (v @ Yp - b0 * yp)
             F = self.fun(t_new, y_new, yp_hat_new)
             # TODO: Is this already l-stable or do we have to add a possible second 
-            # multiplication with LU_real?
+            # multiplication with LUyp_real?
             error_Fabien = self.solve_lu(LU_real, -F)
 
             # # add another Newton step for stabilization
@@ -791,24 +791,24 @@ def eval_collocation_polynomial(xi, C, Y):
 
     def _compute_dense_output(self):
         Q = np.dot(self.Z.T, P)
-        Yp = (A_inv / (self.h_abs * self.direction)) @ self.Z
-        Qp = np.dot(Yp.T, P)
-        # return RadauDenseOutput(self.t_old, self.t, self.y_old, Q)
-        return RadauDenseOutput(self.t_old, self.t, self.y_old, Q, self.yp_old, Qp)
+        # Yp = (A_inv / (self.h_abs * self.direction)) @ self.Z
+        # Qp = np.dot(Yp.T, P)
+        return RadauDenseOutput(self.t_old, self.t, self.y_old, Q)
+        # return RadauDenseOutput(self.t_old, self.t, self.y_old, Q, self.yp_old, Qp)
 
     def _dense_output_impl(self):
         return self.sol
 
 class RadauDenseOutput(DenseOutput):
-    # def __init__(self, t_old, t, y_old, Q):
-    def __init__(self, t_old, t, y_old, Q, yp_old, Qp):
+    def __init__(self, t_old, t, y_old, Q):
+    # def __init__(self, t_old, t, y_old, Q, yp_old, Qp):
         super().__init__(t_old, t)
         self.h = t - t_old
         self.Q = Q
-        self.Qp = Qp
+        # self.Qp = Qp
         self.order = Q.shape[1] - 1
         self.y_old = y_old
-        self.yp_old = yp_old
+        # self.yp_old = yp_old
 
     def _call_impl(self, t):
         x = (t - self.t_old) / self.h
@@ -820,12 +820,13 @@ def _call_impl(self, t):
             p = np.cumprod(p, axis=0)
         # Here we don't multiply by h, not a mistake.
         y = np.dot(self.Q, p)
-        yp = np.dot(self.Qp, p)
+        # yp = np.dot(self.Qp, p)
         if y.ndim == 2:
             y += self.y_old[:, None]
-            yp += self.yp_old[:, None]
+            # yp += self.yp_old[:, None]
         else:
             y += self.y_old
-            yp += self.yp_old
+            # yp += self.yp_old
 
-        return y, yp
+        return y
+        # return y, yp