Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Please login or create an account
Change language to: English - Français - Português - Русский

Please note that the recommended version of Scilab is 6.0.1. This page might be outdated.
See the recommended documentation of this function

Scilab help >> Differential Equations, Integration > ode

ode

常微分方程式ソルバ

呼び出し手順

y=ode(y0,t0,t,f)
[y,w,iw]=ode([type],y0,t0,t [,rtol [,atol]],f [,jac] [,w,iw])
[y,rd,w,iw]=ode("root",y0,t0,t [,rtol [,atol]],f [,jac],ng,g [,w,iw])
y=ode("discrete",y0,k0,kvect,f)

パラメータ

y0

実数ベクトルまたは行列 (初期条件).

t0

実数スカラー (初期時間).

t

実数ベクトル (解を計算する時間).

f

外部 (関数または文字列またはリスト).

type

以下の文字列のどれか: "adams", "stiff", "rk", "rkf", "fix", "discrete", "roots".

rtol, atol

実数定数またはyと同じ大きさの実数ベクトル.

jac

外部 (関数または文字列またはリスト).

w, iw

実数ベクトル.

ng

整数.

g

外部 (関数または文字列またはリスト).

k0

整数 (初期時間).

kvect

整数ベクトル.

説明

ode は,dy/dt=f(t,y) , y(t0)=y0で定義された 明示的なODEシステムを解くための標準関数です. この関数は,種々のソルバ,特に ODEPACKへのインターフェイスです. 解く問題の型と使用される手法は最初のオプション引数type の値に依存します. typeは,以下の文字列のどれかになります:

<not given>:

パッケージODEPACKのlsoda ソルバが デフォルトでコールされます. このソルバは非スティッフな問題用の予測子/修正子 Adams 法とスティッフな問題用の 後退差分法(BDF)を自動的に選択します. まず非スティッフな手法が使用され,使用する手法を 決定するためにデータを動的にモニタします.

"adams":

これは非スティッフな問題用です. ODEPACKパッケージの lsode ソルバがコールされ, この関数は Adams法を使用します.

"stiff":

これはスティッフな問題用です. ODEPACKパッケージの lsode ソルバがコールされ, BDF法が使用されます.

"rk":

4次の適応型 Runge-Kutta (RK4) 法.

"rkf":

Fehlbergの4次と5次の Runge-Kutta法に基づくShampine および Wattsのプログラム (RKF45) が使用されます. この方法は非スティッフおよびややスティッフな問題用で, 微係数の評価コストが低い場合に適しています. この手法は一般にユーザが高い精度を望む場合には使用するべきではありません.

"fix":

"rkf"と同じソルバですがユーザインターフェイスが非常にシンプルで, rtol および atolパラメータのみを ソルバに指定することが可能です. この手法は使用できる中で最も簡単な手法です.

"root":

解を得る機能を有するODE ソルバ. ODEPACKパッケージの lsodar ソルバが使用されます. 使用されているのはlsoda ソルバを改変したもので, 指定したベクトル関数の根を見つけることができます. 詳細は ode_root のヘルプを参照してください.

"discrete":

離散時間シミュレーション. 詳細は ode_discrete のヘルプを参照してください.

このヘルプでは,標準的な陽のODEシステムに関するode の使用法のみについて記述します.

  • odeの最も簡単なコールは次のようになります: y=ode(y0,t0,t,f) ただし y0 は 初期条件のベクトル,t0 は初期時間, t は解yを計算する時間のベクトル, yは解ベクトルの行列y=[y(t(1)),y(t(2)),...] です.

    入力引数 f は1次微分方程式 dy/dt=f(t,y)のRHSを定義します. これは外部,すなわち指定した構文の関数または指定したコール手続きを 有するFortranサブルーチンまたはC関数の名前(文字列),またはリストです:

    • f が Scilab 関数の場合, その構文は ydot = f(t,y)となります. ただし,tは実数スカラー(時間), y は実数ベクトル (状態量), ydotは実数ベクトル (dy/dt)です.

    • f が文字列の場合, FortranサブルーチンまたはC関数の名前が参照されます. すなわち,ode(y0,t0,t,"fex")がコマンドの 場合,サブルーチンfexがコールされます.

      Fortranルーチンは以下の呼び出し手順とする 必要があります: fex(n,t,y,ydot), ただし n は整数, t は倍精度実数, y と ydot は倍精度ベクトルです.

      C 関数は以下のプロトタイプとする必要があります: fex(int *n,double *t,double *y,double *ydot)

      t は時間, y は 状態量, ydotは状態量の微分 (dy/dt)です.

      This external can be build in a OS independant way using ilib_for_link and dynamically linked to Scilab by the link function.

    • 引数f は以下の構造を有するリストと することもできます: lst=list(realf,u1,u2,...un) ただし, realf は次の構文を有するScilab関数です: ydot = f(t,y,u1,u2,...,un)

      この構文により realfの引数にパラメータを使用することができます.

    関数 f はベクトルではなく p x q の行列を返すことができます. この行列表記により, n=p+q 次のODEシステムdY/dt=F(t,Y) を解きます. ただし, Yp x qの行列です. 初期条件 Y0p x qの行列 である必要があり, odeの結果は p x q(T+1) の行列 [Y(t_0),Y(t_1),...,Y(t_T)] となります.

  • オプションの入力パラメータとして次の 解の誤差を指定できます: 相対および絶対推定誤差の閾値である rtol および atol. y(i)の指定誤差は次のようになります: rtol(i)*abs(y(i))+atol(i)

    そして,状態量の全ての要素に関してこの誤差が小さい場合, 積分が行われます. rtol または atol が一つの定数の場合, rtol(i) および atol(i) はこの定数の値に設定されます. rtol および atolのデフォルト値は 多くのソルバではそれぞれrtol=1.d-5 および atol=1.d-7, "rfk" および "fix"では rtol=1.d-3 および atol=1.d-4です.

  • スティッフな問題の場合, RHS関数のヤコビアンを オプションの引数jacとして 指定する方が良いでしょう. これは,外部,すなわち,指定された構文を有する関数, または指定した呼び出し手順を有するFortranサブルーチンまたは C関数の名前(文字列),またはリストです.

    jac が関数の場合, 構文は J=jac(t,y)とする必要があります.

    ただし, t は実数スカラー (時間)および y 実数ベクトル (状態量). 結果の行列 J は df/dx を評価する必要があります. ただし, J(k,i) =dfk/dxifk は fのk番目の要素です.

    jac は文字列で, FortranサブルーチンまたはC関数の名前を指し, 以下の呼び出し手順となります:

    Fortranの場合:

    subroutine fex(n,t,y,ml,mu,J,nrpd) 
    integer n,ml,mu,nrpd
    double precision t,y(*),J(*)

    Cの場合:

    void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)

    jac(n,t,y,ml,mu,J,nrpd). 多くの場合, ml,muおよび nrpdを参照するべきではありません.

    jac がリストの場合, fと同じ同じ構文が適用されます.

  • オプションの引数 w および iw は積分ルーチンにより返される情報を 保存するベクトルです(詳細はode_optional_output 参照). これらのベクトルがodeのRHSで提供される場合, 積分は前回の停止時と同じパラメータで再開されます.

  • より多くのオプションを%ODEOPTIONS変数 によりODEPACKソルバに指定することができます. odeoptionsを参照ください.

// ---------- Simple one dimension ODE (Scilab function external)
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
function ydot=f(t, y),ydot=y^2-y*sin(t)+cos(t),endfunction
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,f)
plot(t,y)

// ---------- Simple one dimension ODE (C coded external)
ccode=['#include <math.h>'
       'void myode(int *n,double *t,double *y,double *ydot)'
       '{'
       '  ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
       '}']
mputl(ccode,TMPDIR+'/myode.c') //create the C file
ilib_for_link('myode','myode.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');//compile
exec(TMPDIR+'/loader.sce') //incremental linking
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,'myode');

// ---------- Simulation of dx/dt = A x(t) + B u(t) with u(t)=sin(omega*t),
// x0=[1;0]
// solution x(t) desired at t=0.1, 0.2, 0.5 ,1.
// A and u function are passed to RHS function in a list. 
// B and omega are passed as global variables
function xdot=linear(t, x, A, u),xdot=A*x+B*u(t),endfunction
function ut=u(t),ut=sin(omega*t),endfunction
A=[1 1;0 2];B=[1;1];omega=5;
ode([1;0],0,[0.1,0.2,0.5,1],list(linear,A,u))

// ---------- Matrix notation Integration of the Riccati differential equation
// Xdot=A'*X + X*A - X'*B*X + C , X(0)=Identity
// Solution at t=[1,2] 
function Xdot=ric(t, X),Xdot=A'*X+X*A-X'*B*X+C,endfunction  
A=[1,1;0,2]; B=[1,0;0,1]; C=[1,0;0,1];
t0=0;t=0:0.1:%pi;
X=ode(eye(A),0,t,ric)

// ---------- Matrix notation, Computation of exp(A)
A=[1,1;0,2];
function xdot=f(t, x),xdot=A*x;,endfunction 
ode(eye(A),0,1,f)
ode("adams",eye(A),0,1,f)

// ---------- Matrix notation, Computation of exp(A) with stiff matrix, Jacobian given
A=[10,0;0,-1];
function xdot=f(t, x),xdot=A*x,endfunction 
function J=Jacobian(t, y),J=A,endfunction 
ode("stiff",[0;1],0,1,f,Jacobian)

参照

  • ode_discrete — 常微分方程式ソルバ, 離散時間シミュレーション
  • ode_root — 求解付きの常微分方程式ソルバ
  • dassl — 微分代数方程式
  • impl — 微分代数方程式
  • odedc — 離散/連続 ODE ソルバ
  • odeoptions — ODEソルバのオプションを設定する
  • csim — 線形システムのシミュレーション (時間応答)
  • ltitr — 離散時間応答 (状態空間)
  • rtitr — 離散時間応答 (伝達行列)
  • intg — 定積分

参考文献

Alan C. Hindmarsh, lsode and lsodi, two new initial value ordinary differential equation solvers, acm-signum newsletter, vol. 15, no. 4 (1980), pp. 10-11.

使用される関数

関連する関数は SCI/modules/differential_equations/src/fortran directory :

lsode.f lsoda.f lsodar.f です.

Scilab Enterprises
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Mon Oct 01 17:40:22 CEST 2012