odeで取得した解を​プロットするためのs​truct→doub​le変換ができない

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

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

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

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

二階常微分方程式をodeによって解き、

その結果をグラフにして可視化しようとしているのですがうまくいきません。

上記URLでいただいた回答を活用した上で、plot()を用いてグラフ化を試みました。

プロットは横軸がｔ（時間）、縦軸がθを表すようにプロットしようとしています。

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

clc
clear
close all

%振動方程式
%仮定
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)
%ze＝ζ
%om=ωn
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'})
%慣性モーメントの算出
%各値の宣言
m1=1    %推進機の質量[kg] 
m2=3    %カウンターウェイトの質量[kg]
m3=0.6  %ふりこの腕の質量[kg]
r1=0.04   %円筒型推進機の半径[m]
r2=0.04   %円筒型カウンターウェイトの周方向半径[m]
L=0.26   %振り子の腕の全長[m] <261[mm]
Lcm1=0.24 %支点からの推進機の距離[m]
Lcm2=0.200 %支点からのカウンターウェイトの距離[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]

ze=0 %減衰係数
g=9.8
k=m2*g*Lcm2-m1*g*Lcm1
om=sqrt(k/I)
F(t)=10^-9 %おおよその理論値を使用
%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));
sol=ode45(@(t,y) myfun(t,y),[0 50],[0 10])

% グラフ化
plot(t,sol)    %エラー発生箇所
xlabel('t')
ylabel('振れ角θ')
title('時間,s')

これに対し以下のようにエラーが出ました。

エラー: plot
データは数値、datetime、duration であるか、double に変換可能な配列でなければなりません。

エラー: thrust_stand_ode (line 78)
plot(t,sol)
 

このエラーに対し、得た解(上記コードでいう所のsol)をdouble型に変換しようとしました。

myfun = @(t,y) double(M(t,y,I,L,om,ze,F));
sol=ode45(@(t,y) myfun(t,y),[0 50],[0 10])%0~50秒の間の現象で、θは10°以内かつできるだけ大きくしたい。
sol=double(sol)　%ここを追加しました。

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

このコードを走らせて、下記のようなエラーが発生しました。

エラー: double
struct から double に変換できません。

エラー: thrust_stand_ode (line 75)
sol=double(sol)

方程式内の文字で数値を定義できていないものがあるのかと思いましたが、

Iθ''+2ζωnθ'+ωn^2θ=FtL/I 　という解きたい方程式の中で、

一通り数値は定義できているのでdouble化は問題なくできるはずと考えております。

（Ftは、F(t)>0の時のF(t)の値です。）

現状、struct型をdoubleに変換する方法について詳しい解説は見つかっていないこともあり進展がありません。

Answer 1

[michio]

sol は構造体でこんなデータになっているはずです。(ステップ数によって具体的な配列サイズは違います。)

>> sol
sol = 
  フィールドをもつ struct:

     solver: 'ode45'
    extdata: [1×1 struct]
          x: [1×19 double]
          y: [2×19 double]
      stats: [1×1 struct]
      idata: [1×1 struct]      

ODEの解は y に入っているので振れ角が1つ目であれば

plot(sol.x, sol.y(1,:))    

両方プロットするなら

plot(sol.x, sol.y)

でOKです。もし時刻 t と解 y を double 型で求めるなら、

[t,y]=ode45(@(t,y) myfun(t,y),[0 50],[0 10])%0~50秒の間の現象で、θは10°以内かつできるだけ大きくしたい。

としてみてください。

plot(t,y)

でプロットできると思います。

[優佑 熊谷]

ｙが二つ出力されると思うのですが、一つ目が振れ角が一つ目、二つ目は何になるのでしょうか。

Answer 2

[Noruji Muto]

お陰様でプロットに成功しました。

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