Starting from version 2.7, customisation of spectral radius (h o ∞ ho_\infty h o ∞ γ \gamma γ
The algorithm is known to be able to conserve energy and momentum.
References:
The implementation follows the original algorithm. The only difference is the step size defined in suanPan is the mean step size of two sub-step sizes, that is Δ t / 2 \Delta{}t/2 Δ t /2
For further discussions on this type of algorithm, one can check the following references.
Syntax
Copy integrator BatheTwoStep (1) [2] [3]
# (1) int, unique integrator tag
# [2] double, spectral radius, \rho_\infty, default: 0
# [2] double, sub-step size, gamma, default: 0.5 Using integrator BatheTwoStep (1) with two optional parameters omitted gives the same results as in versions prior to version 2.7.
The spectral radius h o ∞ ho_\infty h o ∞ inclusive .
The sub-step size γ \gamma γ exclusive .
Theory
The First Sub-step
For trapezoidal rule,
v n + 1 = v n + γ Δ t 2 ( a n + a n + 1 ) , v_{n+1}=v_n+\dfrac{\gamma\Delta{}t}{2}\left(a_n+a_{n+1}\right), v n + 1 = v n + 2 γ Δ t ( a n + a n + 1 ) , u n + 1 = u n + γ Δ t 2 ( v n + v n + 1 ) . u_{n+1}=u_n+\dfrac{\gamma\Delta{}t}{2}\left(v_n+v_{n+1}\right). u n + 1 = u n + 2 γ Δ t ( v n + v n + 1 ) . Then,
u n + 1 = u n + γ Δ t 2 ( v n + v n + γ Δ t 2 ( a n + a n + 1 ) ) , u_{n+1}=u_n+\dfrac{\gamma\Delta{}t}{2}\left(v_n+v_n+\dfrac{\gamma\Delta{}t}{2}\left(a_n+a_{n+1}\right)\right), u n + 1 = u n + 2 γ Δ t ( v n + v n + 2 γ Δ t ( a n + a n + 1 ) ) , Δ u = γ Δ t v n + γ 2 Δ t 2 4 a n + γ 2 Δ t 2 4 a n + 1 . \Delta{}u=\gamma\Delta{}tv_n+\dfrac{\gamma^2\Delta{}t^2}{4}a_n+\dfrac{\gamma^2\Delta{}t^2}{4}a_{n+1}. Δ u = γ Δ t v n + 4 γ 2 Δ t 2 a n + 4 γ 2 Δ t 2 a n + 1 . One could obtain
a n + 1 = 4 γ 2 Δ t 2 Δ u − 4 γ Δ t v n − a n , Δ a = 4 γ 2 Δ t 2 Δ u − 4 γ Δ t v n − 2 a n , a_{n+1}=\dfrac{4}{\gamma^2\Delta{}t^2}\Delta{}u-\dfrac{4}{\gamma\Delta{}t}v_n-a_n,\qquad \Delta{}a=\dfrac{4}{\gamma^2\Delta{}t^2}\Delta{}u-\dfrac{4}{\gamma\Delta{}t}v_n-2a_n, a n + 1 = γ 2 Δ t 2 4 Δ u − γ Δ t 4 v n − a n , Δ a = γ 2 Δ t 2 4 Δ u − γ Δ t 4 v n − 2 a n , v n + 1 = 2 γ Δ t Δ u − v n , Δ v = 2 γ Δ t Δ u − 2 v n . v_{n+1}=\dfrac{2}{\gamma\Delta{}t}\Delta{}u-v_n,\qquad \Delta{}v=\dfrac{2}{\gamma\Delta{}t}\Delta{}u-2v_n. v n + 1 = γ Δ t 2 Δ u − v n , Δ v = γ Δ t 2 Δ u − 2 v n . The effective stiffness is then
K ˉ = K + 2 γ Δ t C + 4 γ 2 Δ t 2 M . \bar{K}=K+\dfrac{2}{\gamma\Delta{}t}C+\dfrac{4}{\gamma^2\Delta{}t^2}M. K ˉ = K + γ Δ t 2 C + γ 2 Δ t 2 4 M . The Second Sub-step
The second step is computed by
v n + 2 = v n + Δ t ( q 0 a n + q 1 a n + 1 + q 2 a n + 2 ) , v_{n+2}=v_n+\Delta{}t\left(q_0a_n+q_1a_{n+1}+q_2a_{n+2}\right), v n + 2 = v n + Δ t ( q 0 a n + q 1 a n + 1 + q 2 a n + 2 ) , u n + 2 = u n + Δ t ( q 0 v n + q 1 v n + 1 + q 2 v n + 2 ) . u_{n+2}=u_n+\Delta{}t\left(q_0v_n+q_1v_{n+1}+q_2v_{n+2}\right). u n + 2 = u n + Δ t ( q 0 v n + q 1 v n + 1 + q 2 v n + 2 ) . The parameters satisfy q 0 + q 1 + q 2 = 1 q_0+q_1+q_2=1 q 0 + q 1 + q 2 = 1
q 1 = ρ ∞ + 1 2 γ ( ρ ∞ − 1 ) + 4 , q 2 = 0.5 − γ q 1 . q_1=\dfrac{\rho_\infty+1}{2\gamma\left(\rho_\infty-1\right)+4},\qquad q_2=0.5-\gamma{}q_1. q 1 = 2 γ ( ρ ∞ − 1 ) + 4 ρ ∞ + 1 , q 2 = 0.5 − γ q 1 . Hence,
a n + 2 = v n + 2 − v n q 2 Δ t − q 0 q 2 a n − q 1 q 2 a n + 1 , a_{n+2}=\dfrac{v_{n+2}-v_n}{q_2\Delta{}t}-\dfrac{q_0}{q_2}a_n-\dfrac{q_1}{q_2}a_{n+1}, a n + 2 = q 2 Δ t v n + 2 − v n − q 2 q 0 a n − q 2 q 1 a n + 1 , v n + 2 = u n + 2 − u n q 2 Δ t − q 0 q 2 v n − q 1 q 2 v n + 1 . v_{n+2}=\dfrac{u_{n+2}-u_n}{q_2\Delta{}t}-\dfrac{q_0}{q_2}v_n-\dfrac{q_1}{q_2}v_{n+1}. v n + 2 = q 2 Δ t u n + 2 − u n − q 2 q 0 v n − q 2 q 1 v n + 1 . The effective stiffness is then
K ˉ = K + 1 q 2 Δ t C + 1 q 2 2 Δ t 2 M . \bar{K}=K+\dfrac{1}{q_2\Delta{}t}C+\dfrac{1}{q_2^2\Delta{}t^2}M. K ˉ = K + q 2 Δ t 1 C + q 2 2 Δ t 2 1 M . 