gnuplotを用いた関数の積分

関数が定義できている場合は、擬似ファイル演算子"+"累次評価演算子(,)を利用し、台形近似を用いた積分値を計算してプロットすることができます。(gnuplot4.3以降が必要)

以下に一つの例を示します。この例ではGauss積分:

を計算しています。

f(x)=exp(-x**2)  # 被積分関数

x0=-30   # 積分区間の下限
x1= 30   # 積分区間の上限
N=1000   # 積分区間の分割数

set xrange [x0:x1]
set samples N+1

plot lastx=0.0,lasty=0.0,integ = 0.0,\
     "+" using 1:(f($1)) w l lw 3 title "gauss function exp(-x^2)", \
     "+" using 1:(dx=$1-lastx, \
     integ = ($0==0 ? 0.0 : integ+dx*(f($1)+lasty)*0.5), \
     lastx=$1, \
     lasty=f($1), \
     integ) w l lw 3 title "integral"

plotコマンドの4行目で積分値integが台形近似で積算されています。5-6行目で現ステップのxf(x)の値は次のステップで積分智の積算に用いるため変数lastxlastyにそれぞれ保存しておきます。最後にintegを評価してプロットします。

このスクリプトの実行例は以下のようになります。下図のように、Gauss関数(exp(-x^2);赤線)をx0からxまで積分した値が緑線でプロットされています。区間の右端では積分値がほぼ√π=1.77245になっていることがわかります。

積分した結果をファイルに保存したい場合は、微分のときと同様にset table "***"を使います。上の例と比べてplotsplotになっていて、usingの中身も少し違っているのに注意してください。

f(x)=exp(-x**2)  # 被積分関数

x0=-30   # 積分区間の下限
x1= 30   # 積分区間の上限
N=1000   # 積分区間の分割数

set xrange [x0:x1]
set samples N+1

set table "integ_func_output.txt"
splot lastx=0.0,lasty=0.0,integ = 0.0,\
     "+" using 1:(f($1)):(dx=$1-lastx, \
     integ = ($0==0 ? 0.0 : integ+dx*(f($1)+lasty)*0.5), \
     lastx=$1, \
     lasty=f($1), \
     integ)
unset table

この例の出力ファイルはこちらです。

被積分関数を変える場合はf(x),x0,x1,Nの値を適宜変更したらOKです。以下にもう一つの例を示します。今度はローレンツ関数の積分:

を行います。

f(x)=1.0/(1.0+x**2)  # 被積分関数

x0=-30   # 積分区間の下限
x1= 30   # 積分区間の上限
N=1000   # 積分区間の分割数

set xrange [x0:x1]
set samples N+1
set key left reverse Left

plot lastx=0.0,lasty=0.0,integ = 0.0,\
     "+" using 1:(f($1)) w l lw 3 title "Lorentz function 1/(1+x^2)", \
     "+" using 1:(dx=$1-lastx, \
     integ = ($0==0 ? 0.0 : integ+dx*(f($1)+lasty)*0.5), \
     lastx=$1, \
     lasty=f($1), \
     integ) w l lw 3 title "integral",\
     atan(x)-atan(x0) w l lw 1

この例では比較のために解析的に求めた積分値atan(x)-atan(x0)もプロットしています。計算値(緑)と解析的な積分値(青の細線)がピッタリ一致していることがお分かり頂けると思います。