solve_ivpでのIntegration step failed.の解決法がわからない。

Question

https://qiita.com/trami/items/a1bdec427fc0420e6f19
このサイトを参考に、以下のような三重振り子のシミュレーションを作ったのですが、0.7秒までしか計算が行われません。
以下の運動方程式では、角速度の一階微分から、角速度、θ、を求めるというようにして角度を計算させようとしています。
しかし、どれだけ調べても原因を見つけることができません。何が原因なのかを教えていただけると幸いです。
どうかよろしくお願いいたします。

from numpy import sin,cos
from scipy.integrate import solve_ivp
import numpy as np
import csv
from fractions import Fraction as fr

Pi  =3.14     #円周率
G   =9.806    #重力加速度
f=open(result.csv,'a',newline="")
writer=csv.writer(f)
th1_DEG =90
th2_DEG =90
th3_DEG =90
om1_DEG =0
om2_DEG =0
om3_DEG =0           
L1      =1
L2      =1
L3      =1
M1      =1
M2      =1
M3      =1
th1 =th1_DEG*Pi/180     #初期の角度をラジアンに変換
th2 =th2_DEG*Pi/180
th3 =th3_DEG*Pi/180
om1 =om1_DEG*Pi/180
om2 =om2_DEG*Pi/180
om3 =om3_DEG*Pi/180

def func(t,state): #運動方程式
    dydx=np.zeros_like(state)
    dydx[0]=state[1]    #dydx[0]=om1
    delta1=state[0]-state[2]    #デルタ１
    delta2=state[2]-state[4]    #デルタ２
    delta3=state[4]-state[0]    #デルタ３
    dydx[1]=((M3*cos(delta2)-fr(1,4)*(M2+4*M3))*(state[3]*state[3]*L3*M3*sin(delta3)-state[3]*state[3]*(L2*(M2+2*M3)*
    sin(delta1)+2*L2*M3*sin(delta2)*cos(delta3))+state[1]*state[1]*2*L1*M3*sin(delta3)*cos(delta3)+(2*M3*cos(delta3)*
    sin(state[4])-(M1+2*M2+2*M3)*sin(state[0]))*G))/(L1*((-M3*sin(delta2)*sin(delta2)-fr(1,4)*M2)*(2*M3*sin(delta3)*
    sin(delta3)+fr(1,2)*M1+2*M2)+fr(1,2)*(M2*cos(delta1)+2*M3*cos(delta1)-2*M3*cos(delta2)*cos(delta3))*(M2*cos(delta1)+
    2*M3*cos(delta1)-2*M3*cos(delta2)*cos(delta3))))
    dydx[2]=state[3]    #dydx[2]=om2
    dydx[3]=(state[1]*state[1]*L1*(2*M3*M3*sin(delta3)*cos(delta2)*cos(delta3)*cos(delta3)+M3*(M2+2*M3)*sin(delta1)*
    cos(delta3)*cos(delta3)+fr(1,2)*M3*(M2+4*M3)*sin(delta3)*cos(delta3)-fr(1,2)*M3*(M1+4*M2+4*M3)-fr(1,4)*(M2+2*M3)*
    (M1+4*M2+4*M3)*sin(delta1)-2*M3*M3*sin(delta3)*cos(delta2)*cos(delta2)*cos(delta3))+state[3]*state[3]*L2*(fr(1,2)*
    M3*(M1+4*M2+4*M3)*sin(delta1)))/(L2*(fr(1,2)*M3*(M1+4*M2+4*M3)*cos(delta2)*cos(delta2)+fr(1,2)*M3*(M2+4*M3)*
    cos(delta3)*cos(delta3)+fr(1,4)*(M2+2*M3)*(M2+4*M3)*cos(delta1)+M3*M3*(cos(delta2)**3)*cos(delta3)-2*M3*M3*cos(delta2)*
    cos(delta2)*cos(delta3)*cos(delta3)-fr(1,8)*(M2+4*M3)*(M1+4*M2+4*M3)-M3*(M2+2*M3)*cos(delta1)*cos(delta2)*cos(delta2)-
    fr(1,4)*M3*(M2+4*M3)*cos(delta2)*cos(delta3)))
    dydx[4]=state[5]    #dydx[4]=om3
    dydx[5]=(state[1]*state[1]*((fr(1,2)*(M2+4*M3)*cos(delta3)-(M2+2*M3)*cos(delta1)*cos(delta2))*(fr(1,2)*L1*(M1+4*M2+4*
    M3)*sin(delta3))+((M2+2*M3)*cos(delta1)*cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*cos(delta2))*(M2+2*M3)*sin(delta2))-
    state[3]*state[3]*L2*((fr(1,2)*(M2+4*M3)*cos(delta3)-(M2+2*M3)*cos(delta1)*cos(delta2))*((M2+2*M3)*sin(delta1)*
    cos(delta3)+fr(1,2)*(M1+4*M2+4*M3)*sin(delta2))-((M2+2*M3)*cos(delta1)*cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*cos(delta2))
    *(M2+2*M3)*sin(delta2)*cos(delta1))+state[5]*state[5]*L3*(((M2+2*M3)*cos(delta1)*cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*
    cos(delta2))*M3*sin(delta2)*cos(delta3)+(fr(1,2)*(M2+4*M3)*cos(delta3)-(M2+2*M3)*cos(delta1)*cos(delta2))*M3*sin(delta3)*
    cos(delta3))-G*((fr(1,2)*(M2+4*M3)*cos(delta3)-(M2+2*M3)*cos(delta1)*cos(delta2))*((M1+2*M2+2*M3)*sin(state[0])*
    cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*sin(state[4]))-((M2+2*M3)*cos(delta1)*cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*cos(delta2))*
    (M2+2*M3)*cos(state[0])*sin(delta2)))/(L3*((fr(1,2)*(M2+4*M3)*cos(delta3)-(M2+2*M3)*cos(delta1)*cos(delta2))*(M3*
    cos(delta3)*cos(delta3)-fr(1,4)*(M1+4*M2+4*M3))-((M2+2*M3)*cos(delta1)*cos(delta3)-fr(1,2)*(M1+4*M2+4*M3)*
    cos(delta2))*(M3*cos(delta2)*cos(delta3)-fr(1,2)*(M2+2*M3)*cos(delta1))))

    return dydx
t_span  =[0.0,10]                                #時間生成
dt      =0.01
t       =np.arange(t_span[0],t_span[1],dt)
writer.writerow([float(dt),float(th1_DEG),float(th2_DEG),float(th3_DEG),float(om1_DEG),float(om2_DEG),float(om3_DEG),float(L1),float(L2),float(L3),float(M1),float(M2),float(M3)])

state   =[th1,om1,th2,om2,th3,om3]
print(len(t))
y   =solve_ivp(func,t_span,state,t_eval=t).y
for tt,th1,th2,th3 in zip(t,y[0,:],y[2,:],y[4,:]):
    x1=L1 *sin(th1)
    y1=-L1*cos(th1)
    x2=x1+L2*sin(th2)
    y2=y1-L2*cos(th2)
    x3=x2+L3*sin(th3)
    y3=y2-L3*sin(th3)
    X1= fr(1,2)*L1*sin(th1)
    Y1=-fr(1,2)*L1*cos(th1)
    X2=x1+fr(1,2)*L2*sin(th2)
    Y2=y1-fr(1,2)*L2*cos(th2)
    X3=x2+fr(1,2)*L3*sin(th3)
    Y3=y2-fr(1,2)*L3*cos(th3)
    Row_list=[float(tt),float(x1),float(y1),float(X1),float(Y1),float(x2)
    ,float(y2),float(X2),float(Y2),float(x3),float(y3),float(X3),float(Y3)]
    writer.writerow(Row_list)
f.close() 

そのあと、

y   =solve_ivp(func,t_span,state,t_eval=t).y

ここのコードを

sol =solve_ivp(func,t_span,state,t_eval=t)
y=sol.y
print(sol.status)

このように書き換えると、-1が出力されました。調べてみると、Integration step failedというエラーのようなのですが、いかにして解決するべきかがわかりません。

Answer 1

式が可積分ではないので積分できなくてエラーになっていたように思います。解決はできていませんが、運動方程式の部分でエラーが出ているということは分かりました。
ルンゲクッタ法で近似値を出そうと思います。