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 件のコメント:
コメントを投稿