この質問は、以下のURLの質問の延長線上にあります。
プログラムの内容もほとんど同じです。
コーディングの背景、詳しい内容を知りたい方は下記URLをご覧ください。
二階常微分方程式を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)
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'})
m1=1
m2=3
m3=0.6
r1=0.04
r2=0.04
L=0.26
Lcm1=0.24
Lcm2=0.200
l=0.04
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
ze=0
g=9.8
k=m2*g*Lcm2-m1*g*Lcm1
om=sqrt(k/I)
F(t)=10^-9
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])
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に変換する方法について詳しい解説は見つかっていないこともあり進展がありません。
