計算結果のプロットをしていて, 時間積分を見たいことがある.
カラム形式のテキストファイルに書き出したデータからx軸とy軸のデータのあるカラムを抜き出して, yの値をxに対して数値積分するプログラムを別に作ってパイプにしてやれば良いのだが, gnuplotでもできるらしいというのを見つけてやってみた.
今までやっていたのは, 例えば"data"の第2, 第3カラムがx, yの場合
plo "<integral 2 3 0 50 data" us 1:2 w lp pn 15 ti "Integral-1"
integral <colx> <coly> <xstart> <xend> <data_file>
というのが外部のプログラムで, これでdataのデータを処理した結果を読み込んでプロットするという方式.
このサイトなどを参考に, gnuplotの中で済ませる方法を考えてみた.
gnuplotを20年以上使っていて今までよく知らなかったのだが,
(ほにゃらら, ほにゃらら, ..., y=x)
などのように( )の中で","で区切られた式を連ねておくと, それらを順に計算して最期の計算値が出力となるらしい. それで, 複数の式から成る計算をひとつの関数として定義したり, using ... の中で使えるということらしい.
こういう感じで数値積分(台形法)を関数にしておく.
# numerical integral of colf data in terms of colx data
integ(colx,colf) = (dx=($0==0? 0.:column(colx)-x0), \
fint=($0==0? finti : fint+dx*0.5*(column(colf)+f0)), \
x0=column(colx), f0=column(colf), fint)
# colx : column no. of x data
# colf : column no. of f data
# finiti : initial value of fint: need to be set in advance
ちなみに, $0 (またはcolumn(0))はデータセットの中の現在のデータの順番を返す. "$0==0?"は最初のデータのときに真となる.
上で"integral"を使ってやったのと同じことは, こうすれば良い.
plo finti=0., "data" using 2:(integ(2,3)) w lp pn 15 ti "Integral-2"
