学生向け

数値積分の簡単な方法を教えてください

と訊かれました（もう自分でできるはずなんだけど）。いつもどおり、低級言語で実装すれば数値計算の勉強になります。答えが欲しいだけなら、Maximaのような高級言語を使いましょう（Mathematicaは使える人が少なそうなので、フリー・ソフトウェアのMaximaでまず実装することにしました）。

例として、f(x)=1+2 Pi sin(2 Pi x)を0<=x<=1の範囲で数値積分します。

/* romberg積分 */
(%i1) romberg(sin(x),x,0,1);
(%o1) 0.45969769038987

/* 別の方法（この問題についてはこっちのほうがいい） */
/* Newton-Cotesの8次多項式による求積法（ルーチンは適応型） */

(%i2) load("qq");
(%o2) C:/PROGRA~1/MAXIMA~1.0/share/maxima/5.10.0/share/numeric/qq.lisp
(%i3) quanc8(sin(x),x,0,1);
(%o3) 0.45969769413186

/* もちろんこの問題は、記号的に解けるわけで */
(%i4) integrate(sin(x),x,0,1);
(%o4) 1-cos(1)

/* 数値にすると */
(%i5) bfloat(%);
(%o5) 4.596976941318603b-1

ブラックボックスで価格も高いけど高機能なMathematicaでは次のようになります。

In[1]:= NIntegrate[Sin[x], {x, 0, 1}]

Out[1]= 0.459698

In[2]:= NIntegrate[Sin[x], {x, 0, 1},Method->Trapezoidal]

NIntegrate::slwcon:
   Numerical integration converging too slowly; suspect one of the following:
    singularity, value of the integration being 0, oscillatory integrand, or
    insufficient WorkingPrecision. If your integrand is oscillatory try using
    the option Method->Oscillatory in NIntegrate.

Out[2]= 0.459698

In[3]:= Integrate[Sin[x], {x, 0, 1}]

Out[3]= 1 - Cos[1]

In[4]:= %//N

Out[4]= 0.459698

指定できる数値積分法はいろいろあります。（参照：数値積分の実装）