Simpson則

寄稿:東條遼平

定積分の近似

Simpson則(シンプソン則)を用いて定積分を求めてみます.

simpson1.png

このようなf(x)を積分区間[a,b]で積分するとします. このとき十分大きな数nを考え,積分区間をn等分します. PCでは無限を考えることができないために有限個にしか分割できませんが, nを大きくすることで誤差を小さくすることができます.

ここでn等分された区間のうちの1つを拡大して見てみます.

simpson2.png

分割された内の1つの区間[x0,x2]を考え,区間の幅をhとします. [x0,x2]の中間値をx1とし,f(x0),f(x1),f(x2)をそれぞれ y0,y1,y2とします.

よく積分の説明で,矩形を考えて幅をゼロに近付けるというのがありますが, それができるのは無限という概念が存在するからで,有限の値しか 扱えないPCで矩形を考えては大きなすき間が空いてしまい, 誤差が大きくなる,または分割数をかなり大きくしなければならなくなります. そのためこの(x0,y0),(x1,y1),(x2,y2)の三点を通る2次曲線で近似します. そうすることですき間を極めて小さくできると考えられます.

img87.png

とおいてx1を中心とした二次式を考えます.こうすると計算が楽になります.

img88.png
img89.png

ここで,

img90.png

より,

img91.png

となり,

img92.png

これより[x0,x2]の範囲で積分すると,

img93.png
img94.png
img95.png

ここで積分区間をn分割したときのi番目の区間を積分した値は,

img96.png
img97.png

となり,積分区間[a,b]で積分すると,

img98.png

ソースコード

以上長くなりましたが,実装するのは簡単です.

double Simpson(double (*f)(double x), double x0, double x1, int num)
{
  double h;
  double sum = 0.0, x;
  int i;

  h = (x1 - x0);
  if(num) h /= num;

  x = x0;

  for(i=0; i<num; i++){
    sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h);
    x += h;
  }

  return sum*h/6;
}

使いまわしがきくように関数ポインタを使用しました. 関数ポインタとは関数をアドレスで渡して使う方法で,動的に使用する関数を変更することが出来る という便利な方法です.main関数より,

double f1(double x)
{
  return x;
}

double f2(double x)
{
  return 4 / (1 + x * x);
}

という2つの関数を順に積分しています.関数名だけを書けば,それが関数のアドレスとなりますので,

Simpson(f1, 2,  5, 1000);

のようにすると,

for(i=0; i<num; i++){
  sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h);
  x += h;
}

の部分は

for(i=0; i<num; i++){
  sum += f1(x) + 4.*f1(x + h/2) + f1(x + h);
  x += h;
}

と同じ処理をします.

Valid XHTML 1.1! home > コンピュータ > プログラミング >
リロード   新規 編集 凍結 差分 添付 複製 改名   トップ 一覧 検索 最終更新 バックアップ   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.5_1 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-jiPowered by PHP 5.2.17HTML convert time to 0.359 sec.