MATLAB Answers

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

62 views (last 30 days)
Noruji Muto
Noruji Muto on 2 Feb 2020
Commented: Noruji Muto on 5 Feb 2020
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'が未定義です」という表示も出ていたのですが、
色々弄っているうちに無くなっていました。

  0 Comments

Sign in to comment.

Answers (1)

michio
michio on 3 Feb 2020
すいません、詳細は追えていないんですが、
最後の 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 で変換するようにしています。

  3 Comments

Noruji Muto
Noruji Muto on 5 Feb 2020
ご回答ありがとうございます。
おかげさまで次のステップに進むことができました。
余談
odeでは@(t,y)といった感じの「@」がつく文をちらほら見かけますが、
調べてみても各種@がつく文について解説しているwebページが見つかりません。
もし何かodeに使用される@がつく文の勉強に役立つサイト等ご存じでしたら
ご教示いただけると幸いです。
michio
michio on 5 Feb 2020
よかったです。
@ がつくのは無名関数と呼ばれるものです。
ode で勾配や fminsearch で目的関数を入力引数として与える際によく登場します。
@ の後の @(t,y) を省略できたりするのが混乱のもとだと感じますが、便利ですので是非上のページ見てみてください。
省略例ですが、と ode45 の第一引数に勾配を計算する関数を与えますが、ode45 では この関数の入力として t, y を想定しています。ですので、fun が もともと t,y を入力として持つ関数であれば
sol=ode45(@fun...
と与えることができますが、これは
sol=ode45(@(t,y) fun(t,y)...
と同じ意味になります。
Noruji Muto
Noruji Muto on 5 Feb 2020
ご丁寧にありがとうございます。
参照してみます。

Sign in to comment.

Sign in to answer this question.