gnuplot で積分

gnuplot積分したい。gnuplot関数型言語であるから、ふつうに再帰で書いてみよう。

f(x)=cos(x)

int_delta(x1,x2)=(f(x1)+f(x2))*(x2-x1)/2
integrate_(x1,x2,n,i)=(i==n?0:integrate_(x1,x2,n,i+1))+int_delta(x1+(x2-x1)*(i-1)/n,x1+(x2-x1)*i/n)
integrate(x1,x2,n)=integrate_(x1,x2,n,1)
g(x)=integrate(0,x,100)

set xrange [-3.1416*2:3.1416*2]
set terminal png
set output "int1.png"
plot f(x) title "cos(x)", g(x) title "integrate 0 -> x, cos(x)"

f:id:crckyl:20150619191823p:plain

上手く動いたようだ。問題は精度である。x = 200 あたりを見てみよう。

set xrange [200:206]
set output "int2.png"
plot f(x) title "cos(x)", \
     g(x) title "integrate 0 -> x, cos(x)", \
     sin(x) title "sin(x)"

f:id:crckyl:20150619192407p:plain

台形近似だからだんだん小さくなってくるよね。分割数を増やしてみよう。

gnuplot> g(x)=integrate(0,x,1000)
gnuplot> set output "int3.png"
gnuplot> replot
         recursion depth limit exceeded

となって詰まっていたのだが、sum [i=1:n] すればよいことに気付いた(gnuplot 4.5 あたりで導入されたらしい)。

integrate2(x1,x2,n)=sum [i=1:n] int_delta(x1+(x2-x1)*(i-1)/n,x1+(x2-x1)*i/n)
g(x)=integrate2(0,x,1000)
set output "int3.png"
replot

f:id:crckyl:20150619192658p:plain