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

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

Answer 1

[Musashi Ito]

ソルバーについてあまり詳しくないですが、ode15s のソルバーの選択が必須でしょうか？

実験的に他のソルバーも試して、得られた解の値が想定に近いものを選択する方法はいかがでしょうか。

[Kaz]

回答ありがとうございました。

ode15sの他に使用できそうなode23tbソルバーを試してみます。