Up to date benchmarks again.

examples/daes/brenan.py

@@ -47,7 +47,6 @@ def F(t, y, yp):
     start = time.time()
     # method = "BDF"
     method = "Radau"
-    # method = "PSIDE"
     sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method)
     end = time.time()
     print(f"elapsed time: {end - start}")

examples/odes/one_element_timoshenko_rod.py

@@ -212,7 +212,7 @@ def f(t, vy):
     # dae solution
     ##############
     start = time.time()
-    sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval, stages=3)
+    sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval, stages=5)
     # sol = solve_ivp(f, t_span, y0, method=method, t_eval=t_eval, atol=atol, rtol=rtol)
     end = time.time()
     t = sol.t

scipy_dae/integrate/_dae/benchmarks/knife_edge/knife_edge.py

@@ -81,7 +81,6 @@ def sol_true(t):
 
 def run_knife_edge():
     # exponents
-    # m_max = 
     m_max = 32
     ms = np.arange(m_max + 1)
 

scipy_dae/integrate/_dae/benchmarks/kvaerno/kvaerno.py

@@ -32,8 +32,7 @@ def true_sol(t):
 
 def run_kvaerno():
     # exponents
-    # m_max = 32
-    m_max = 10
+    m_max = 32
     ms = np.arange(m_max + 1)
 
     # tolerances and initial step size

scipy_dae/integrate/_dae/benchmarks/weissinger/weissinger.py

@@ -29,7 +29,7 @@ def run_weissinger():
     # tolerances and initial step size
     rtols = 10**(-(4 + ms / 4))
     atols = rtols
-    h0s = 1e-2 * rtols
+    h0s = 1e2 * rtols
 
     # time span
     t0 = np.sqrt(0.5)

scipy_dae/integrate/_dae/radau.py

@@ -6,10 +6,8 @@
 from .dae import DaeSolver
 
 
-DAMPING_RATIO_ERROR_ESTIMATE = 0.01 # Hairer (8.19) is obtained by the choice 1.0. 
+DAMPING_RATIO_ERROR_ESTIMATE = 0.05 # Hairer (8.19) is obtained by the choice 1.0. 
                                     # de Swart proposes 0.067 for s=3.
-MIN_FACTOR = 0.2  # Minimum allowed decrease in a step size.
-MAX_FACTOR = 10  # Maximum allowed increase in a step size.
 KAPPA = 1 # Factor of the smooth limiter
 
 # UNKNOWN_VELOCITIES = False
@@ -205,16 +203,11 @@ def solve_collocation_system(fun, t, y, h, Z0, scale, tol,
         if dW_norm_old is not None:
             rate = dW_norm / dW_norm_old
 
-        if (rate is not None and (rate >= 1 or rate ** (newton_max_iter - k) / (1 - rate) * dW_norm > tol)):
-            break
-
-        # if (rate is not None and rate >= 1.0):
-        #     # print(f"rate >= 1")
+        # if (rate is not None and (rate >= 1 or rate ** (newton_max_iter - k) / (1 - rate) * dW_norm > tol)):
         #     break
-        #     # if n_bad_iter > 5:
-        #     #     break
-        #     # else:
-        #     #     n_bad_iter += 1
+
+        if (rate is not None and rate >= 1.0):
+            break
 
         # # TODO: Why this is a bad indicator for divergence of the iteration?
         # if (rate is not None and rate ** (newton_max_iter - k) / (1 - rate) * dW_norm > tol):
@@ -230,20 +223,6 @@ def solve_collocation_system(fun, t, y, h, Z0, scale, tol,
             converged = True
             break
 
-        # # inexact simplified Newton method (https://doi.org/10.1137/S0036142999360573)
-        # # kappa2 = 0.1
-        # kappa2 = 0.01
-        # if (
-        #     dW_norm == 0 
-        #     or rate is not None and (
-        #         rate / (1 - rate) * dW_norm < tol 
-        #         # or dW_norm_old is not None and rate / (1 - rate) * dW_norm < kappa2 * rate**2 * dW_norm_old
-        #         or dW_norm_old is not None and rate / (1 - rate) * dW_norm < kappa2 * rate**2 / (1 - rate) * dW_norm_old
-        #     )
-        # ):
-        #     converged = True
-        #     break
-
         dW_norm_old = dW_norm
 
     return converged, k + 1, Y, Yp, Z, rate
@@ -299,12 +278,12 @@ def solve_collocation_system2(fun, t, y, h, Yp0, scale, tol,
         if dY_norm_old is not None:
             rate = dY_norm / dY_norm_old
 
-        if (rate is not None and (rate >= 1 or rate ** (newton_max_iter - k) / (1 - rate) * dY_norm > tol)):
-            break
-
-        # if (rate is not None and rate >= 1.0):
+        # if (rate is not None and (rate >= 1 or rate ** (newton_max_iter - k) / (1 - rate) * dY_norm > tol)):
         #     break
 
+        if (rate is not None and rate >= 1.0):
+            break
+
         # # TODO: Why this is a bad indicator for divergence of the iteration?
         # if (rate is not None and rate ** (newton_max_iter - k) / (1 - rate) * dY_norm > tol):
         #     break
@@ -357,9 +336,6 @@ def predict_factor(h_abs, h_abs_old, error_norm, error_norm_old, s):
     with np.errstate(divide='ignore'):
         factor = min(1, multiplier) * error_norm ** (-1 / (s + 1))
 
-    # # nonsmooth limiter
-    # factor = max(MIN_FACTOR, min(factor, MAX_FACTOR))
-
     # smooth limiter
     factor = 1 + KAPPA * np.arctan((factor - 1) / KAPPA)
 
@@ -656,8 +632,12 @@ def _step_impl(self):
                     Z0 = np.zeros((s, y.shape[0]))
             else:
                 if UNKNOWN_VELOCITIES:
-                    # note: this only works when we exrapolate the derivative of the collocation polynomial but do not use the sth order collocation polynomial for the derivatives as well.
+                    # note: this only works when we exrapolate the derivative 
+                    # of the collocation polynomial but do not use the sth order 
+                    # collocation polynomial for the derivatives as well.
                     Yp0 = self.sol(t + h * C)[1].T
+                    # Z0 = self.sol(t + h * C)[0].T - y
+                    # Yp0 = (1 / h) * A_inv @ Z0
                 else:
                     Z0 = self.sol(t + h * C)[0].T - y
             scale = atol + np.abs(y) * rtol
@@ -813,13 +793,15 @@ def _step_impl(self):
         return step_accepted, message
 
     def _compute_dense_output(self):
-        # Q = np.dot(self.Z.T, self.P)
-        # h = self.t - self.t_old
-        # Yp = (self.A_inv / h) @ self.Z
-        Z = self.Y - self.y_old
-        Q = np.dot(Z.T, self.P)
-        Zp = self.Yp - self.yp_old
+        Q = np.dot(self.Z.T, self.P)
+        h = self.t - self.t_old
+        Yp = (self.A_inv / h) @ self.Z
+        Zp = Yp - self.yp_old
         Qp = np.dot(Zp.T, self.P)
+        # Z = self.Y - self.y_old
+        # Q = np.dot(Z.T, self.P)
+        # Zp = self.Yp - self.yp_old
+        # Qp = np.dot(Zp.T, self.P)
         return RadauDenseOutput(self.t_old, self.t, self.y_old, Q, self.yp_old, Qp)
 
     def _dense_output_impl(self):