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

Question

Kaz 2021-2-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/732248-ode15s

评论： Kaz 2021-2-25

采纳的回答： Musashi Ito

在 MATLAB Online 中打开

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