odedc

引数

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 [ t_k , t_k+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オブジェクト, 外部関数またはルーチン