MATLAB Answers

ばね駆動の二重振り子​をode15sソルバ​ーで数値計算させると​失敗します

9 views (last 30 days)
Kaz
Kaz on 1 Feb 2021
Commented: Kaz on 25 Feb 2021
MATLAB R2020aにて、以下の条件および図の二重振り子の運動を数値計算しようとしています。
  • 系はx軸方向に並進運動する
  • 駆動力はばね(二重振り子のそれぞれの角度, に比例する量)
  • それぞれの質点に対して、速度に比例する抵抗係数を乗じている
  • ラグランジアンを用いて運動方程式を導いた
  • MATLABソルバーは、ode15s
困っている点は、次の通りです。
  • 相対許容誤差、絶対許容誤差を変えても、計算が失敗してしまう
  • 計算が失敗しない場合でも「相対許容誤差」について書かれたページの「ヒント」の項目にあるように、許容誤差の値を変えて数値計算させると、得られた解の値が大きく異なる
許容誤差あるいはソルバーをどのように選べば良いでしょうか。
運動方程式を記述したスクリプトは以下です。
function dx = odefun(t, x, kk1, kk2) % 支点が動く粘性抵抗付き二重バネ振子,無次元化, 抵抗とばねの比
%% 値
kk1 ; kk2 ; % ばね定数は可変
l1 = 10*10^(-6) ; % 腕の長さ
l2 = 10*10^(-6) ;
m0 = 10^(-10) ; % 質点の質量
m1 = l1*10^(-6) ;
m2 = l2*10^(-6) ;
k1 = 10^(-kk1)/l1 ; % ばね定数
k2 = 10^(-(kk1+kk2-2))/l2 ;
Crot = 10^(0) ; % 抵抗係数1
Ctrans = 10^(-7) ; % 抵抗係数2
thetaf1 = 20*pi/180 ; % 終端角度
thetaf2 = 40*pi/180 ;
L = 10^(-4) ; % 無次元化パラメータ
L1 = l1/L ;
L2 = l2/L ;
M = m0+m1+m2 ;
M12 = (m1+m2)/M ; % 無次元化した質量
M2 = m2/M ;
C1 = (Crot*L)/(k1*M)^0.5 ; % 無次元化した抵抗係数1
C2 = ( Ctrans+(l1+l2)*Crot )*(L/(k1*M)^0.5) ;
%% 行列と変数の定義
dx = zeros(size(x)) ;
dx(1) = x(4) ; % x(1)=position
dx(2) = x(5) ; % x(2)=theta_1 根元の角度
dx(3) = x(6) ; % x(3)=theta_2 先端の角度
%% 質量行列の定義
m11 = 1 ;
m12 = -M12*L1*sin(x(2)) ;
m13 = -M2*L2*sin(x(3)) ;
m22 = M12*(L1)^2 ;
m23 = M2*L1*L2*cos(x(3) - x(2)) ;
m33 = M2*(L2)^2 ;
%% 力(水平面より下方に触れようとしたらばね定数の値を大きくする)
kcc = 1.5 ;
if x(2) < 0
kc1 = kcc ;
else
kc1 = 1 ;
end
if x(3) - x(2) < 0
kc2 = kcc ;
else
kc2 = 1 ;
end
%% 無次元化したばねの力
K1 = ( thetaf1 - x(2) ) ;
K2 = ( kc2*k2/(kc1*k1) )*( thetaf2 - x(3) ) ;
%% 無次元化した抵抗力
Rx = -( (C2+C1*(L1+L2))*x(4) - C1*( (L1+L2)*L1*x(5)*sin(x(2)) + (L2^2)*x(6)*sin(x(3)) ) ) ;
R1 = -C1*( -(L1+L2)*L1*x(4)*sin(x(2)) + (L1+L2)*(L1^2)*x(5) + L1*(L2^2)*x(6)*cos(x(3) - x(2)) ) ;
R2 = -C1*( -(L2^2)*x(4)*sin(x(3)) + (L2^3)*x(6) + L1*(L2^2)*x(5)*cos(x(3) - x(2)) ) ;
%% 運動方程式を構成する行列
M = [ m22*m33-m23^2 -(m12*m33-m13*m23) m12*m23-m13*m22 ;
-(m12*m33-m23*m13) m11*m33-m13^2 -(m11*m23-m13*m12) ;
m12*m23-m22*m13 -(m11*m23-m12*m13) m11*m22-m12^2 ] ./ ( m11*m22*m33 + 2*m12*m13*m23 - m22*m13^2 - m11*m23^2 - m33*m12^2 ) ; % 逆行列
A = [ M12*L1*(x(5)^2)*cos(x(2)) + M2*L2*(x(6)^2)*cos(x(3)) + Rx ;
M2*L1*L2*(x(6)^2)*sin(x(3) - x(2)) + K1 + R1 ;
-M2*L1*L2*(x(5)^2)*sin(x(3) - x(2)) + K2 + R2 ] ;
%% 運動方程式
dx(4:6) = M*A ;
end
このスクリプトを実行するために、次のスクリプトを用意しました
clear all;
initCond = [ 0 0 0 0 0 0 ];
opts = odeset( 'RelTol', 1e-6, 'AbsTol', 1e-4 ) ; % 絶対誤差相対誤差の設定
for kk1 = 10:1:20
for kk2 = 1:1:3
solsA = ode15s( @odefun, [0 10^7], initCond, opts, kk1, kk2 ) ;
A = solsA.y(1,end) ;
end
end

Accepted Answer

Musashi Ito
Musashi Ito on 15 Feb 2021
ソルバーについてあまり詳しくないですが、ode15s のソルバーの選択が必須でしょうか?
実験的に他のソルバーも試して、得られた解の値が想定に近いものを選択する方法はいかがでしょうか。
  1 Comment
Kaz
Kaz on 25 Feb 2021
回答ありがとうございました。
ode15sの他に使用できそうなode23tbソルバーを試してみます。

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!