質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

342閲覧

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

Tanaka.sot

総合スコア18

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

1グッド

0クリップ

投稿2019/11/17 05:31

編集2019/11/18 13:29

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

python

1from numpy import sin,cos 2from scipy.integrate import solve_ivp 3import numpy as np 4import csv 5from fractions import Fraction as fr 6 7Pi =3.14 #円周率 8G =9.806 #重力加速度 9f=open(result.csv,'a',newline="") 10writer=csv.writer(f) 11th1_DEG =90 12th2_DEG =90 13th3_DEG =90 14om1_DEG =0 15om2_DEG =0 16om3_DEG =0 17L1 =1 18L2 =1 19L3 =1 20M1 =1 21M2 =1 22M3 =1 23th1 =th1_DEG*Pi/180 #初期の角度をラジアンに変換 24th2 =th2_DEG*Pi/180 25th3 =th3_DEG*Pi/180 26om1 =om1_DEG*Pi/180 27om2 =om2_DEG*Pi/180 28om3 =om3_DEG*Pi/180 29 30def func(t,state): #運動方程式 31 dydx=np.zeros_like(state) 32 dydx[0]=state[1] #dydx[0]=om1 33 delta1=state[0]-state[2] #デルタ1 34 delta2=state[2]-state[4] #デルタ2 35 delta3=state[4]-state[0] #デルタ3 36 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)* 37 sin(delta1)+2*L2*M3*sin(delta2)*cos(delta3))+state[1]*state[1]*2*L1*M3*sin(delta3)*cos(delta3)+(2*M3*cos(delta3)* 38 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)* 39 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)+ 40 2*M3*cos(delta1)-2*M3*cos(delta2)*cos(delta3)))) 41 dydx[2]=state[3] #dydx[2]=om2 42 dydx[3]=(state[1]*state[1]*L1*(2*M3*M3*sin(delta3)*cos(delta2)*cos(delta3)*cos(delta3)+M3*(M2+2*M3)*sin(delta1)* 43 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)* 44 (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)* 45 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)* 46 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)* 47 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)- 48 fr(1,4)*M3*(M2+4*M3)*cos(delta2)*cos(delta3))) 49 dydx[4]=state[5] #dydx[4]=om3 50 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* 51 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))- 52 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)* 53 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)) 54 *(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)* 55 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)* 56 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])* 57 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))* 58 (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* 59 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)* 60 cos(delta2))*(M3*cos(delta2)*cos(delta3)-fr(1,2)*(M2+2*M3)*cos(delta1)))) 61 62 return dydx 63t_span =[0.0,10] #時間生成 64dt =0.01 65t =np.arange(t_span[0],t_span[1],dt) 66writer.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)]) 67 68state =[th1,om1,th2,om2,th3,om3] 69print(len(t)) 70y =solve_ivp(func,t_span,state,t_eval=t).y 71for tt,th1,th2,th3 in zip(t,y[0,:],y[2,:],y[4,:]): 72 x1=L1 *sin(th1) 73 y1=-L1*cos(th1) 74 x2=x1+L2*sin(th2) 75 y2=y1-L2*cos(th2) 76 x3=x2+L3*sin(th3) 77 y3=y2-L3*sin(th3) 78 X1= fr(1,2)*L1*sin(th1) 79 Y1=-fr(1,2)*L1*cos(th1) 80 X2=x1+fr(1,2)*L2*sin(th2) 81 Y2=y1-fr(1,2)*L2*cos(th2) 82 X3=x2+fr(1,2)*L3*sin(th3) 83 Y3=y2-fr(1,2)*L3*cos(th3) 84 Row_list=[float(tt),float(x1),float(y1),float(X1),float(Y1),float(x2) 85 ,float(y2),float(X2),float(Y2),float(x3),float(y3),float(X3),float(Y3)] 86 writer.writerow(Row_list) 87f.close()

そのあと、

python

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

ここのコードを

python

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

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

退会済みユーザー👍を押しています

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

自己解決

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

投稿2019/11/25 09:28

編集2019/12/22 16:17
Tanaka.sot

総合スコア18

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問