風来坊@真幸福知

出雲市 福知寺 行事ご案内 坐禅会 日記 テクノロジー English Zazen session technology Korean 한국어 참선 묵조선 텍놀로지

gnuplotでデータの積分

計算結果のプロットをしていて, 時間積分を見たいことがある.

カラム形式のテキストファイルに書き出したデータから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"

 

f:id:cheonghongsa:20191028204401p:plain