Stiff noonlinear Timoshenko rod and stiff van der Pol equation.

examples/one_element_timoshenko_rod.py

@@ -0,0 +1,252 @@
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy_dae.integrate import solve_dae, consistent_initial_conditions
+from scipy.integrate import solve_ivp
+
+
+def smoothstep2(x, x_min=0, x_max=1):
+    x = np.clip((x - x_min) / (x_max - x_min), 0, 1)
+    return 6 * x**5 - 15 * x**4 + 10 * x**3
+
+
+EPS = 1e-6
+
+t0 = 0
+t1 = 10
+
+# geometry of the rod
+length = 1 # [m]
+width = 0.01 # [m]
+# width = 1e-3 # [m]
+density = 8.e3  # [kg / m^3]
+
+# material properties
+E = 260.0e9 # [N / m^2]
+G = 100.0e9 # [N / m^2]
+# E = 260.0e6 # [N / m^2]
+# G = 100.0e6 # [N / m^2]
+
+
+# # [m] -> [mm]
+# length *= 1e3
+# # cross section properties
+# width *= 1e3
+# density *= 1e-9
+
+# # [s] -> [ms]
+# t1 *= 1e3
+
+# # [N / m^2] -> [kN / mm^2]
+# E *= 1e-9
+# G *= 1e-9
+
+A = width**2
+I = width**4 / 12
+
+EA = E * A
+GA = G * A
+
+M = np.diag([A, A, I]) * density * length**3 / 3
+
+
+def A_IB(phi):
+    sphi, cphi = np.sin(phi), np.cos(phi)
+    return np.array([
+        [cphi, -sphi],
+        [sphi,  cphi],
+    ])
+
+def df_int(xi, vq):
+    x, y, phi = vq
+
+    A_IB_ = A_IB(phi * xi)
+    n = A_IB_ @ np.diag([E * A, G * A]) @ (
+        A_IB_.T @ vq[:2] / length - np.array([1, 0])
+    )
+    m = E * I * phi / length
+
+    K1 = np.array([
+        [1,  0, 0],
+        [0,  1, 0],
+        [y, -x, 1]
+    ])
+
+    return K1 @ np.array([*n, m])
+
+
+def f_int(vq):
+    # one point quadrature
+    x = np.array([0])
+    w = np.array([2])
+
+    # # two point quadrature
+    # x = np.array([-np.sqrt(1 / 3), np.sqrt(1 / 3)])
+    # w = np.array([1, 1])
+
+    # transform from [-1, 1] on [0, 1]
+    b = 1
+    a = 0
+    w = (b - a) / 2 * w
+    x = (a + b) / 2 + (b - a) / 2 * x
+
+    return np.sum([df_int(xi, vq) * wi for (xi, wi) in zip(x, w)], axis=0)
+
+    # # exact internal forces
+    # x, y, phi = vq
+
+    # K1 = np.array([
+    #     [1,  0, 0],
+    #     [0,  1, 0],
+    #     [y, -x, 1]
+    # ])
+
+    # sphi, cphi = np.sin(phi), np.cos(phi)
+
+    # if abs(phi) > EPS:
+    #     cphi2 = 0.5 + sphi * cphi / (2 * phi)
+    #     sphi2 = 0.5 - np.sin(2 * phi) / (4 * phi)
+    #     sphi_cphi = sphi**2 / (2 * phi)
+
+    #     n0 = EA * np.array([
+    #         sphi,
+    #         1 - cphi
+    #     ]) / phi
+    # else:
+    #     n0 = EA * np.array([
+    #         # https://www.wolframalpha.com/input?i=taylor+series+sin%28x%29+%2F+x
+    #         1 - phi**2 / 6 + phi**4 / 120,
+    #         # https://www.wolframalpha.com/input?i=taylor+series+%281+-+cos%28x%29%29+%2F+x
+    #         phi / 2 - phi**3 / 24 + phi**5 / 720
+    #     ])
+
+    # K_phi = np.array([
+    #     [EA * cphi**2 + GA * sphi**2,     (EA - GA) * cphi * sphi],
+    #     [    (EA - GA) * cphi * sphi, EA * sphi**2 + GA * cphi**2],
+    # ])
+
+    # n = K_phi @ vq[:2] - n0
+    # m = E * I * phi / length
+
+    # return K1 @ np.array([*n, m])
+
+
+def m(t):
+    m_max = E * I / length * 2
+    return m_max * (
+        smoothstep2(t, 0, 4)
+        - smoothstep2(t, 4, 4.1)
+    )
+    # if t < 4:
+    #     return t * m_max
+    # else:
+    #     return 0.0
+
+def f_ext(t, vq):
+    return np.array([
+        0,
+        0,
+        m(t)
+    ])
+
+
+def F(t, vy, vyp):
+    vq, vu = vy[:3], vy[3:]
+    vqp, vup = vyp[:3], vyp[3:]
+
+    return np.concatenate([
+        vqp - vu,
+        M @ vup + f_int(vq) - f_ext(t, vq),
+    ])
+
+def f(t, vy):
+    vq, vu = vy[:3], vy[3:]
+
+    q_dot = vu
+    u_dot = -np.linalg.solve(M, f_int(vq) - f_ext(t, vq))
+
+    return np.concatenate([q_dot, u_dot])
+
+
+if __name__ == "__main__":
+    # time span
+    t_span = (t0, t1)
+    t_eval = np.linspace(t0, t1, num=int(1e4))
+    t_eval = None
+
+    # method = "BDF"
+    method = "Radau"
+    # method = "RK23"
+
+    # initial positions
+    q0 = np.array([length, 0, 0], dtype=float)
+
+    # r_OP1 = np.array([0.76738219, 0.54021608]) * 1e3
+    # phi1 = 1.22664684
+    # q0 = np.array([*r_OP1, phi1])
+
+    # initial velocities
+    u0 = np.array([0, 0, 0], dtype=float)
+
+    y0 = np.concatenate([q0, u0])
+
+    q_dot0 = u0
+    u_dot0 = np.array([0, 0, 0], dtype=float)
+    yp0 = np.concatenate([q_dot0, u_dot0])
+
+    F0 = F(t0, y0, yp0)
+    print(f"F0: {F0}")
+
+    yp0 = np.zeros_like(y0)
+    print(f"y0: {y0}")
+    print(f"yp0: {yp0}")
+    y0, yp0, F0 = consistent_initial_conditions(F, t0, y0, yp0)
+    print(f"y0: {y0}")
+    print(f"yp0: {yp0}")
+    print(f"F0: {F0}")
+
+    # solver options
+    atol = rtol = 1e-3
+
+    ##############
+    # dae solution
+    ##############
+    start = time.time()
+    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
+    y = sol.y
+    success = sol.success
+    status = sol.status
+    message = sol.message
+    print(f"message: {message}")
+    print(f"elapsed time: {end - start}")
+    print(f"nfev: {sol.nfev}")
+    print(f"njev: {sol.njev}")
+    print(f"nlu: {sol.nlu}")
+
+    # visualization
+    fig, ax = plt.subplots(4, 1)
+
+    ax[0].set_xlabel("t")
+    ax[0].set_xlabel("x")
+    ax[0].grid()
+    ax[0].plot(t, y[0], "-ok")
+
+    ax[1].set_xlabel("t")
+    ax[1].set_xlabel("y")
+    ax[1].grid()
+    ax[1].plot(t, y[1], "-ok")
+
+    ax[2].set_xlabel("t")
+    ax[2].set_xlabel("phi")
+    ax[2].grid()
+    ax[2].plot(t, y[2], "-ok")
+
+    ax[3].set_xlabel("t")
+    ax[3].set_xlabel("h")
+    ax[3].grid()
+    ax[3].plot(t[1:], np.diff(t), "-ok")
+
+    plt.show()

examples/van_der_pol.py

@@ -5,20 +5,23 @@
 from scipy_dae.integrate import solve_dae, consistent_initial_conditions
 
 
-"""Stiff van der Pol equation, see mathworks.
-References:
------------
-mathworks: https://de.mathworks.com/help/matlab/math/solve-stiff-odes.html
+"""Stiff van der Pol equation, see [1]_.
+
+References
+----------
+.. [1] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
+       Stiff and Differential-Algebraic Problems", Sec. IV.8.
 """
 
-mu = 1e3
+
+eps = 1e-6
 
 def rhs(t, y):
     y1, y2 = y
 
     yp = np.zeros(2, dtype=y.dtype)
     yp[0] = y2
-    yp[1] = mu * (1 - y1 * y1) * y2 - y1
+    yp[1] = ((1 - y1 * y1) * y2 - y1) / eps
 
     return yp
 
@@ -34,14 +37,14 @@ def f(t, z):
 if __name__ == "__main__":
     # time span
     t0 = 0
-    t1 = 3e3
+    t1 = 2
     t_span = (t0, t1)
 
     # method = "BDF"
     method = "Radau"
 
     # initial conditions
-    y0 = np.array([2, 0], dtype=float)
+    y0 = np.array([2, -0.66], dtype=float)
     yp0 = rhs(t0, y0)
     z0 = np.concatenate((y0, yp0))
 
@@ -54,11 +57,11 @@ def f(t, z):
     print(f"fnorm: {fnorm}")
 
     # solver options
-    atol = rtol = 1e-4
+    atol = rtol = 7e-5
 
     t_eval = np.linspace(t0, t1, num=int(1e3))
     t_eval = None
-    first_step = 1e-3
+    first_step = 1e-6
 
     ####################
     # reference solution
@@ -102,7 +105,7 @@ def f(t, z):
     # visualization
     fig, ax = plt.subplots(4, 1)
 
-    t_eval = np.linspace(t0, t1, num=int(1e2))
+    t_eval = np.linspace(t0, t1, num=int(1e3))
     y_eval = sol.sol(t_eval)
 
     ax[0].plot(t, y[0], "ok", label=f"y ({method})", mfc="none")

scipy_dae/integrate/_dae/radau.py

@@ -7,10 +7,9 @@
 
 
 NEWTON_MAXITER = 6  # Maximum number of Newton iterations.
-NEWTON_MAXITER_EMBEDDED = 2  # Maximum number of Newton iterations for embedded method.
-# DAMPING_RATIO_ERROR_ESTIMATE = 1.0
-# DAMPING_RATIO_ERROR_ESTIMATE = 0.8
-DAMPING_RATIO_ERROR_ESTIMATE = -0.8
+NEWTON_MAXITER_EMBEDDED = 1  # Maximum number of Newton iterations for embedded method.
+DAMPING_RATIO_ERROR_ESTIMATE = -0.8 # 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.
 
@@ -658,6 +657,8 @@ def _step_impl(self):
         factor = predict_factor(h_abs, h_abs_old, error_norm, error_norm_old, s)
         factor = min(MAX_FACTOR, safety * factor)
 
+        # TODO: That's what Hairer proposes. We should also use the method of??? here that ensures a smooth change of the sep size, see ???.
+        # if not recompute_jac and 1.0 <= factor <= 1.2:
         if not recompute_jac and factor < 1.2:
             factor = 1
         else: