Home>

### python - i don't know how to solve integration step failed with solve_ivp

https://qiita.com/trami/items/a1bdec427fc0420e6f19
Using this site as a reference, I made a simulation of a triple pendulum as shown below, but it can only be calculated up to 0.7 seconds.
In the following equation of motion, the angle is calculated by obtaining the angular velocity, θ, from the first derivative of the angular velocity.
However, no matter how much you investigate, you cannot find the cause. I would appreciate it if you could tell me what caused it.

``````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 # Pi
G = 9.806 #Gravity acceleration
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 #convert initial angle to radians
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): # equation of motion
dydx = np.zeros_like (state)
dydx [0] = state [1] #dydx [0] = om1
delta1 = state [0] -state [2] #Delta 1delta2 = state [2] -state [4] #Delta 2
delta3 = state [4] -state [0] # delta3
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] #Time generation
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 ()``````

Then

``y = solve_ivp (func, t_span, state, t_eval = t) .y``
the code here
``````sol = solve_ivp (func, t_span, state, t_eval = t)
y = sol.y
print (sol.status)``````

When rewritten like this, -1 was output. When I examine it, it looks like an error called Integration step failed, but I don't know how to solve it.