RubyのRangeを積分範囲としてブロックで与えられた関数を数値積分するライブラリ　使用例: (0.0..Math::PI).simpson(20){|x| Math::sin(x)}

00README.md

RubyのRangeを積分範囲としてブロックで与えられた関数を数値積分するライブラリ

下に置いてある短いライブラリ https://gist.github.com/t-nissie/b6ef8d39229a2534498b/raw/integrations.rb はRubyのRangeを積分範囲としてブロックで与えられた関数を数値積分します。 台形公式か シンプソンの公式か が選べます。

使い方

Rangeのメソッドになってます。 引数として分割数を ブロックとして関数を 与えます。 https://gist.github.com/t-nissie/b6ef8d39229a2534498b/raw/integrations.rb の使い方は次の通りです。

#!/usr/bin/env ruby
##
require './integrations.rb'
p (0.0..Math::PI  ).trapezoidal(2000){|x| Math::sin(x)}
p (0.0..Math::PI  ).simpson(    2000){|x| Math::sin(x)}

https://gist.github.com/t-nissie/b6ef8d39229a2534498b/raw/integrations.rb の後半部分や https://gist.github.com/t-nissie/b6ef8d39229a2534498b/raw/part.rb なども参照のこと。

例の積分

数値積分の値の再現

2018年ごろに少し話題になっていた岩波数学公式Iのp.240の数値積分に挑戦する。

朝日新聞1988年3月23日の科学欄を見ると、 阪大で物性物理を研究していた斎藤基彦さんが公式の間違いを見つけたらしい。 この「公式」は岩波の数学公式や丸善の数学大公式集にも間違いのまま掲載されていたが、 元は1867年のアムステルダムの公式集で、さらに1790年の本が大元だったようだ。 ［朝日新聞の新聞記事のコピーは略。元tweetを参照のこと。］

— Paul Painlevé@JPN (@Paul_Painleve) April 19, 2018

$$ \int_{0}^{1} \frac{1}{\log|\log x|} dx = -0.154479641320 $$

という定積分なのだが、曲者。被積分関数を $f(x) = \frac{1}{\log|\log x|}$ としてプロットすると、

となって、 $1/e$ で発散しているので直接数値積分はできない。 座標 $(1/e,0)$ の点のまわりで回転対称に近いので、 $x=1/e$ で $f(x)$ をひっくり返したもの $f(2/e-x)$ を足して、

$$ \int_{0}^{1/e} [f(x)+f(2/e-x)] dx + \int_{2/e}^{1} f(x) dx $$

がおそらく求める積分。 ただしこれも一筋縄ではいかず、 $x=0$ , $x=1$ 近傍で被積分関数は大きく変化するので数値積分の分割を細かくする必要がありそうだし、 $x=1/e$ 近傍では「大きな数−大きな数」になっているので、 区間 $(\frac{1}{e}-10^{-2}, \frac{1}{e}+10^{-2})$ では級数展開

$$f(x)+f(2/e-x) = -\frac{e^2}{12}(x-1/e)^2 - \frac{2e^4}{45}(x-1/e)^4 -\frac{859e^6}{30240}(x-1/e)^6 + O((x-1/e)^8)$$

のほうを数値積分する。級数展開は https://wolframalpha.com が教えてくれた。自分では導けていない。 $(\frac{1}{e}-10^{-2}-10^{-3},\frac{1}{e}-10^{-2}+10^{-3})$ での $f(x)+f(2/e-x)$ と級数展開との差が下のexpansion.png。

上に書いたようにx=0とx=1の近傍だけ分割を細かくして数値積分しているのが https://gist.github.com/t-nissie/b6ef8d39229a2534498b/raw/part.rb

#!/usr/bin/env ruby
##
include Math
require './integrations.rb'
def f(x)
  1.0/log(-log(x))
end
def g(x)
  (x-1/E).abs<1e-2 ? -(x-1/E)**2*E*E/12 -2*(x-1/E)**4*E**4/45 -859*(x-1/E)**6*E**6/30240: f(x)+f(2/E-x)
end
Fmt   = "%20.17f\n"
printf Fmt, i2=(2/E..1.0-1e-2).simpson(    10000){|x| f(x)}
printf Fmt, i3=(1.0-1e-2..1.0).simpson(100000000){|x| f(x)}
printf Fmt, i4=(0.0..1e-4).simpson(    100000000){|x| g(x)}
printf Fmt, i5=(1e-4..1/E).simpson(       100000){|x| g(x)}
puts '--------------------'
printf Fmt, i2+i3+i4+i5
puts '===================='
printf "%20.17f %s\n", -0.154479641320, 'Iwanami'

Ruby 3.1.2p20と2022年ごろの最新のマシンの１コアで35秒くらいかかる。

$ /usr/bin/time ruby --jit part.rb
-0.13449639088600002
-0.00183072931408726
-0.00004147904805336
-0.01811104207072666
--------------------
-0.15447964131886729
====================
-0.15447964132000000 Iwanami
       32.25 real        34.35 user         0.24 sys

変数変換とかもっとエレガントな方法があるのかも。

底を調整して定積分を０にする

底を $b$ として

$$ \int_{0}^{1} \frac{1}{\log|\log_b x|} dx = 0 $$

を $b$ について解きたい。 $x=1/b$ において主値をとる超曲者。 だいたい $b=2.317$ 。 高速かつ正確に数値積分ができるならニュートン法か何かで $b$ を決めることができるが、 https://wolframalpha.com 先生をもってしても $x=1/b$ のまわりの級数展開がわからない。 $g(x) = \frac{1}{\log|\log_b x|}$ として

$$g(x)+g(2/b-x) = 1 - \log b - c_2(x-1/b)^2 - c_4(x-1/b)^4 - c_6 (x-1/b)^6 + O((x-1/b)^8)$$

らしいが、なぜ $x=1/b$ で $1 - \log b$ なのか？ $c_2, c_4, c_6$ は解らないのか？

まとめ

とりあえずRubyで積分をやってみた。 Rangeのメソッドになっていてかっこよいのだが大した機能もないし精度保証もない。 実はJuliaなら既存のパッケージで簡単に精度を保証しながら数値積分ができる → https://gist.github.com/t-nissie/cd1bd788db8aad647566653171c2dfa8 （『岩波数学公式I』の例の数値積分をJuliaのQuadGKパッケージを使って精度保証する）

expansion.gp

#!/usr/bin/env gnuplot
# Author: Takeshi NISHIMATSU
# Usage: gnuplot expansion.gp
##
set terminal png
set output "expansion.png"
set key left top
set samples 2000
s=1e-3
set xrange [exp(-1)-1e-2-s:exp(-1)-1e-2+s]
set xtics exp(-1)-1e-2-s, s

f(x)= 1.0/log(-log(x))
expansion(x)=-(x-exp(-1))**2*exp(2)/12 -2*(x-exp(-1))**4*exp(4)/45 -859*(x-exp(-1))**6*exp(6)/30240

plot f(x)+f(2*exp(-1)-x)-expansion(x)

set output

expansion.png

integrations.rb

#!/usr/bin/env ruby
# integrations.rb is a collection of numerical integration methods
# Gist: https://gist.github.com/t-nissie/b6ef8d39229a2534498b
# Author: Takeshi Nishimatsu
# Copying:
#    Copyright (c) 2016 by Takeshi Nishimatsu
#    integrations.rb is distributed in the hope that
#    it will be useful, but WITHOUT ANY WARRANTY.
#    You can copy, modify and redistribute integrations.rb,
#    but only under the conditions described in
#    the GNU General Public License version 3 (the "GPLv3").
##
class Range
  def trapezoidal(n=10)
    n  = n.to_i
    x0 = self.first.to_f
    xn = self.last.to_f
    h  = (xn - x0) / n
    sum =  yield(x0)/2
    (1...n).each{|i| sum += yield(x0+i*h)}
    sum += yield(xn)/2
    sum*h
  end
  def simpson(n=10)
    m  = n.to_i / 2
    x0 = self.first.to_f
    xn = self.last.to_f
    h  = (xn - x0) / m / 2
    sum   =    yield(x0)
    (1...m).each{|i|
      sum += 4*yield(x0+(2*i-1)*h)
      sum += 2*yield(x0+ 2*i   *h)
    }
    sum   += 4*yield(xn-h)
    sum   +=   yield(xn)
    sum*h/3
  end
end
if $0 == __FILE__
  p (0.0..Math::PI  ).trapezoidal(    ){|x| Math::sin(x)}
  p (0.0..Math::PI  ).trapezoidal('10'){|x| Math::sin(x)}
  p (0.0..Math::PI  ).trapezoidal(2000){|x| Math::sin(x)}
  p (0.0..Math::PI  ).simpson(    2000){|x| Math::sin(x)}
  p (0.0..Math::PI  ).simpson(    2001){|x| Math::sin(x)}
  p (0.0..Math::PI/2).simpson(    2000){|x| Math::cos(x)}
end

logabslogx5-5.gp

#!/usr/bin/env gnuplot
# Time-stamp: <2022-09-19 19:10:16 takeshi>
# Author: Takeshi NISHIMATSU
# Usage: gnuplot logabslogx5-5.gp
##
set terminal postscript eps enhanced color dashed "Times-Roman" 20
set encoding iso_8859_1   # for longer minus
set output "logabslogx5-5.eps"
#set key left
set xrange [-1/exp(1):1]
set xtics ("-1/e" -1/exp(1), "0" 0, "1/e" 1/exp(1), "2/e" 2/exp(1), "1" 1)
m=5
set yrange [-m:m]
set ytics 1
set samples 10000
set grid
f(x)   = 1/log(-log(x))
plot f(x) w l lw 4, -f(2/exp(1)-x) t "-f(2/e-x)" w l lw 4, f(x)+f(2/exp(1)-x) t "f(x)+f(2/e-x)" w l lw 4
set output
#!epstopdf "logabslogx5-5.eps"
#Local variables:
#  compile-command: "gnuplot logabslogx5-5.gp"
#End:

logabslogx5-5.png

loss.gp

#!/usr/bin/env gnuplot
# Author: Takeshi NISHIMATSU
# Usage: gnuplot loss.gp
##
set terminal postscript eps enhanced color dashed "Times-Roman" 16
set encoding iso_8859_1   # for longer minus
set output "loss.eps"
set key right bottom spacing 2
set xrange [1/exp(1)-3e-4:1/exp(1)+3e-4]
set xtics ("1/e-3{/Symbol \264}10^{-4}" 1/exp(1)-3e-4,\
           "1/e-2{/Symbol \264}10^{-4}" 1/exp(1)-2e-4,\
           "1/e-1{/Symbol \264}10^{-4}" 1/exp(1)-1e-4,\
           "1/e" 1/exp(1),\
           "1/e+1{/Symbol \264}10^{-4}" 1/exp(1)+1e-4,\
           "1/e+2{/Symbol \264}10^{-4}" 1/exp(1)+2e-4,\
           "1/e+3{/Symbol \264}10^{-4}" 1/exp(1)+3e-4)
m=1e-7
set yrange [-m:m]
set samples 50000
#set grid
f(x)   = 1/log(-log(x))
plot f(x)+f(2/exp(1)-x)        t "1/log |log x|+1/log |log(2/e-x)|" w l lw 1,\
     -(x-exp(-1))**2*exp(2)/12 t "-e^2(x-1/e)^2/12" w l lw 2
set output
#!epstopdf "loss.eps"
#Local variables:
#  compile-command: "gnuplot loss.gp"
#End:

loss.png

part.rb

#!/usr/bin/env ruby
##
include Math
require './integrations.rb'
def f(x)
  1.0/log(-log(x))
end
def g(x)
  (x-1/E).abs<1e-2 ? -(x-1/E)**2*E*E/12 -2*(x-1/E)**4*E**4/45 -859*(x-1/E)**6*E**6/30240: f(x)+f(2/E-x)
end
Fmt   = "%20.17f\n"
printf Fmt, i2=(2/E..1.0-1e-2).simpson(    10000){|x| f(x)}
printf Fmt, i3=(1.0-1e-2..1.0).simpson(100000000){|x| f(x)}
printf Fmt, i4=(0.0..1e-4).simpson(    100000000){|x| g(x)}
printf Fmt, i5=(1e-4..1/E).simpson(       100000){|x| g(x)}
puts '--------------------'
printf Fmt, i2+i3+i4+i5
puts '===================='
printf "%20.17f %s\n", -0.154479641320, 'Iwanami'