nishiru3の日記

備忘録です。ネットのゴミ。

Perlで移流方程式

PerlPythonみたいなモジュールがないかと思ってましたが、
「PDL」というものを使うと、描画までいけるみたいです。
「いや、pythonで良いじゃん!」というのもありますが・・・。
ソース自体は、リファレンスとか使えば、もうちょっと短くなるかと思います。

どうもPDLだとScilabなんかも使えるみたいです。
結構、更新されてるみたいですね。

公式?
http://pdl.perl.org/


インストール
http://d.hatena.ne.jp/perlcodesample/20121214/1355388288

移流方程式は下記のとおりとしました。

\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} = 0
移流速度uは一定値の線形移流です。

移流項は1次風上、時間積分オイラーの陽積分です。
面白くはないですが、ちょっとしたアルゴリズムを書いてみて、すぐ可視化する場合、
pythonのような使い方もできるようです。あと、CPANにはもっと多くのモジュールがあるようなので、
さらに便利にできるかもしれません。

use strict;
use warnings;
use PDL::Lite;
use PDL::NiceSlice;
use PDL::Graphics::PLplot;
my $dx = 1.0;
my $dt = 1.0;
my $u  = 0.2;
my @f;
my @fini;
my @fn;
my @x;
for my $i(0..99){
    if ($i >10 && $i < 20){
	$f[$i] = 1.0;
    } else {
	$f[$i] = 0.0;  
    }$fn[$i] = $f[$i];
    $fini[$i] = $f[$i];
    $x[$i] = $i;
}
# 100ステップ
for my $n (1..100){
    # 1次風上+eulerの陽解法
    for my $i(1..99){
	$fn[$i] = $f[$i] - $u/$dx*$dt*($f[$i]-$f[$i-1]);
    }
    # 変数の更新
    for my $i(1..99){
  	$f[$i] = $fn[$i];
    }
}
# 数値データの出力
for my $i (0..99){
    print $i.",".$f[$i]."\n";
}
#
# 描画
my $pl = PDL::Graphics::PLplot->new(DEV => 'png', FILE => 'adv.png');
my $x1 = pdl->zeros(100);
for my $i(0..99){
    $x1($i) .= $i;
}
my $y1 = pdl->zeros(100);
my $y2 = pdl->zeros(100);
for my $i(0..99){
    $y1($i) .= $f[$i];
    $y2($i) .= $fini[$i];
}
$pl->xyplot($x1,$y1);
$pl->xyplot($x1,$y2);

$pl->close;

f:id:nishiru3:20141022090245p:plain