python 微分方程式の計算時にエラーと警告が出ます

前提・実現したいこと

プログラム中の,vector_solove関数内に記述した微分方程式をsolve_ivpを使って計算させようと考えています．

発生している問題・エラーメッセージ

プログラム内の，ηの値を1より大きくするとエラーも警告もなく動くのですが，
1未満にすると下記に示すエラーと警告が出ます．
（やりたいのはηの値を1以下にいたいと考えています）

PS C:\Users\□□□□\Desktop\○○○○> python mass-spring-FTSM.py
mass-spring-FTSM.py:62: RuntimeWarning: invalid value encountered in double_scalars
  nextu = (m/k*c2)*((c1*(f_x1x2(x1,x2))+0.02*np.sin(t)-(k/m)*u)+c2*(df_x1x2(x2,u)+0.02*np.cos(t))+κ*np.sign(σ(x1,x2,u,t))*abs(σ(x1,x2,u,t))**η+ɤ*σ(x1,x2,u,t)+χ*np.sign(σ(x1,x2,u,t))+(c1*(x2-0.02*np.cos(t))+c2*(f_x1x2(x
1,x2)+0.02*np.cos(t)-(k/m)*u))*(λ+μ(x2)*η*s1(x1,x2,t)**(η-1)))
mass-spring-FTSM.py:54: RuntimeWarning: invalid value encountered in sign
  return -(k/m)*x2-μ(x2)*g*np.sign(u)
mass-spring-FTSM.py:62: RuntimeWarning: invalid value encountered in sign
  nextu = (m/k*c2)*((c1*(f_x1x2(x1,x2))+0.02*np.sin(t)-(k/m)*u)+c2*(df_x1x2(x2,u)+0.02*np.cos(t))+κ*np.sign(σ(x1,x2,u,t))*abs(σ(x1,x2,u,t))**η+ɤ*σ(x1,x2,u,t)+χ*np.sign(σ(x1,x2,u,t))+(c1*(x2-0.02*np.cos(t))+c2*(f_x1x2(x
1,x2)+0.02*np.cos(t)-(k/m)*u))*(λ+μ(x2)*η*s1(x1,x2,t)**(η-1)))
mass-spring-FTSM.py:61: RuntimeWarning: invalid value encountered in sign
  nextx2 = -(k/m)*x1 + d(t) - μ(x2)*g*np.sign(x2) - (k/m)*u
mass-spring-FTSM.py:51: RuntimeWarning: invalid value encountered in sign
  return -(k/m)*x1-μ(x2)*g*np.sign(x2)
mass-spring-FTSM.py:48: RuntimeWarning: invalid value encountered in sign
  return c1*(x2-0.02*np.cos(t))+c2*(-(k/m)*x1+d(t)-μ(x2)*g*np.sign(x2)-(k/m)*u+0.02*np.cos(t))
Traceback (most recent call last):
  File "mass-spring-FTSM.py", line 72, in <module>
    v = solve_ivp(vector_solove, t_span, [x1_f, x2_f, u_f ], method = 'RK45', t_eval = t)
  File "C:\Users\□□□□\AppData\Local\Programs\Python\Python37\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 574, in solve_ivp
    ts = np.hstack(ts)
  File "C:\Users\□□□□\AppData\Local\Programs\Python\Python37\lib\site-packages\numpy\core\shape_base.py", line 340, in hstack
    return _nx.concatenate(arrs, 1)
ValueError: need at least one array to concatenate

該当のソースコード

python
1import matplotlib
2#matplotlib.use('Agg')
3import matplotlib.pyplot as plt
4from scipy.integrate import odeint
5from scipy.integrate import solve_ivp
6import numpy as np
7import pylab
8from openpyxl import Workbook
9
10ms = 0.26
11mc = 0.85
12k = 177
13g = 9.8
14m = mc + (1/3)*ms
15e_integ = 0
16c1 = 1
17c2 = 0.05
18λ = 10
19κ = 2.5
20ɤ = 5.5
21χ = 5
22β = 2
23η = 0.7
24
25u_rec = []
26t_rec = []
27
28def μ(x2):
29    if x2==0:
30        return 0.24
31    else:
32        return 0.11
33
34def d(t):
35    if t < 20:
36        return 0
37    elif t >= 20:
38        return 5*np.sin(3.14*t)
39    
40def σ(x1,x2,u,t):
41    return s2(x1,x2,t,u)+λ*s1(x1,x2,t)+μ(x2)*s1(x1,x2,t)
42    #return 1
43
44def s1(x1,x2,t):
45    return c1*(x1-0.02*np.sin(t))+c2*(x2-0.02*np.cos(t)) 
46
47def s2(x1, x2, t, u):
48    return c1*(x2-0.02*np.cos(t))+c2*(-(k/m)*x1+d(t)-μ(x2)*g*np.sign(x2)-(k/m)*u+0.02*np.cos(t))
49
50def f_x1x2(x1,x2):
51    return -(k/m)*x1-μ(x2)*g*np.sign(x2)
52
53def df_x1x2(x2, u):
54    return -(k/m)*x2-μ(x2)*g*np.sign(u)
55
56
57def vector_solove(t, state):
58
59    x1, x2, u = state
60    nextx1 = x2
61    nextx2 = -(k/m)*x1 + d(t) - μ(x2)*g*np.sign(x2) - (k/m)*u
62    nextu = (m/k*c2)*((c1*(f_x1x2(x1,x2))+0.02*np.sin(t)-(k/m)*u)+c2*(df_x1x2(x2,u)+0.02*np.cos(t))+κ*np.sign(σ(x1,x2,u,t))*abs(σ(x1,x2,u,t))**η+ɤ*σ(x1,x2,u,t)+χ*np.sign(σ(x1,x2,u,t))+(c1*(x2-0.02*np.cos(t))+c2*(f_x1x2(x1,x2)+0.02*np.cos(t)-(k/m)*u))*(λ+μ(x2)*η*s1(x1,x2,t)**(η-1)))
63    
64    return nextx1, nextx2, nextu
65    
66x1_f = 0.0
67x2_f = 0.0
68x3_f = 0.0
69u_f = 0.0
70t_span = [0.0, 40.0]
71t = np.arange(0.0, 40.0, 0.01)
72v = solve_ivp(vector_solove, t_span, [x1_f, x2_f, u_f ], method = 'RK45', t_eval = t)
73
74
75t_vec = v.t
76x1_vec = v.y[0,:]
77x2_vec = v.y[1,:]
78
79
80t = np.arange(0.0, 40.0, 0.01)
81traj_p = []
82traj_v = []
83for i in t:
84    traj_v.append(0.02*np.cos(i))
85    traj_p.append(0.02*np.sin(i))
86
87fig, ax = plt.subplots(2, 1, figsize=(10,6))
88
89
90ax[0].set_xlabel('Time[s]')
91ax[0].set_ylabel('Position[m]')
92ax[0].grid()
93ax[0].set_ylim(-0.03, 0.08)
94ax[0].set_xlim(0.0, 40.0)
95ax[0].plot(t_vec, x1_vec, ls = '-', label = 'PID Controler')
96ax[0].plot(t, traj_p, ls = '--', label = 'Desired Path')
97ax[0].legend()
98
99ax[1].set_xlabel('Time[s]')
100ax[1].set_ylabel('Velocity[m/s]')
101ax[1].grid()
102ax[1].set_ylim(-0.1, 0.2)
103ax[1].set_xlim(0.0, 40.0)
104ax[1].plot(t_vec, x2_vec, ls = '-', label = 'PID Controler')
105ax[1].plot(t, traj_v, ls = '--', label = 'Desired Velocity')
106ax[1].legend()
107plt.savefig('Simulation results of position and velocity of mass–spring system_FTSM.png')
108plt.show()
109
110

試したこと

vector_solove関数の中の微分方程式の一部をコメントアウトしながら実行してみると，nextuの最後の(η-1)乗がないとちゃんと動きますが，(η-1)があると上記で示したエラーと警告が出ます．

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

moootoko

2019/10/11 02:33 編集

エラーの情報の2段落目(?)にある mass-spring-FTSM.py:54: RuntimeWarning: invalid value encountered in sign で、符号を返すメソッドで不正な入力値だと返されているようですが、 numpy.signの引数にnumpy.nanが入っているということはないでしょうか。実数が入れば必ず正常に値を返すメソッドなので、変な値が入っているかもと思いました。