Title : Math::ODE
|
Location : Home > Languages > Perl > Package Title : Math::ODE |
|
Math::ODE - N-次の常微分方程式を解く
本モジュールはできるだけ少ない労力で N-次の常微分方程式を解く。現在では IVP（初期値問題）のみをサポートしているが、 BVP（境界値問題）は近い将来追加されるであろう。N-次の方程式を解くために、MATLAB における処理のように最初のターンでN+1次の方程式を解かねばならない。
use Math::ODE;
# ファイルにデータを格納する新しいオブジェクトを生成し、
# 区間 [0,10]、初期値 y(t0) = 0 の条件で所与の方程式を解く。
my $o = new Math::ODE ( file => '/home/user/data',
step => 0.1,
initial => [0],
DE => [ \&DE1 ],
t0 => 1,
tf => 10 );
$o->evolve();
# 方程式 y' = 1/$t を解く。
# $t は独立変数のスカラー
# $y は従属変数の配列参照
sub DE1 { my ($t,$y) = @_; return 1/$t; }
$o->evolve()
$o->t0 から $o->tf まで４次のルンゲ＝クッタ法を利用して方程式を解く。
# Example 1: y'' + y = 0, y(0) = 0, y'(0) = 1 を解く。
# 解は: y = sin(x)
use Math::ODE;
my $o = new Math::ODE ( file => 'data',
verbose => 1,
step => 0.1,
initial => [0,1],
DE => [ \&DE0, \&DE1 ],
t0 => 0,
tf => 10 );
$o->evolve;
# gnuplot でデータをプロットするなら: plot 'data' using 1:2, sin(x)
# y'' + y = 0
sub DE0 { my ($t,$y) = @_; return $y->[1]; }
sub DE1 { my ($t,$y) = @_; return -$y->[0]; }
y'' + y = 0 を解くために２つの要素からなるベクトル y を想定する。１つめの要素は y0 = y であり、２つめの要素は y1 = y0' である。これらの変数を用いて方程式を書き換えると y1' + y0 = 0 になる。y1' について解けば y1' = -y0 を得る。ベクトル y' は y0' = y1 及び y1' = -y0 という要素を持つ。これらは配列参照を利用することを除いて DE* サブルーチンで行っていることと同じである。DE0 サブルーチンは y0' に、DE1 サブルーチンは y1' に対応する。これは容易にN次方程式や連立方程式に一般化することができる。２次の連立方程式の例が example/de4 にある。そこにある他の例もよく見てほしい。
initial 配列参照は従属変数ベクトルの各要素の初期条件に対応している。[0,1] は $o->t0 （たまたま０だが）における最初の要素の値が 0 であり、$o->t0 （これは１０）における２番目の要素の値が 1 であることを意味している。
$o->DE($F)
解くべき方程式を $F に設定する。$F はコード参照の配列参照でなければならない。さもなければマズイことが起こるかも。$F が $o->initial と同じサイズでなければ、あなたは自分を恥じるべきだ。プログラムは終了するだろう。
$o->initial($arrayref)
$arrayref に初期条件を設定する。引数が指定されなければ現在の arrayref の初期条件を返す。$o->DE と同じサイズでなけれなならない。さもなければ自己生体解剖を始める。
$o->step($s)
$s に対するステップサイズ（間隔の幅）を設定する。正の値でなければならない。デフォルトは 0.01 である。
引数が指定されなければ現在のステップサイズを返す。
$o->t0($a)
区間の初期値（左端）の値を $a に設定する。引数が指定されなければ現在の値を返す。
$o->tf($b)
区間の最終値（右端）の値を $b に設定する。引数が指定されなければ現在の値を返す。
$o->file($somefile)
データを $somefile に格納する。引数が指定されなければデータが格納されているファイル名を返す。ファイルが指定されていない場合には標準出力（STDOUT）に出力される。データファイルは、$N を方程式の数として $N + 1 列あり、直接 gnuplot に出力できる。最初の列は独立変数であり、残りは被説明変数ベクトルの n 番目の要素を最初に通過する点である。データファイルからグラフ化する例は配布ソースの example/ ディレクトリにある。
$o->verbose($number)
デバッグのための冗長モードを設定する。デフォルトは 1 で、evolve() を使用した際に NaN または Inf が出力されたのみ場合にのみ警告を発するものである。0 に設定すると undef を返す。shooting method で境界値問題を解くためにコードを書く際には 0 に設定しておいたほうがよいだろう。初期条件を探索する際に、いちいち警告が表示されるのを好みはしないだろうから。
また、 2 に設定すると表示されるメッセージは次のようになる。
Exiting RK4 with t=9.9 ,$y1 = ['-0.544013766248772833','-0.839075464413064726'];
独立変数 $t が増加するたびに、$t に対する４次 Runge-Kutta 法による計算値が出力される。
Jonathan Leto, <jonathan@leto.net>
Boyce, DiPrima "Elementary Differential Equations" 5th Ed.
Orwant, Hietaniemi, Macdonald "Mastering Algorithms with Perl" Ch. 16.
Copyright (c) 2001 by Jonathan Leto. All rights reserved.
本パッケージはフリーソフトウェアであり、Perl 本体と同等の条件で修正／再配布してもよい。
|
|
Updated : 2006/12/13 |