ライブライscipyのodeの使い方

Question

前提・実現したいこと

質問を見てくださりありがとうございます．
質問内容が少し違いますが，odeに問題があると思っており，題名が思いつかなかったので上記の題名にしています．

微分方程式をodeで解いたものをグラフに出力したもの(プログラム内の u )と計算させてから余分なものを引いて出力させたもの(プログラム内の　u_rec )をグラフに出力させたのですが（図1：uのグラフ，図2：u_recのグラフ）同じグラフが出ると思ったのですが違うグラフが出力されます．

uとu_recの出力をExcelに出力すると，確かに全体的にuとu_recの値が違っていましたが，uとu_recが同じ値になるようにプログラムを組んだのですが，どこが間違っているのかわかりません．

図1

図2

図1よりも図2はガタガタしたグラフが出力されます

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

該当のソースコード

import matplotlib.pyplot as plt
import scipy.integrate as sp
import numpy as np
import random


ms = 0.26
mc = 0.85
k = 177
g = 9.8
m = mc + (1/3)*ms

u_rec = []
t_rec = []

def d(t):

    if t < 20:
        return 0
    elif t >= 20:
        return 5*np.sin(3.14*t)
def μ(x2):
    if x2==0:
        return 0.24
    else:
        return 0.11
# SMC制御系
def f(x1,x2):
    return -(k/m)*x1-μ(x2)*g*np.sign(x2)

def df(x2,x3):
    return -(k/m)*x2-μ(x3)*g*np.sign(x3)

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

def du(x1,x2,x3,t):
    return 7*(x2-0.02*np.cos(t))+0.01*(x3+0.02*np.sin(t))+30*(x1-0.02*np.sin(t))

def du_g(x1,x2,x3,t):
    return 7*(x2-0.02*np.cos(t))+0.01*(x3+0.02*np.sin(t))+30*(x1-0.02*np.sin(t))

def heaviside(t):
    if(t < 20):
        return 0
    elif(t >= 20):
        return 1.5

def system(t, x):
    # 戻り値用のリスト
    y = [0]*4
    
    #ダイナミクスの記述
    dx1 = x[1]                     
    dx2 = -(k/m)*x[0] -μ(x[1])*g*np.sign(x[1]) - (k/m)*(7* (x[0]-0.02*np.sin(t)) + 0.1*(x[1]-0.02*np.cos(t))+ 30*x[2]) + d(t)
    dx3 = x[0] - 0.02*np.sin(t)
    dxu = 7*(x[1]-0.02*np.cos(t)) + 0.1*(x[2] + 0.02*np.sin(t)) + 30*(x[0]-0.02*np.sin(t))
    
    # 計算結果を返す
    y[0] = dx1
    y[1] = dx2
    y[2] = dx3
    y[3] = dxu
    
    return y

def simulation(x0, end, step):
    x1 = []
    x2 = []
    t = []
    u = []
    
    ode =  sp.ode(system)
    ode.set_integrator('dopri5', method='bdf', atol=1.e-2)
    ode.set_initial_value(x0, 0)
    t.append(0)
    x1.append(x0[0])
    x2.append(x0[1])
    u.append(0)
    while ode.successful() and ode.t < end - step:
        ode.integrate(ode.t + step)
        t.append(ode.t)
        x1.append(ode.y[0])
        x2.append(ode.y[1])
        u.append(ode.y[3])#積分して取り出したu
    return x1, x2, u, t

# パラメータ
x0 = [0.0, 0.0, 0.0, 0.0] # 初期値
end = 40        # シミュレーション時間
step = 0.01     # 時間の刻み幅

# シミュレーション
x1, x2, u, t = simulation(x0, end, step)

#u_recのリストを作成
for i in range(len(t)):
     u_rec.append(-(m/k)*(x2[i]+(k/m)*x1[i]-d(i*0.01)+μ(x2[i])*g*np.sign(x2[i])))

# 結果をグラフ表示

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))


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, ls = '-',color = 'b', label = 'PID Controler')
ax[0].plot(t, traj_p, ls = '--',color ='r', 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, ls = '-', color='b', label = 'PID Controler')
ax[1].plot(t, traj_v, ls = '--', color='r', label = 'Desired Velocity')
ax[1].legend()
#plt.savefig('position and Velocity PID+d.png')
plt.show()


fig, ax = plt.subplots(figsize=(10,4))

ax.set_xlabel('Time[s]')
ax.set_ylabel('Base Position[m]')
ax.grid()
ax.set_ylim(-0.08, 0.18)
ax.set_xlim(0.0, 40.0)
ax.plot(t, u_rec, ls = '-', color='b', label = 'PID Controler')
#ax.plot(t, u, ls = '-', color='b', label = 'PID Controler')
plt.show()

試したこと

ライブラリscipyのodeについて色々調べてみましたがわかりあせんでした．

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

Comments

[tanishi_a]

`u` と `u_rec` の中身が同じはずというのであれば、

1. `u` と `u_rec` それぞれファイルに出して、diff を取る
2. 差の出てる行を確認
3. 差の出てる行で式と値が想定通りなのか確認

という感じで狭めて行くのが良さそうに思えますが。

[NR4200]

excelに出力したところ全体的に違う値が出ていました．

[NR4200]

しかし，プログラム上では同じ値が出るようにしているのですが，どこが間違っているのかわかりません

[tanishi_a]

数式のことは全くわからないですが、3.14と書いているところは、math.piとかnp.piを使わないのは何か意味があるのですか。

[NR4200]

特に意味はありませんが，np.piに変えたとしても変わりません

[tanishi_a]

u_rec がギザギザしているのが問題と言っているのですよね。

ギザギザさせているのは
u_rec.append(-(m/k)*(x2[i]+(k/m)*x1[i]-d(i*0.01)+μ(x2[i])*g*np.sign(x2[i])))

の中の
np.sign(x2[i])

の中の
x2[i]

のようですが、
そうしたら何か辿れたりしますか？

[NR4200]

ギザギザしているのが問題というよりも，uの出力グラフとu_recのグラフが違っているのが問題です．（例えば，図1と図2を比べた時に5秒のところでは図1に対して図２がほぼ直角になっているなどが問題です）．

プログラム上のsystem関数の中のdxuと

system関数内のdx2を出力してから

for i in range(len(t)):
u_rec.append(-(m/k)*(x2[i]+(k/m)*x1[i]-d(i*0.01)+μ(x2[i])*g*np.sign(x2[i])))

上記プログラムで余分なものを引いたときと同じになるようにプログラムを組んだつもりです．

説明がわかりにくくてすいません．
ちなみにx2[i]は問題なかったです．

[tanishi_a]

問題ないというのは、どう確認しているのですか？

[NR4200]

プログラム中
x1, x2, u, t ,u_s = simulation(x0, end, step)
の戻り値x2とnp.sign(x2[i])をエクセルに保存し

np.sign(x2[i])がx2の値によって-1，1，0の値を出力しているかの確認

エクセルに保存した数値でグラフを作成したx2のグラフとmatplotlibで出力されるグラフの比較し同じグラフだったため問題ないと判断しました．

やり方が違うのであればtanishi_aさんのやり方を教えていただければ幸いです．
よろしくお願いします．

[tanishi_a]

ありがとうございます。

答えは持ってないのですが、
グラフはいったん置いておいて、その式(appendで足してるやつ)のどこが違っているか、誤差が出てるのかとか式を分解？して検算していくしかないのかなと思ってます。
1行ずつprintして確認するイメージです。

[NR4200]

tanishi_aさんへ

printしたりしながら色々といじっりました．
プログラム内ode.set_integratorの引数atolを1.e-4にしたらギザギザが少し収まったのですが，
atolについて何か知っていることありますか？

ちなみにpythonのホームページでは以下のように示されていました．
rtol, atol　float, optional

The input parameters rtol and atol determine the error control performed by the solver. The solver will control the vector, e, of estimated local errors in y, according to an inequality of the form max-norm of (e / ewt) <= 1, where ewt is a vector of positive error weights computed as ewt = rtol * abs(y) + atol. rtol and atol can be either vectors the same length as y or scalars. Defaults to 1.49012e-8.

訳してみたも意味が分かりませんでした

[tanishi_a]

とある行で、引数の指定によって、
想定の値と出てくる結果の誤差が少なくなる、ということですよね？
絞れてきたように見えるので、そこだけに絞り込んだ内容を質問にして、具体的な値を出して質問文を書き直すか新しい質問として出すかすると、わかる人がいるかもしれないな、と思います。
（いつも一般的なことしか言えてなくてすみません）

[NR4200]

tanishi_aさんへ

貴重なご意見ありがとうございます．tanishi_aさんのご意見にはいつも助けられています．

atolについて新しい質問として挙げてみます．