symfunをfloat(single, double)に変換する方法

10 views (last 30 days)
Noruji Muto
Noruji Muto on 19 Feb 2020
Commented: michio on 2 Mar 2020
お世話になります。
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とすることもできない事を前提としているので、どうしたらいいのか分からずにいます。

Accepted Answer

michio
michio on 21 Feb 2020
ひとまず、いまコメントアウトされている行
%myfun = @(t,y) double(M(t,y,I,L,om,ze,F));
を使用すればエラーは回避できるようですが、いかがでしょう?
  2 Comments
Noruji Muto
Noruji Muto on 2 Mar 2020
返信が遅くなり申し訳ありません。
それで回避できました。
ありがとうございました。
michio
michio on 2 Mar 2020
よかったです。
コメント頂きありがとうございました。

Sign in to comment.

More Answers (0)

Categories

Find more on プログラミング in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!