Removed unnecessary function calls in BDF and Radau. Minor cleanup of…

… Radau as well.

examples/one_element_timoshenko_rod.py

@@ -17,8 +17,8 @@ def smoothstep2(x, x_min=0, x_max=1):
 
 # geometry of the rod
 length = 1 # [m]
-width = 0.01 # [m]
-# width = 1e-3 # [m]
+# width = 1e-2 # [m]
+width = 1e-3 # [m]
 density = 8.e3  # [kg / m^3]
 
 # material properties

pyproject.toml

@@ -1,6 +1,6 @@
 [build-system]
 build-backend = "mesonpy"
-requires = ["meson-python>=0.15.0"]
+requires = ["meson-python>=0.16.0"]
 
 [project]
 name = "scipy_dae"

scipy_dae/integrate/_dae/base.py

@@ -224,7 +224,6 @@ def fun(t, y, yp):
         self.fun = fun
         self.fun_single = fun_single
         self.fun_vectorized = fun_vectorized
-        self.f = self.fun(self.t, self.y, self.yp)
 
         self.direction = np.sign(t_bound - t0) if t_bound != t0 else 1
         self.status = 'running'
@@ -276,8 +275,11 @@ def _validate_jac(self, jac, sparsity):
             else:
                 sparsity_y, sparsity_yp = None, None
 
-            def jac_wrapped(t, y, yp, f):
+            def jac_wrapped(t, y, yp):
                 self.njev += 1
+
+                f = self.fun_single(t, y, 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)
@@ -295,8 +297,8 @@ def jac_wrapped(t, y, yp, f):
                     t, yp, f, threshold, self.jac_factor_yp, sparsity_yp)
 
                 # # test better Jacobian approximation
-                # method = "2-point"
-                # # method = "3-point"
+                # # 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)
@@ -306,7 +308,7 @@ def jac_wrapped(t, y, yp, f):
 
                 return Jy, Jyp
 
-            Jy, Jyp = jac_wrapped(t0, y0, yp0, self.f)
+            Jy, Jyp = jac_wrapped(t0, y0, yp0)
 
         elif callable(jac):
             Jy, Jyp = jac(t0, y0, yp0)

scipy_dae/integrate/_dae/bdf.py

@@ -264,7 +264,6 @@ def _step_impl(self):
         D = self.D
         y = self.y
         yp = self.yp
-        f = self.f
 
         max_step = self.max_step
         min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)
@@ -279,10 +278,6 @@ def _step_impl(self):
         else:
             h_abs = self.h_abs
 
-        # self.hs.append(h_abs)
-        # self.orders.append(self.order)
-        # print(f"- t: {t:.3e}; h: {h_abs:.3e}; order: {self.order}")
-
         atol = self.atol
         rtol = self.rtol
         order = self.order
@@ -331,7 +326,7 @@ def _step_impl(self):
                 if not converged:
                     if current_jac:
                         break
-                    Jy, Jyp = self.jac(t, y, yp, f)
+                    Jy, Jyp = self.jac(t, y, yp)
                     LU = None
                     current_jac = True
 
@@ -371,19 +366,14 @@ def _step_impl(self):
         self.Jyp = Jyp
         self.LU = LU
 
-        f_new = self.fun(t_new, y_new, yp_new)
-        self.f = f_new
-
         # Update differences. The principal relation here is
         # D^{j + 1} y_n = D^{j} y_n - D^{j} y_{n - 1}. Keep in mind that D
         # contained difference for previous interpolating polynomial and
         # d = D^{k + 1} y_n. Thus this elegant code follows.
-        # print(f"D before:\n{D}")
         D[order + 2] = d - D[order + 1]
         D[order + 1] = d
         for i in reversed(range(order + 1)):
             D[i] += D[i + 1]
-        # print(f"D after:\n{D}")
 
         if self.n_equal_steps < order + 1:
             return True, None
@@ -407,8 +397,8 @@ def _step_impl(self):
         # choose order with largest factor
         delta_order = np.argmax(factors) - 1
 
-        # choose with smallest error
-        # TODO: This choice is advertaised in Shampine2002 but experminets 
+        # choose order with smallest error
+        # TODO: This choice is advertised in Shampine2002 but experiments 
         # indicate it is not worth it
         # delta_order = np.argmin(error_norms) - 1