積分の数値計算

忘れないうちにメモ. 台形則により関数 f(x) の定積分(広義積分を含む)を計算するときに収束しやすい形にする方法.

1. 積分区間を [-1,1] にするように置換積分を施し, ∫_[y=-1,1]g(y)dyなる形にする.
   • 積分範囲が閉区間[a,b]であるものは y=(2x-a-b)/(b-a), dy=2dx/(b-a) と置換積分.
   • 積分範囲が半開区間[0,+∞)であるものは x=(1+y)/(1-y), dx=2dy/(1-y)² と置換積分.
   • 積分範囲が両開区間(-∞,+∞)であるものは x=y/(1-y²), dx=(1+y²)dy/(1-y²)² と置換積分.
2. 端点特異性 (端点における高階微分係数が発散する場合に台形則による積分計算の精度が悪くなる) の問題を避けるため, y=tanh((π/2)sinh(z)), dy=(π/2)cosh(z)dz/(cosh²((π/2)sinh(z))) なる変換をして, ∫_{[z=-∞, ∞]}h(z)dz とする.

定積分I=∫_[x=a,b]f(x)dx の台形則による計算は, I_N=h((f(a)+f(b))/2+Σ_[n=1,N-1]f(a+nh)) により計算 (h=(b-a)/N として, 必要な精度が得られるまで徐々に N を増やす) するが, 両端点が無限大に発散している場合は h を指定してI_h=hΣ_{[n=-∞, ∞]}f(nh) を計算する. 実際には, Σ_{[n=-∞, ∞]}f(nh)=f(0)+Σ_{[j=1, ∞]}f(jh)+Σ_{[k=1, ∞]}f(-kh)
と分解して, 正方向, 負方向それぞれの無限級数の精度が十分である (f(x)の絶対値がそれまでの有限和の絶対値と比べて十分に小さい) と判定されたときに計算を打ち切る.

例えば, ∫[x=-1,1]2√(1-x²)dx は ∫[z=-∞,∞]πcosh(z)/cosh³((π/2)sinh(z))dz
と変形することで台形則にとっては都合がよくなる. もっとも, この積分は円周率を計算したいがために数値積分をしているのに, 式の変形後の被積分関数に円周率が出てくるのでは無意味だという意見が出るのは仕方がないことではあるが.