odedc
離散/連続 ODE ソルバ
呼び出し手順
yt = odedc(y0, nd, stdel, t0, t, f)
引数
- y0
実数列ベクトル (初期条件),
y0=[y0c;y0d]
ただしy0d
はnd
個の要素を有します.- nd
整数,
y0d
の次元- stdel
1個または2個の要素を有する実数ベクトル,
stdel=[h, delta]
(デフォルト値はdelta=0
).- t0
実数のスカラー (初期時間).
- t
実数 (行) ベクトル,
yt
が計算される時間- f
外部,すなわち 関数または文字列または 呼び出し手順を指定するリスト:
yp=f(t,yc,yd,flag)
.- f
外部ルーチン すなわち,関数または文字列またはリスト, 呼び出し手順は以下:
yp=f(t,yc,yd,flag)
- リスト
この形式の外部ルーチンは関数にパラメーラを 指定する際に使用されます. 以下のようにします:
list(f, p1, p2,...)
ただし,関数
f
の呼び出し手順は,yp = f(t, yc, yd, flag, p1, p2,...)
f
はここでも(t, yc, yd, flag, p1, p2,...)
の関すとして関数値を返し,p1, p2,...
は関数パラメータです.- 文字列
Cまたはfortranルーチンの名前を参照する必要があり, <
f_name
> が指定された名前と 仮定します.Fortranの呼び出し手順は以下となります
<f_name>(iflag, nc, nd, t, y, ydp)
double precision t, y(*), ydp(*)
integer iflag, nc, nd
Cの呼び出し手順は以下となります
void <f_name> (int *iflag, int *nc, int *nd, double *t, double *y, double *ydp)
FortranおよびCのケースは共に, 入力引数は以下となります:
iflag
=0
または1
nc
= 連続状態量の数yc
nd
= 離散状態量の数yd
t
= 時間y
=[yc; yd; param]
. パラメータは odedcコール(y = odedc([y0c; y0d], nd, stdel, t0, t, list('fexcd', param)))
で指定された外部引数を取得する際に使用できます.出力
ydp
と同様に, このルーチンはydp[0:nc-1] = d/dt ( yc(t) )
(iflag=0
の場合) およびydp[0:nd-1] = yd(t+)
(iflag=1
の場合) を計算する必要があります.
説明
y = odedc([y0c;y0d], nd, [h,delta], t0, t, f)
は,
離散/連続混合系の解を計算します.
離散系の状態量 ydk
は以下のように
区間内は定数となる時間関数yd(t)
として
組み込まれます:
以下の t
について yd(t) = ydk
[tk = h*(delta+k) , tk+1 = h*(delta+k+1) )
.
シミュレーションされる方程式は以下のようになります for t in [ tk , tk+1 ):
dyc/dt = f(t, yc(t), yd(t), 0) yc(t0) = y0c
そして, 時刻 tk
において
離散変数 yd
が以下のように更新されます:
yd(tk+)
= f(yc(tk-), yd(tk-), 1)
.
yd(t)
の定義に基づき,
最後の式から以下が導かれることに注意してください
ydk
= f(tk, yc(tk-), yd(tk-1), 1)
(yc
は時間連続:
yc(tk-) = yc(tk)
).
f
をコールする際のパラメータは以下のように固定されています:
ycd = f(t,yc,yd,flag)
;
この関数は, flag=0
の場合に
微係数ベクトル yc
,
flag=1
の場合に
yd
の更新値を返す必要があります.
ycd = dot(yc)
は,
flag=0
の場合にはyc
と同じ次元の
ベクトルである必要があり,
flag=1
の場合にはyd
と同じ次元の
ベクトルである必要があります.
t
は解y
が計算される時刻を
表すベクトルです.
y
は
ベクトルy = [y(t(1)), y(t(2)),...]
です.
この関数は ode
関数と同じオプションパラメータを指定して
コールすることができます
(nd
および stdel
はコール時に2番目および3番目のパラメータとして指定されます).
特に積分フラグ,許容誤差を設定することができます.
オプションのパラメータは odeoptions
関数により
設定できます.
外部ルーチンは動的にリンクすることが可能です (link参照).
例
//切替入力を有する線形システム deff('xdu=phis(t,x,u,flag)','if flag==0 then xdu=A*x+B*u; else xdu=1-u;end'); x0=[1;1];A=[-1,2;-2,-1];B=[1;2];u=0;nu=1;stdel=[1,0];u0=0;t=0:0.05:10; xu=odedc([x0;u0],nu,stdel,0,t,phis);x=xu(1:2,:);u=xu(3,:); nx=2; plot2d(t, x', [1:nx], '161') plot2d2('onn',t',u',[nx+1:nx+nu],'000'); //Fortran外部関数 (fydot2.f参照): norm(xu-odedc([x0;u0],nu,stdel,0,t,'phis'),1) //サンプル値フィードバック // // | xcdot=fc(t,xc,u) // (システム) | // | y=hc(t,xc) // // // | xd+=fd(xd,y) // (フィードバック) | // | u=hd(t,xd) // deff('xcd=f(t,xc,xd,iflag)',... ['if iflag==0 then ' ' xcd=fc(t,xc,e(t)-hd(t,xd));' 'else ' ' xcd=fd(xd,hc(t,xc));' 'end']); A=[-10,2,3;4,-10,6;7,8,-10];B=[1;1;1];C=[1,1,1]; Ad=[1/2,1;0,1/20];Bd=[1;1];Cd=[1,1]; deff('st=e(t)','st=sin(3*t)') deff('xdot=fc(t,x,u)','xdot=A*x+B*u') deff('y=hc(t,x)','y=C*x') deff('xp=fd(x,y)','xp=Ad*x + Bd*y') deff('u=hd(t,x)','u=Cd*x') h=0.1;t0=0;t=0:0.1:2; x0c=[0;0;0];x0d=[0;0];nd=2; xcd=odedc([x0c;x0d],nd,h,t0,t,f); norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd1')) // 計算高速化 (fydot2.f参照) plot2d([t',t',t'],xcd(1:3,:)'); scf(2); plot2d2("gnn",[t',t'],xcd(4:5,:)'); scf(0);
参照
- ode — 常微分方程式ソルバ
- link — dynamic linker
- odeoptions — ODEソルバのオプションを設定する
- csim — 線形システムのシミュレーション (時間応答)
- external — Scilabオブジェクト, 外部関数またはルーチン
Report an issue | ||
<< ode_root | 微分方程式 | odeoptions >> |