二階常微分方程式をo​deで解こうとする際​に発生する入力引数の​不足

ode45を用いて二階常微分方程式を解こうとしていますが、

なぜか入力引数が不足しているという表示が出てきてしまいます。

入力した文字に対する数値の個数はあっているはずですし、色々弄ってみてはいるものの同じエラーが出てきて困っています。

下の図のようなシーソー型の振り子について、

二階常微分方程式で表される、振り子の振動を表す方程式を解こうとしています。

図左側の機械からは推力が発生し、振り子には重力方向に力が加わって変位が発生し、

図右側のカウンターウェイト（振り子の変位を復元するための重り）によって、

発生した振り子の振れ角θ（図にはありませんが…）は復元され、元の位置に戻ります。

この時推力F(t)は,時間t=0の時は0で、t>0の場合はFtなる推力を定常的に発生させます。

振動方程式は以下のように表現されます。

当初ラプラス変換で解こうとしましたが、頓挫したのでodeの使用を試みています。

以下に現状のコードを記載します。

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)  %θ'の宣言
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    %k
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)'に対応する値の範囲として
sol=ode45(M,[0 50],[0 10],[I],[L],[om],[0],[10^-9]);%ここでエラー発生

エラーが表示されている時、以下のような表示がでます。

入力引数が不足しています。

エラー:
symengine>@(t,Y,I,L,om,ze,F)[Y(2);-1.0./I.^2.*(-L.*F(t)+I.*om.^2.*Y(1)+I.*om.*ze.*Y(2).*2.0)]

エラー: odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

エラー: ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

エラー: thrust_stand_ode (line 73)
sol=ode45(M,[0 50],[0 10],[I],[L],[om],[0],[10^-9]);

元々、「関数または'F'が未定義です」という表示も出ていたのですが、

色々弄っているうちに無くなっていました。

Answer 1

[michio]

すいません、詳細は追えていないんですが、

最後の F(t) = 10^-9 で F(t) を定義した後に、

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]);

とすればエラーは避けられました。

エラーの原因としては symbolic object からの出力を ode45 に入る前に double 型に変換するところ、そして ode45 の第一引数を t,y を入力とする関数として定義する必要がある（他にもかつて方法はあったと思いますが・・とりえあず現時点での ode45 のドキュメンテーションではそのような使い方が紹介されています）点かと思います。

M の入力として F は symbolic である必要がありそうだったので、M の出力を double で変換するようにしています。

[Noruji Muto]

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

おかげさまで次のステップに進むことができました。

余談

odeでは@(t,y)といった感じの「@」がつく文をちらほら見かけますが、

調べてみても各種@がつく文について解説しているwebページが見つかりません。

もし何かodeに使用される@がつく文の勉強に役立つサイト等ご存じでしたら

ご教示いただけると幸いです。

[michio]

よかったです。

@ がつくのは無名関数と呼ばれるものです。

https://jp.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html

ode で勾配や fminsearch で目的関数を入力引数として与える際によく登場します。

@ の後の @(t,y) を省略できたりするのが混乱のもとだと感じますが、便利ですので是非上のページ見てみてください。

省略例ですが、と ode45 の第一引数に勾配を計算する関数を与えますが、ode45 では この関数の入力として t, y を想定しています。ですので、fun が もともと t,y を入力として持つ関数であれば

sol=ode45(@fun...    

と与えることができますが、これは

sol=ode45(@(t,y) fun(t,y)...

と同じ意味になります。

[Noruji Muto]

ご丁寧にありがとうございます。

参照してみます。