Hermite interpolation works for Radau but seems to be worse for extra…

…polation.

examples/particle_on_circular_track.py

@@ -75,8 +75,8 @@ def sol_true(t):
     t_eval = np.linspace(t0, t1, num=int(5e2))
     # t_eval = None
 
-    method = "BDF"
-    # method = "Radau"
+    # method = "BDF"
+    method = "Radau"
 
     # initial conditions
     y0, yp0 = sol_true(t0)

examples/robertson.py

@@ -56,7 +56,8 @@ def fun_composite(t, z):
     t0 = 0
     t1 = 1e7
     t_span = (t0, t1)
-    t_eval = np.logspace(-6, 7, num=1000)
+    # t_eval = np.logspace(-6, 7, num=1000)
+    t_eval = np.logspace(-6, 7, num=500)
 
     # method = "BDF"
     method = "Radau"
@@ -85,6 +86,7 @@ def fun_composite(t, z):
     print(f"elapsed time: {end - start}")
     t_scipy = sol.t
     y_scipy = sol.y
+    yp_scipy = np.array([f(ti, yi) for ti, yi in zip(t_scipy, y_scipy.T)]).T
     success = sol.success
     status = sol.status
     message = sol.message
@@ -104,6 +106,7 @@ def fun_composite(t, z):
     print(f"elapsed time: {end - start}")
     t = sol.t
     y = sol.y
+    yp = sol.yp
     success = sol.success
     status = sol.status
     message = sol.message
@@ -115,17 +118,28 @@ def fun_composite(t, z):
     print(f"nlu: {sol.nlu}")
 
     # visualization
-    fig, ax = plt.subplots()
-
-    ax.set_xlabel("t")
-    ax.plot(t, y[0], "-ok", label="y1 DAE" + f" ({method})", mfc="none")
-    ax.plot(t, y[1] * 1e4, "-ob", label="y2 DAE" + f" ({method})", mfc="none")
-    ax.plot(t, y[2], "-og", label="y3 DAE" + f" ({method})", mfc="none")
-    ax.plot(t_scipy, y_scipy[0], "xr", label="y1 scipy Radau", markersize=7)
-    ax.plot(t_scipy, y_scipy[1] * 1e4, "xy", label="y2 scipy Radau", markersize=7)
-    ax.plot(t_scipy, y_scipy[2], "xm", label="y3 scipy Radau", markersize=7)
-    ax.set_xscale("log")
-    ax.legend()
-    ax.grid()
+    fig, ax = plt.subplots(2, 1)
+
+    ax[0].plot(t, y[0], "-ok", label="y1 DAE" + f" ({method})", mfc="none")
+    ax[0].plot(t, y[1] * 1e4, "-ob", label="y2 DAE" + f" ({method})", mfc="none")
+    ax[0].plot(t, y[2], "-og", label="y3 DAE" + f" ({method})", mfc="none")
+    ax[0].plot(t_scipy, y_scipy[0], "xr", label="y1 scipy" + f" ({method})", markersize=7)
+    ax[0].plot(t_scipy, y_scipy[1] * 1e4, "xy", label="y2 scipy" + f" ({method})", markersize=7)
+    ax[0].plot(t_scipy, y_scipy[2], "xm", label="y3 scipy" + f" ({method})", markersize=7)
+    ax[0].set_xlabel("t")
+    ax[0].set_xscale("log")
+    ax[0].legend()
+    ax[0].grid()
+
+    ax[1].plot(t, yp[0], "-ok", label="yp1 DAE" + f" ({method})", mfc="none")
+    ax[1].plot(t, yp[1], "-ob", label="yp2 DAE" + f" ({method})", mfc="none")
+    ax[1].plot(t, yp[2], "-og", label="yp3 DAE" + f" ({method})", mfc="none")
+    ax[1].plot(t_scipy, yp_scipy[0], "xr", label="yp1 scipy" + f" ({method})", markersize=7)
+    ax[1].plot(t_scipy, yp_scipy[1], "xy", label="yp2 scipy" + f" ({method})", markersize=7)
+    ax[1].plot(t_scipy, yp_scipy[2], "xm", label="yp3 scipy" + f" ({method})", markersize=7)
+    ax[1].set_xlabel("t")
+    ax[1].set_xscale("log")
+    ax[1].legend()
+    ax[1].grid()
 
     plt.show()

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)
@@ -69,10 +69,10 @@ def fun_composite(t, z):
     # print(f"fnorm: {fnorm}")
 
     # solver options
-    atol = rtol = 1e-4
+    atol = rtol = 1e-6
 
     t_eval = np.linspace(t0, t1, num=int(1e3))
-    t_eval = None
+    # t_eval = None
     first_step = 1e-3
 
     ####################

scipy_dae/integrate/_dae/bdf.py

@@ -485,6 +485,8 @@ def _call_impl(self, t):
                 y2 += self.D[j, :, None] * factor_y2 / (factorial(j) * self.h**j)
                 yp2 += self.D[j, :, None] * factor_yp2 / (factorial(j) * self.h**j)
 
+        assert np.allclose(y, y2)
+
         # compute y recursively which enables the computation of yp as well using Horner's rule, see
         # https://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture08.pdf
         # https://en.wikipedia.org/wiki/Horner's_method#Python_implementation
@@ -501,8 +503,6 @@ def _call_impl(self, t):
         # print(f"|y - y2|: {np.linalg.norm(y - y2)}")
         # # print(f"|yp - yp2|: {np.linalg.norm(yp - yp2)}")
 
-        assert np.allclose(y, y2)
-
         # print(f"|y - y3|: {np.linalg.norm(y - y3)}")
         print(f"|yp2 - yp3|: {np.linalg.norm(yp2 - yp3)}")
         assert np.allclose(y, y3)

scipy_dae/integrate/_dae/radau.py

@@ -393,7 +393,7 @@ class RadauDAE(DaeSolver):
            solution of implicit delay differential equations", Journal of 
            Computational and Applied Mathematics 185, 261-277, 2006.
     """
-    def __init__(self, fun, t0, y0, yp0, t_bound, stages=3,
+    def __init__(self, fun, t0, y0, yp0, t_bound, stages=5,
                  max_step=np.inf, rtol=1e-3, atol=1e-6, 
                  continuous_error_weight=0.0, jac=None, 
                  jac_sparsity=None, vectorized=False, 
@@ -638,6 +638,42 @@ def _step_impl(self):
 
         return step_accepted, message
 
+    # def _compute_dense_output(self):
+    #     h = self.t - self.t_old
+    #     Y = self.y_old + self.Z
+    #     Yp = (self.A_inv / h) @ self.Z
+    #     # Zp = Yp - self.yp_old
+
+    #     # Compute the inverse of the Vandermonde matrix to get the 
+    #     # interpolation matrix P.
+    #     c_hat = np.array([0, *self.C])
+    #     c_hat = self.C
+    #     N = len(c_hat)
+    #     vander = np.vander(c_hat, N=2 * N, increasing=True)#.T
+    #     # P1 = vander[2:, 2:]
+    #     # P2 = vander[1:-1, 1:-1]
+
+    #     V = np.zeros((2 * N, 2 * N))
+    #     V[:N, :] = vander
+    #     V[N:, 1:] = vander[:, :-1] * (np.arange(2 * N - 1) + 1) #/ h
+
+    #     P = np.linalg.inv(V)
+
+    #     # rhs = np.zeros((2 * N, Y.shape[1]))
+    #     # rhs[0] = self.y_old
+    #     # rhs[1:N] = Y
+    #     # rhs[N] = self.yp_old * h
+    #     # rhs[N+1:] = Yp * h
+
+    #     # only use stage values
+    #     rhs = np.zeros((2 * N, Y.shape[1]))
+    #     rhs[:N] = Y
+    #     rhs[N:] = Yp * h
+
+    #     Q = P @ rhs
+
+    #     return RadauDenseOutputHermite(self.t_old, self.t, Q)
+
     def _compute_dense_output(self):
         Q = np.dot(self.Z.T, self.P)
         h = self.t - self.t_old
@@ -674,4 +710,29 @@ def _call_impl(self, t):
             y = np.squeeze(y)
             yp = np.squeeze(yp)
 
+        return y, yp
+
+
+class RadauDenseOutputHermite(DenseOutput):
+    def __init__(self, t_old, t, Q):
+        super().__init__(t_old, t)
+        self.h = t - t_old
+        self.Q = Q
+        self.order = Q.shape[0] - 1
+
+    def _call_impl(self, t):
+        x = (t - self.t_old) / self.h
+        x = np.atleast_1d(x)
+
+        p = x**np.arange(self.order + 1)[:, None]
+        dp = np.zeros_like(p)
+        c = np.arange(1, self.order + 1)[:, None]
+        dp[1:] = (c / self.h) * (x**(c - 1))
+
+        y = np.dot(self.Q.T, p)
+        yp = np.dot(self.Q.T, dp)
+        if t.ndim == 0:
+            y = np.squeeze(y)
+            yp = np.squeeze(yp)
+
         return y, yp