lsoda の和訳 (未完) - ドレッシングのような

lsoda: ODE (常微分方程式) の求解系

Description

1階常微分方程式 (ODEs) の硬い系もしくは柔らかい系について初期値問題を解くことについて，R の関数 lsoda は，同じ名前の Fortran ODE ソルバへのインタフェースを提供する．Fortran 版の lsoda は Linda R. Petzold と Alan C. Hindmarsh によって書かれた．ODE の系は R の関数 (もちろん，外部コードの呼び出しのために .C，.Fortran，.Call などを使用できる) として記述されるか，動的にロードされるコンパイル済みコードで定義される．パラメータのベクトルは ODE に渡される．そのため，ソルバを ODE についてのモデル化パッケージの一部として使用できる．また，optim，nls，nlm，そして nlme のような，R の非線形モデルのための任意の適切なモデル化ツールを使用することで，ソルバをパラメータ推定のためにも使用できる．

Usage

lsoda(y, times, func, parms, rtol, atol, tcrit=NULL, jacfunc=NULL, verbose=FALSE, dllname=NULL, hmin=0, hmax=Inf, ...)

Arguments

`y`

ODE 系の初期値．もし y が名前属性を持つならば，その名前は出力行列のラベルとして使用される．

`times`

y について明示的に推定される時刻．times の第一要素は初期時刻でなければならない．

`func`

時刻 t で ODE 系 (モデル定義) の微分係数を計算するユーザ供給関数，もしくは動的にロードされる共有ライブラリのコンパイル済関数名を与える文字列のどちらか．

もし func がユーザ供給関数であれば，yprime = func(t, y, parms) として呼び出せなければならない．t は積分の現在時点，y は ODE 系の変数が持つ現在の推定値である．もし初期値 y が名前属性を持てば，その名前は func 内でも有効である．parms はパラメータのベクトルである (これも名前属性を持てる，大きな系では望まれる)．

func の返り値はリストであるべきで，その第一要素はt に関する y の微分を含むベクトル，第二要素は times の各点で必要となる大域的な値のベクトル (名前属性を持てる) である．

もし func が文字列なら，dllname は (拡張子を除いた) 共有ライブラリ名を与える必要がある．その共有ライブラリは lsoda() が呼び出される前にロードされている必要がある．より多くの情報は Details を参照すること．

`parms`

func 内で使用される任意のパラメータで，その関数を書き換えること無く変更できるべきである．

`rtol`

相対許容誤差．スカラか，もしくは y と同じ長さの配列である．Details を見よ．

`atol`

絶対許容誤差．スカラか，もしくは y と同じ長さの配列である．Details を見よ．

`tcrit`

Fortran ルーチンとしての lsoda は，そのターゲット (times ベクトルの時刻点) を飛び越え，そして望まれる時刻点に対して値を補間する．もし，(たぶん特異点などによって) 積分すべきでない時刻があるならば，それは tcrit で提供されるべきである．ソルバは，tcrit よりも早い times の最終点で停止して戻るため，tcrit を通り越してtimes の時刻を含むことが (誤りでないにもかかわらず) 意味を成さない事に注意すること．

`jacfunc`

もし NULL でなければ，微分方程式系 dydot(i)/dy(j) の Jacobian を計算する R の関数，もしくは Jacobian を計算する `dllname' の関数名またはサブルーチン名を与える文字列である (このオプションの詳細について Details を見よ)．いくつかの事情で，系が硬いならば，jacfunc を供給することは計算を高速化できる．R の jacfunc に対する呼び出しシーケンスは func に対するものと同じである．jacfunc は，第((i - 1)*length(y) + j)番目の値が dydot(i)/dy(j) であるようなベクトルを返さなければならない．すなわち，行列 dydot/dy を返す．ここで第 i 番目の行は $dy_i/dt$ の $y_j$ ，すなわち (R および Fortran が行列を保持する方法である) カラムによる微分である．

`verbose`

論理値であり，TRUE のときは ODE ソルバからより詳細な出力がなされるべきである．現在は何もしない．

`dllname`

(拡張子を除いた) 共有ライブラリ名を与える文字列である．これは，func と jacfunc で参照されるすべてのコンパイル済み関数を含んでいる．

`hmin`

積分ステップサイズの最適な最小値．特殊な状況では，このパラメータは精度のコストと共に計算を高速化する．もしあなたがこの理由を知らないのであれば，hminを使用しないこと！

`hmax`

積分ステップサイズの最適な最大値．(外部入力を伴う) 非自律系について，最大値が必要となる可能性がある．その他の場合では，シミュレーションが短い外部事象を無視する可能性がある．

`...`

追加引数．これによって，この関数が一般化される．

Details

すべての難しい仕事は Fortran のサブルーチンである lsoda によって行なわれる．その文書は詳細のために参照されるべきである (ソースファイル `src/lsoda.f' のコメントとして含まれている)．これは，Netlib の 1997年2月24日版の lsoda に基づいている．以下の誤差制御についての説明は，その文書から改変したものである (上記の入力引数 rtol と atol)：

入力パラメータ rtol と atol はソルバによって処理される誤差制御を決定する．ソルバは， $\mathbf{y}$ の推定局所誤差であるベクトル $\mathbf{e}$ を， $(\mathbf{e}/\mathbf{ewt})$ の最大ノルム ${} \leq 1$ という形の不等式に従って制御する．ここで， $\mathbf{ewt}$ は正の誤差荷重ベクトルである．rtol と atol の値は，すべてが非負となるべきである． $\mathbf{ewt}$ の形は以下で与えられる：

$\mathbf{rtol} \times \mathrm{abs}(\mathbf{y}) + \mathbf{atol}$

ここで，ふたつのベクトルの乗算は要素-対-要素で行う．

もし精度の要求がマシンの能力を越えているならば，Fortran のサブルーチン lsoda はエラーコードを返す．いくつかの条件の下で，R の関数 lsoda は，解を得るために，リーズナブルな精度の還元を試みる．もしそれが行われるならば，警告を出力する．

モデルは，R のものと同様に，C もしくは Fortran のコンパイル済コードで定義できる．C に対して，func の呼び出しシーケンスは以下の例のようにならなければならない：

void myderivs(int *neq, double *t, double *y, double *ydot)
{
    ydot[0] = -k1*y[0] + k2*y[1]*y[2];
    ydot[1] = k3 * y[1]*y[1];
    ydot[2] = -ydot[0]-ydot[2];
}

ここで，*neq は方程式の数，*t は独立変数の値，y は状態変数の現在値を持つ長さ *neq の倍精度配列へのポインタ，そして ydot は計算された微分値を保持する配列へのポインタである．

この例では，パラメータは C コードのグローバル変数として保持され，以下のように宣言される．

static double params[3];

#define k1 params[0] のようにすることで，#define ステートメントがコードをより読み易くする．

これは，R のコードから呼び出されるコンパイル済 C 関数へパラメータを渡す唯一の方法である．このメカニズムを用いる関数は，共有ライブラリと同じ名前を持つ関数を付随しなければならない．その関数は，単一の関数ポインタ引数を持つ (以下の宣言を見よ)．その引数で渡される関数は，パラメータの値をグローバル変数へコピーするために，倍精度値で倍精度配列を埋める．次の例では，ライブラリは `mymod.so' と名付けられ，次の関数はパラメータベクトルの初期化のために要求される：

void mymod(void (* odeparms)(int *, double *))
{
    int N=3;
    odeparms(&N, parms);
}

ここで mymod は odeparms を呼び出すだけである．odeparms は，パラメータベクトルの次元を持つ int へのポインタと，パラメータ値を持つ配列へのポインタを伴って呼び出される．

Model は Fortan でも定義できる．例えば：

subroutine myderivs (neq, t, y, ydot)
double precision t, y, ydot, parms(3)
integer neq
dimension y(3), ydot(3)
common /myparms/parms

ydot(1) = -parms(1)*y(1) + parms(2)*y(2)*y(3)
ydot(3) = parms(3)*y(2)*y(2)
ydot(2) = -ydot(1) - ydot(2)

return
end

Fortran では，パラメータは共有ブロックに保持できる．この場合，モデル関数定義を持つファイルは，そのファイルと同じ名前のサブルーチンも持たなければならない：

subroutine mymod(odeparms)
external odeparms
integer N
double precision parms(3)
common /myparms/parms

N = 3
call odeparms(N, parms)
return
end

モデルがコンパイル済コードで定義されているとき，ODE の直接的な解として返される量の供給は無い (R のコードで定義されたモデルとは異なる)．

(以下未訳のため省略)