2008年5月3日土曜日

シンプソンの公式の誤差

定積分の数値計算に用いるシンプソンの公式は良く知られているが、その誤差については、あまり知られていない。結論だけ記載しておきます。(出典:微分積分学 吉田洋一著:培風館)

m*h^5*G/90=(b-a)^5*G/(2880*m^4)
G:d^4*f(x)/(dx^4)の最大値
m:閉区間[a,b]の分割数/2
h=(b-a)/(2*m)

この式を眺めただけでどの程度分割すればよいかが推測できます。別トピックで示したExcelVBAのプログラムの分割数は余分に分割していることが推測されます。

下記は菅沼研究室のURL(リンクに記載済み)から引用したC++のプログラム例
http://www।sist.ac.jp/~suganuma/kougi/other_lecture/SE/num/Simpson/C++/Simpson.txt より

/********************************/
/* シンプソン則 */
/* coded by Y.Suganuma */
/********************************/
#include
#include

double snx(double);
double simpson(double, double, int, double(*)(double));

main()
{
double y, pi2 = 2.0 * atan(1.0);

y = simpson(0.0, pi2, 100, snx);

printf("result %f\n", y);

return 0;
}

/********************/
/* 関数値の計算 */
/********************/
double snx(double x)
{
return sin(x);
}

/***********************************************************/
/* 数値積分(シンプソン則) */
/* x1 : 下限 */
/* x2 : 上限 */
/* n : 分割数 */
/* fun : 関数名(f) */
/* return : 積分値 */
/***********************************************************/
double simpson(double x1, double x2, int n, double(*fun)(double))
{
double s1 = 0.0;
double s2 = 0.0;
double s, x, h, h2;

h = (x2 - x1) / n;
h2 = 2.0 * h;

for (x = x1+h; x < x =" x1+h2;" s =" h">

0 件のコメント: