symfunをfloat（single, double）に変換する方法

お世話になります。

odeを用いて二階微分振動方程式を解き、解析を行おうとしているのですが、

エラーが発生して進めずにいます。

この質問は以下の質問に関連していますが、全体の流れは同じであるものの、

数値など細かい点は変更があります。（計算の流れそのものはあまり変わりありません）

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

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

clc
clear
close all
%-----慣性モーメントの算出-----
%推力測定器の振動方程式
%仮定
syms  ze om I F(t) L s(t)
assume([I L s(t)]>0);
assume(L<0.45);
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'});
%各値の宣言
F(t)=10^-9; %おおよその理論値を使用         %おそらくここが悪さをしている
%F(t)=doubel(F(t))
g=9.8;%重力加速度
m1=2.2 ;  %推進機の質量[kg] 
L=0.45  ; %振り子の腕の全長[m] <0.45[m]
m2=0.5*L;%振り子の腕の質量（1mあたり0.5kg ( https://jp.misumi-ec.com/vona2/detail/110302683830/?KWSearch=%e3%82%a2%e3%83%ab%e3%83%9f%e3%83%95%e3%83%ac%e3%83%bc%e3%83%a0&searchFlow=results2products ) ）
r1=0.04;   %円筒型推進機の半径[m]
Lcm1=0.26; %支点から推進機の距離[m]
Lcm2=L-Lcm1 ;%支点からばねまでの距離[m]
a=0.03;%支点から推進機がついていない側への腕の長さ（仮定）
%慣性モーメントの算出
%推進機部分の慣性モーメント単体
%推進機の形状は円筒型とする。
J1=(1/2)*m1*r1^2+m1*Lcm1^2;
%振り子の腕の慣性モーメント単体
J2=(m2*L^2)/12;
%結果として全体の慣性モーメントは
I=J1+J2 ;%[kg・m^2]

j=0.001;%転がり摩擦係数(仮定
T=(m1+m2)*g*j;%転がり抵抗としてのトルク
sita_ss=deg2rad(3)%定常振れ角として欲しい角度
baneT=(m1*g*Lcm1*sin(sita_ss)+F(t)*Lcm1)/(Lcm2^2*(sita_ss))%トルクの平衡を保つために必要な
Lm1=0
Lm1=sym(Lm1)
Lm1=((m1*g)*(Lcm1+a)+(m2*g)*(L-a))/((m1+m2)*g);
k1=((m1)*g*sin(sita_ss)+F*(Lcm1))/(Lcm2^2*sita_ss);%ばねに必要なばね定数
k2=F(t)*L/sita_ss;%定常振れ角として設定した角度から測定器の復元トルクを得る。
%k1=bane*Lcm2^2*sita(t);%ばねによるトルク
p=((m1*g*Lcm1)+(F(t)*Lcm1))-(k1+baneT);%推進機稼働時のトルク平衡式
%p=double(p)
%ze=0.02 %減衰係数
ze=0;
ze=double(ze)

k2=k1-m1*g*Lcm1;  %測定器そのものの復元トルク
om=sqrt(k2/I);
om=double(om)
de1=F(t)/k2;%静的たわみ
%de1=double(de1)

%腕のたわみ設計（両端支持梁として）
%腕にはミスミのアルミフレーム20*20cmを使用( https://jp.misumi-ec.com/vona2/detail/110302683830/?KWSearch=%e3%82%a2%e3%83%ab%e3%83%9f%e3%83%95%e3%83%ac%e3%83%bc%e3%83%a0&searchFlow=results2products )
P=F(t)+(m1+m2)*g;
Id=0.0000000074 ;%断面二次モーメント[m^4]
E=69972000000 ;  %ヤング率[N/m2](ミスミ曰く定数らしい)
de2=(P*(Lcm1^2)*(Lcm2^2))/(3*E*Id*L);

%sita_ss=double(sita_ss)
%myfun = @(t,y) double(M(t,y,I,L,om,ze,F));
myfun = @(t,y) (M(t,y,I,L,om,ze,F));
opts = odeset('Refine',10);
%t=double(t)
sol=ode45(@(t,y) myfun(t,y),[0 100],[0 0]);%%ここでエラー発生
%[t,y]=ode45(@(t,y) myfun(t,y),[0 10],[0 0],opts)
%plot(sol.x, sol.y(1,:))   
plot(sol.x, sol.y(1,:))  ;
xlabel('t');
ylabel('振れ角θ');
title('時間,s');

このコードを走らせたとき、

以下のようなエラーが発生します

エラー: odearguments (line 113)
入力は float、つまり single か double でなければなりません。

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

エラー: standform2_ode (line 81)
sol=ode45(@(t,y) myfun(t,y),[0 100],[0 0]);

おそらく指定した値であるF(t)が悪さをしているのだろうと予想できるのですが、

t≦0の時にF(0)=0, t>0の時にF(t)=10^-9となるように条件指定が必要なので,率直にF=10^-9とすることもできない事を前提としているので、どうしたらいいのか分からずにいます。

Answer 1

[michio]

ひとまず、いまコメントアウトされている行

%myfun = @(t,y) double(M(t,y,I,L,om,ze,F));

を使用すればエラーは回避できるようですが、いかがでしょう？

[Noruji Muto]

返信が遅くなり申し訳ありません。

それで回避できました。

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

[michio]

よかったです。

コメント頂きありがとうございました。