pythonライブラリscipyのodeintを使ってグラフを表示するとグラフが正しく表示されない

Question

前提・実現したいこと

図1のようにな摩擦のあるマスバネ系の運動方程式を微分方程式でしめして，図1に示す台車の位置と速度ををグラフにしようとしています．
運動方程式は
mx''= -Kx - μmgsgn(x')-ku
(x''は時間における2階微分，x'は1階微分を表します，sgn(.)は符号関数)

SciPy の odeモジュールは１階の微分方程式しか解けないので，上式を
x = x1
x1' = x2の関係を利用して，

x1' = x2
x2' = -(k/m)x1 - μgsgn(x2) - (k/m)u

に式変形してからプログラム上のvector関数内に記述しています．
（プログラム上ではu=((7*(x1-0.02np.sin(t)))+(0.1(x2-0.02np.cos(t)))+30e_integ) ）- (m/k)d(t)が代入されています．）
図1

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

摩擦をいれなければ（プログラム中のvector関数の中のコメントアウトしているnext2を使ったとき）図2のようにグラフが出るのですが
摩擦（プログラム中のvector関数の中のコメントアウトしていないnext2を使ったとき）を入れると図3のように常に0のグラフが出ます．

また警告が出ているので，図4にコマンドプロンプトを示しておきます．

該当のソースコード

import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
import pylab


ms = 0.26
mc = 0.85
k = 177
g = 9.8
μs = 0.24
μk = 0.11

m = mc + (1/3)*ms
e_integ = 0

def μ(x2):
    if x2==0:
        return 0.24
    else:
        return 0.11

def d(t):

    if t < 20:
        return 0
    elif t >= 20:
        return 5*np.sin(3.14*t)

def vector(state, t, g, k, m, e_integ):

    x1, x2 = state

    e_integ += x1 - 0.02*np.sin(t) 
    nextx1 = x2
    nextx2 = -(k/m)*x1 - (k/m)*((7*(x1-0.02*np.sin(t)))+(0.1*(x2-0.02*np.cos(t)))+30*e_integ) + d(t) - μ(x2)*g*np.sign(x2)
    #nextx2 = -(k/m)*x1 - (k/m)*((7*(x1-0.02*np.sin(t)))+(0.1*(x2-0.02*np.cos(t)))+30*e_integ) + d(t) 
   

    return nextx1, nextx2
    

x1_f = 0
x2_f = 0

t = np.arange(0.0, 40.0, 0.01)

v = odeint(vector, [x1_f, x2_f ], t, args=(g, k, m, e_integ))
print(v)
x1_vec = v[:,0] 
x2_vec = v[:,1]
'''
traj_p = []
traj_v = []
for i in t:
    traj_v.append(0.02*np.cos(i))
    traj_p.append(0.02*np.sin(i))
'''
fig, ax = plt.subplots(2, 1, figsize=(10,6))
#fig.tight_layout()

ax[0].set_xlabel('Time[s]')
ax[0].set_ylabel('Position[m]')
ax[0].grid()
ax[0].set_ylim(-0.03, 0.08)
ax[0].set_xlim(0.0, 40.0)
ax[0].plot(t, x1_vec, ls = '-', label = 'PID Controler')
#ax[0].plot(t, traj_p, ls = '--', label = 'Desired Path')
ax[0].legend()

ax[1].set_xlabel('Time[s]')
ax[1].set_ylabel('Velocity[m/s]')
ax[1].grid()
ax[1].set_ylim(-0.1, 0.2)
ax[1].set_xlim(0.0, 40.0)
ax[1].plot(t, x2_vec, ls = '-', label = 'PID Controler')
#ax[1].plot(t, traj_v, ls = '--', label = 'Desired Velocity')
ax[1].legend()
plt.show()

試したこと

pythonライブライブscipyのodeintが間違っているのかと思い調べてみましたが，間違った使い方はしていませんでした．
また正しい結果が図3のようになるのかと思い色々調べてみましたが，ほとんど図２と同じ結果になるそうです

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

図２

図3

図4

Comments

[t_obara]

これ、実行するとWARNINGが出ませんか？もし出ているのであれば、質問を編集してその旨記載した方が回答が得られやすいかと思います。もちろん、そのメッセージでググってみて、ご自身の問題と関連があるのかを確かめられるのもお勧めいたします。

[NR4200]

アドバイスありがとうございます．
確かに警告が出ていたので図4として載せます．

[nandymak]

＞そのメッセージでググってみて、ご自身の問題と関連があるのかを確かめられるのもお勧めいたします。
関連はなかったのでしょうか？

Answer 1

摩擦項の符号をプラスにすれば解決するのではないでしょうか。以下のコードにしたところ、それっぽいグラフが出ました(物理に基づかない回答ですみません)。

nextx2 = -(k/m)*x1 - (k/m)*((7*(x1-0.02*np.sin(t)))+(0.1*(x2-0.02*np.cos(t)))+30*e_integ) + d(t) + μ(x2)*g*np.sign(x2)

Comments

[NR4200]

回答してくださってありがとうございます．

恐らく，摩擦項の符号を逆にしてしまうと物理的に成り立たないので＋にすることはできないと思います．

[NR4200]

早急に回答してくださったのに申し訳ありません