[基]数値積分と乱数（第六章）

私は数学が苦手なのです。


Rubyは、我々やMathematicaのような積分のテクニックをもっていません。
なので、台形公式、Simpson公式、Monte Carlo法などで近似的に数値積分をおこないます。
わかっておかなければならないのは、これらは近似にすぎず、誤差が場合によってはかなり生まれてしまうことです。


台形公式

この数値積分を行う関数をtrapezoid(xs,xe,n)とすると、定義はこうなります。
def  trapezoid(xs,xe,n)
  deltax = (xe-xs)*1.0 / n
  sum = (f(xs)+f(xe)) / 2.0
  for  i  in 1..(n-1)
    sum = sum + f(xs+i*deltax)
  end
  deltax * sum
end
f(x)を別に定義すれば、いろんな関数を台形公式で数値積分できるのです。


Simpson公式
こちらは、より正確にするため、直線ではなく二次曲線で近似しようぜ！ってもの。

これをsimpson(xs,xe,n)として定義すると、以下のようになります。
def  simpson(xs,xe,n)
  deltax=(xe-xs)/(2.0*n)
  sum=f(xs)+4*f(xs+deltax)+f(xe)
  for  i  in 1..(n-1)
    sum=sum+2*(xs+2*i*deltax)+4*f(xs+(2*i+1)*deltax)
  end
  deltax*sum/3.0
end


Monte Carlo法
Monte Carlo法は、乱数を用いた数学的問題の解法の総称で、積分以外の問題にも用いられます。
それらの問題は、主に決定論的問題（積分含む）と非決定論的問題にわけられます。
前者は、厳密な解が存在し、乱数でそれを近似的に求めるもの、後者は、自然科学・社会科学におけるシミュレーションのように、確率的な変動を含む、厳密な解が存在しないものを指します。
積分において、Monte Carlo法は、「積分区間、その区間におけるf(x)の値域が、0以上有限値以下におさまっている高次元の関数」に対して使われます。


例として円周率を求めてみましょう。

まあこうやって、ランダムにn個点をうちまくって、そのうち何割が関数より下（または上）にくるのかな？って感じで、積分する・・・ていうか、面積を出します。
で、上の場合の積分を行う関数montecarlo(n)はこちら。
def  montecarlo(n)
  m=0
  for  i  in 1..n
    x=rand()
    y=rand()
    if  x*x+y*y < 1.0
      m=m+1
    end
  end
  m*1.0/n
end


rand()が乱数で、[0,1)の数をランダムに出力します。
乱数が含まれているので、montecarlo(n)は、nを同じ数にしても、毎回違う値になります。
そして、最終的に円周率を出すにはこうですね
def  pi(n)
  4.0*montecarlo(n)
end


ちなみに、参考までに、使い回しできるmonte carlo法の関数はこちら。
積分区間[xs,xe]、その区間のf(x)の最大値がmaxの場合の関数です。
また、xs≧0かつ、f(x)≧0（xs≦x≦xe）でなくては使えません。
def  montecurlo_super(xs,xe,max,n)
  m = 0
  for  i  in 1..n
    x = rand() * (xe-xs)
    y = rand() * max
    if  y<f(x+xs)
      m=m+1
    end
  end
  m*(xe-xs)*max*1.0/n
end