連立常微分方程式解くルンゲクッタ書く

やって意味があるのかはまだ自分でもよくわからないんですけど、研究関係で、連立常微分方程式を解くことになりました。
なんか定性的にでも動きが分かるといいよなーとか思っています。

だいたいこういうところに出てくる常微分方程式というのはヤケにひねくれていて、閉じた形で表せないことが多い。
もうちょっと素直になればみんなに好かれるのになぁと思うんですけど、ツンデレぽいので仕方ありませんね。
でそういうツンデレどもは解析的に解けないので、数値的に計算する。
連立常微分方程式はどう計算していくのかなーとかあんまし知らなくて、皆さんにお聞きしたところオイラー法が簡単で良いらしい。
よっしゃーオイラー法でひねくれた性格を真っ直ぐにしてやるぜ!!!!って書いたんですけど、なんかスゲー勢いで発散していきます。
久々にinfとか見た！

もしかして、この方程式はすげーセンシティブで、誤差がちょっと溜まったら発散する超ひねくれ者なんじゃ、みたく怖くなったので、誤差の小さいルンゲクッタ法で解くことにした。
したはいいんですけど、時間がかかりそうな計算でもないので、Perlで書いて楽することにしたよ!!
実はCPANにMath::RungeKuttaとかあるんですけど無視する。
後から分かったんですけど、上で言ってる発散する問題は、実のところ式そのものが間違っていて、正しいやつ使ってたらオイラー法でも余裕で解けたはず。

作ったオレオレモジュールを使いますと、たとえばLotka-Volterraの微分方程式はこんな感じで解けるようになる。
ちなみにぼくの研究とLotca-Volterraは何の関係もないですね。

use strict;
use RungeKutta::LotkaVolterra;

my $rd = RungeKutta::LotkaVolterra->new
    ( filename => 'result.dat',
      stepsize => 1e-3,              
      repetition => 1e5,            # ステップ数
      initial  => [ 1000, 100 ],    # 初期条件
    );

$rd->calculate();
$rd->print_results();

プロットしたらこんな風になるよ。ロトカヴォルテラは、食うヤツと食われるヤツの多さをモデル化した、こわい式だー。

よく覚えてないですけど、ロトカヴォルテラはわりかしカッコいい式の形で表せたような気がするよ。

上のスクリプトの中にはLotka-Volterraモデルっつーか式なんてどこにも出てないので、あほか、みたく思ったことでしょう。
解くべき方程式は、オレオレモジュールのRungeKuttaを継承したクラスの_initメソッドで定義する感じにしてみました。
式とかパラメータとかは新潟大学様のヤツをお借りしましたよー。

RungeKutta::LotcaVolteraはこんな感じ。エラー処理とかは消しまくった。

package RungeKutta::LotkaVolterra;
use strict;
use base 'RungeKutta';

use constant a => 0.01;
use constant b => 0.05;
use constant c => 0.0001;

sub _init {
    my ($self, %args) = @_;

    # 解く微分方程式はここに書く!!!!!!!!!
    my $eq = [
	      sub {
		  my ($t, $x, $y) = @_;
		  return a * $x - c * $x * $y;
	      },
	      sub {
		  my ($t, $x, $y) = @_;
		  return c * $x * $y - b * $y;
	      }
	     ];

    $self->SUPER::_init( %args,
			 equations => $eq,
			 initial   => $args{initial}
		       );
}

1;

そのスーパクラスのRungeKutta。

package RungeKutta;

use IO::File;
use base 'Class::Accessor::Fast';
use strict;
use Carp;

__PACKAGE__->mk_accessors( qw(fh equations stepsize repetition initial_condition results) );

sub new {
    my ($class, %args) = @_;

    my $self = bless {}, $class;
    $self->_init( %args );
    return $self;
}

sub _init {
    my ($self, %args) = @_;
    
    my $filename = $args{filename};
    my $fh = IO::File->new( "> $filename" )
	or croak "cannot open $filename: $!";
    $self->fh( $fh );
    $self->stepsize( $args{stepsize} );
    $self->repetition( $args{repetition} );
    $self->initial_condition( $args{initial} );
    $self->equations( $args{equations} );
}

sub calculate {
    my $self = shift;

    my @results;
    my @values = @{ $self->initial_condition() };

    for( my $t = 0, my $cnt = 0; $cnt++ <= $self->repetition(); $t += $self->stepsize() ) {
	@values = $self->next_value( $t, \@values, $self->stepsize() );
	push @results, [ $t, @values ];
    }
    $self->results( \@results );
    return \@results;
}

sub print_results {
    my $self = shift;

    my $fh = $self->fh();
    my $format = "%10.5f " x 4;

    my $result = $self->results();
    while( my $e = shift @$result ) {
	printf $fh "$format\n", @$e;
    }
}

sub next_value {
    my ($self, $t, $arguments, $dt) = @_;

    my $f = $self->equations();
    my @num = ( 0 .. @$f - 1 );
    my @args;

    my @k1 = map { $dt * $f->[$_]->( $t, @$arguments ) } @num;

    @args  = map { $arguments->[$_] + $k1[$_] / 2 } @num;
    my @k2 = map { $dt * $f->[$_]->( $t + $dt / 2, @args ) } @num;

    @args  = map { $arguments->[$_] + $k2[$_] / 2 } @num;
    my @k3 = map { $dt * $f->[$_]->( $t + $dt / 2, @args ) } @num;
    
    @args  = map { $arguments->[$_] + $k3[$_] } @num;
    my @k4 = map { $dt * $f->[$_]->( $t, @args ) } @num;

    map { $arguments->[$_] + ( $k1[$_] + 2 * $k2[$_] + 2 * $k3[$_] + $k4[$_] ) } @num;
}

1;