関数が定義できている場合は、擬似ファイル演算子"+"
と累次評価演算子(,)
を利用し、台形近似を用いた積分値を計算してプロットすることができます。(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行目で現ステップのxとf(x)の値は次のステップで積分智の積算に用いるため変数lastxとlastyにそれぞれ保存しておきます。最後にintegを評価してプロットします。
このスクリプトの実行例は以下のようになります。下図のように、Gauss関数(exp(-x^2);赤線)をx0からxまで積分した値が緑線でプロットされています。区間の右端では積分値がほぼ√π=1.77245になっていることがわかります。
積分した結果をファイルに保存したい場合は、微分のときと同様にset table "***"
を使います。上の例と比べてplot
がsplot
になっていて、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)
もプロットしています。計算値(緑)と解析的な積分値(青の細線)がピッタリ一致していることがお分かり頂けると思います。