gnuplotを用いた関数の積分（再帰定義）

再帰定義を利用することで、積分を含む関数を定義することができます。この方法を使うと関数が「定義」できるので、積分関数をフィッティングに使える可能性があります。これが累次評価を利用する方法に対する長所です。ただし、再帰呼び出し回数に上限があることに注意しなければなりません。

再帰定義を用いた関数の積分

ある関数func(x)をx0からx1まで積分するための関数integ(x1,x0)は以下のようにして定義することができます。

integ_sub(x1,x0,n,N)=(n <= 1 ? \
     (func(x0)+func(x0+(x1-x0)/N))*(x1-x0)/N*0.5 : \
     integ_sub(x1,x0,n-1,N) + (func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5) 
N=99                              # 区間の分割数
integ(x1,x0)=integ_sub(x1,x0,N,N) # 積分を実行する関数

上記のうち、integ_sub()関数が再帰定義を用いたループの本体です。nが1より大きい場合はnを一つ減らしたinteg_sub(x1,x0,n-1,N)に、グラフの小区間を台形近似した面積(func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5を足して行きます。これをn=1になるまで繰り返すことで、台形の面積の和が得られるわけです。また、Nは区間の分割数を表しており、エラーが出ない範囲でなるべく大きくする方が計算精度は良くなるはずです。

例えば、被積分関数をfunc(x)=sin(x)とした例を以下に示します。

integ_sub(x1,x0,n,N)=(n <= 1 ? \
     (func(x0)+func(x0+(x1-x0)/N))*(x1-x0)/N*0.5 : \
     integ_sub(x1,x0,n-1,N) + (func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5)
N=99                              # 区間の分割数
integ(x1,x0)=integ_sub(x1,x0,N,N) # 積分を実行する関数
func(x)=sin(x)                    # 被積分関数

plot integ(x,0), -cos(x)+1

この例の実行例は以下の通りで、サイン関数の0からxまでの積分integ(x,0)が理論値-cos(x)+1にほぼ等しくなっていることが分かります。

また、plotの部分を以下のように変更して積分の相対誤差をプロットすることができます。

plot (integ(x,0)+cos(x)-1)/(1.0-cos(x))

この実行結果は以下のようになっており、この積分範囲内では誤差は0.1%以下であることが分かります。

しかし、積分区間をあまり大きくしてしまうと積分の計算精度が悪くなります。例えば、上記例のplotコマンドの前にset xrange [100:110]を付け加えると実行結果は以下のようになります。積分と理論値が一致しなくなっていることがお分かり頂けると思います。

積分を含む関数のフィッティング

再帰定義を用いて関数の積分をすると、積分を含む関数をデータに対してフィッティングすることができます。以下ではその一例として通常金属の電気抵抗データへの理論曲線のフィッティングのやり方を示します。

電気抵抗に関する理論によると、金属の電気抵抗の温度依存性は以下の式で与えられます。

ここで、Aは比例係数、ρ0は残留抵抗と呼ばれ試料中の不純物の量などに依存する量、そしてTDはデバイ温度と呼ばれる物質固有のパラメータです。ある金属の電気抵抗の温度依存性のデータがあれば、そこに上記式をフィッティングしてデバイ温度を求めることができるはずです。しかしながら、上の式は右辺に積分を含むため、通常の方法ではgnuplot上でフィッティングをすることができません。そこで、先に述べた再帰定義による積分のテクニックを使います。

関数の定義部分は以下の通りになります。

integ_sub(x1,x0,n,N)=(n <= 1 ? \
     (func(x0)+func(x0+(x1-x0)/N))*(x1-x0)/N*0.5 : \
     integ_sub(x1,x0,n-1,N) + (func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5)
integ(x1,x0)=integ_sub(x1,x0,N,N)
N=98                                                   # 区間の分割数
func(x)=(x==0.0 ? 0.0 : x**5*exp(-x)/(1.0-exp(-x))**2) # 被積分関数
resist(x)=integ(1.0/x,0)*(x**5)                        # 電気抵抗（x=T/TD）

上の例で、resist(x)関数は、電気抵抗の式のうちA、ρ0以外の部分を返す関数になっています。ただしその引数xは温度Tをデバイ温度TDで割ったものにしなくてはなりません。

まず、積分が正しく実行されていることを確かめるために、resist(x)関数を、他のプログラムを使ってまじめに数値積分したデータ（Debye_resistance.dat; ダウンロードはこちら）とともにプロットしてみます。

integ_sub(x1,x0,n,N)=(n <= 1 ? \
     (func(x0)+func(x0+(x1-x0)/N))*(x1-x0)/N*0.5 : \
     integ_sub(x1,x0,n-1,N) + (func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5)
integ(x1,x0)=integ_sub(x1,x0,N,N)
N=98                                                   # 区間の分割数
func(x)=(x==0.0 ? 0.0 : x**5*exp(-x)/(1.0-exp(-x))**2) # 被積分関数
resist(x)=integ(1.0/x,0)*(x**5)                        # 電気抵抗（x=T/TD）

set xrange [0:10]
plot resist(x), "Debye_resistance.dat" u 1:2

すると、以下のようにresist(x)関数は、ほぼまじめな数値積分を再現していることが分かります。

plot部分を以下のように変更して相対誤差をプロットしてみると、相対誤差はだいたい0.01%以下に収まっていることが分かります。

plot "Debye_resistance.dat" u 1:($2-resist($1))/($2)

従って、それほど精度が必要無い用途であれば、resist(x)関数をフィッティングに用いても大丈夫そうです。

フィッテングに使うデータとして、私が昔測定した、白金温度計の電気抵抗の温度依存性データ（Pt.dat；ダウンロードはこちら）を使ってみます。このデータは以下にグラフで示しているようなデータです。横軸は絶対温度（ケルビン）、縦軸が電気抵抗（オーム）になっています。

このデータに対して以下のようにフィッティングをかけて結果を表示してみます。

integ_sub(x1,x0,n,N)=(n <= 1 ? \
     (func(x0)+func(x0+(x1-x0)/N))*(x1-x0)/N*0.5 : \
     integ_sub(x1,x0,n-1,N) + (func(x0+(x1-x0)/N*(n-1))+func(x0+(x1-x0)/N*n))*(x1-x0)/N*0.5)
integ(x1,x0)=integ_sub(x1,x0,N,N)
N=98                                                   # 区間の分割数
func(x)=(x==0.0 ? 0.0 : x**5*exp(-x)/(1.0-exp(-x))**2) # 被積分関数
resist(x)=integ(1.0/x,0)*(x**5)                        # 電気抵抗（x=T/TD）

#fittingパラメータの初期値
TD=200.0
A=100.0
R0=0.9

#fittingの実行
fit A*resist(x/TD)+R0 "Pt.dat" via A,TD,R0

#結果のプロット
plot "Pt.dat", A*resist(x/TD)+R0

上のようにうまくフィッティングができていそうなことが分かります。相対誤差は下のように低温でやや大きくなっています。

plot "Pt.dat" using 1:($2-A*resist($1/TD)-R0)/(A*resist($1/TD)+R0)

ちなみに、フィッティングで得られたデバイ温度は230ケルビンで、白金のデバイ温度の文献値（230ケルビン）と非常に良い一致を示しています。

注意！ 上記例はたまたまうまくいった例ですが、何度も言うように再帰定義を用いた積分は場合によって誤差が非常に大きくなります。このテクニックは十分注意して使って下さい。