ode45で得る結果の刻み幅を、より小さくしたい

ode45で微分方程式を解き、

刻み幅をより細かくしてその結果をプロットしようとしているのですがうまくいきません。

（細かくなっていない。）

この質問は、以下のURLの質問の延長線上にあります。

プログラムの内容もほとんど同じです。

コーディングの背景、詳しい内容を知りたい方は下記URLをご覧ください。

https://jp.mathworks.com/matlabcentral/answers/503228-ode

以下のＵＲＬ2つを参考にして、

odesetの"Refine"によって刻み幅を小さくしようとしましたが、どうもプロットされるグラフが細かくなっていないようです。

https://jp.mathworks.com/help/matlab/ref/odeset.html

https://jp.mathworks.com/matlabcentral/answers/462566-ode-45-odeset-nonnegative?s_tid=srchtitle

以下に現在のコードを示します。

clc
clear
close all

%-----慣性モーメントの算出-----
%ただしこの時、推進機やカウンターウェイト類を質点、振り子の腕を角柱,と仮定する。

%振り子の腕についての慣性モーメント
%3辺の長さを決定する
%振り子の回転方向は、bの方向に可動であるものとす
%単位は[m]とする

%減衰係数
%ze=(c/2)*sqrt(1/(I*k))
%固有振動数
%om=(c/2)*sqrt(1/(I*k))

%推力測定器の振動方程式
%仮定
syms  ze om I F(t) L s(t)
assume([I L s(t)]>0)
assume([L]<0.261)
assume([F(t)]>=0)


syms sita(t)
assume(abs(sin(sita(t))-sita(t))<=0.01 & abs(1-(1/2)*sita(t)^2)<=0.01)

%方程式の宣言
%Iθ''+2ζωnθ'+ωn^2θ=FtL/I
ds2=diff(sita,t,2)
ds1=diff(sita,t,1)
eqn = ds2+2*ze*om*ds1+(om^2)*sita(t) == F(t)*L/I

[V] = odeToVectorField(I*ds2+2*ze*om*ds1+om^2*sita(t)==F(t)*L/I)

%初期値の設定
cond1=F(0)==0
cond2=sita(0)==0
cond3=ds1(0)==0
cond4=ds2(0)==0

M = matlabFunction(V,'vars', {'t','Y','I','L','om','ze','F'})
%慣性モーメントの算出
%各値の宣言
g=9.8
m1=2.5%推進機とマウントの質量[kg] 
m2=3    %カウンターウェイトの質量[kg]
m3=0.6  %ふりこの腕の質量[kg]
r1=0.04   %円筒型推進機の半径[m]
r2=0.04   %円筒型カウンターウェイトの周方向半径[m]
L=0.26   %振り子の腕の全長[m] <261[mm]
Lcm1=0.155 %支点からの推進機の距離[m]
Lcm2=L-Lcm1 %支点からのカウンターウェイトの距離[m]
l=0.04    %カウンターウェイトの直線方向の長さ[m]
%推進機部分の慣性モーメント単体
%推進機の形状は円筒型とする。
J1=(1/2)*m1*r1^2+m1*Lcm1^2
%カウンターウェイトの慣性モーメント
J2=((m2*(3*r2^2+(l/2)^2))/12)+m2*Lcm2^2
%振り子の腕の慣性モーメント単体
J3=(m3*L^2)/12
%結果として全体の慣性モーメントは
I=J1+J2+J3 %[kg・m^2]

j=[0.01 0.001 0.0001]%支点のころがり摩擦係数
%F(t)=10^-9 %おおよその理論値を使用
F(t)=10^-9
p=((m1*g*Lcm1)+(F(t)*Lcm1))-((m2*g*Lcm2)+((m1+m2+m3)*g*j))
%↑支点の摩擦力があっても振れることが可能かをトルクのつり合い域から求める
p=double(p)
ze=0.00000000001 %減衰係数

%k=m2*g*Lcm2-m1*g*Lcm1
k=0.049
om=sqrt(k/I)

%F(t)=double(F(t))
%微小角度は10度以内とする( 引用：http://www.ll.em-net.ne.jp/~m-m/reference/smallAngle/smallAngleApprox.htm )
%'t','Y(つまりθ)','I','L','om','ze''F(t)'に対応する値の範囲として
myfun = @(t,y) double(M(t,y,I,L,om,ze,F));

%以下、問題となっている箇所
opts = odeset('Refine',10)
sol=ode45(@(t,y) myfun(t,y),[0 50],[0 10])

% グラフ化
plot(sol.x, sol.y(1,:))
xlabel('t')
ylabel('振れ角θ')
title('時間,s')

Answer 1

[michio]

odeset: https://jp.mathworks.com/help/matlab/ref/odeset.html

の Refine の項目を見ると、以下の記載がありました。

Refine は length(tspan) > 2 の場合や ODE ソルバーが解を構造体として返す場合には適用されません。

ですので、構造体ではなく数値を返すよう ode45 を実行して見てください。

例：

[t,y] = ode45(@(t,y) myfun(t,y),[0 50],[0 10])

[Noruji Muto]

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

以下のように修正してみました。（関連個所のみ抜粋）

myfun = @(t,y) double(M(t,y,I,L,om,ze,F));
opts = odeset('Refine',10)

 tspan=[0 10]
 y0=[0 0] %θの初期値
sol=ode45(@(t,y) myfun(t,y),tspan,y0)
%  tspan=[0 10]
%  y0=[0 0] %θの初期値
 [t,y]=ode45(@(t,y) myfun(t,y),tspan,y0)

% グラフ化
%plot(sol.x, sol.y(1,:))
plot(t,y(:,1))
xlabel('t')
ylabel('振れ角θ')
title('時間,s')

ｙについて、元々[0 10]としていましたが、[ ] の中身は初期値であるというのをどこかで見かけたので、初期値として[0 0]に修正しました。

（この認識正しいでしょうか？）

odesetを使用する前後のグラフを比較してみました。

↓（前）

↓後

どうもデータ点数・滑らかさが変わっていないようです。

何かコーディングがおかしいのでしょうか？

[michio]

パッと気が付きませんでしたが、ode45 を実行するときに オプションも引数に加えてあげる必要があります。

opts = odeset('Refine',10)
tspan=[0 10]
y0=[0 0] %θの初期値
[t,y]=ode45(@(t,y) myfun(t,y),tspan,y0,opts)

これで結果が変わると見ているんですが、データ点数を確認するときは、t や y の配列サイズを確認したほうが確実ですので、そこに注目してみてください。

[Noruji Muto]

その方法で解決しました。

かなり滑らかになりましたね。

ありがとうございました！