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

18 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'が未定義です」という表示も出ていたのですが、
色々弄っているうちに無くなっていました。

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
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.

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!