シンプソン法による数値積分

このVCSSLプログラムでは、数値積分シリーズの連載の第3回として、前回に引き続き、積分の値を数値的に求めます。 アルゴリズムとしては、前回からもう一歩進歩して、 微小領域の面積をより高い精度で求めて足しあげる、シンプソン法という方法を使用します。 矩形法（長方形近似）や台形法（台形近似）とのコード上での違いは、今回もたったの一行なのですが、 それで計算精度が台形法からさらに劇的に向上するのがポイントです。

なお、今回のプログラムは、あくまでアルゴリズム解説のためのサンプルコード的なものなので、 リッチなGUIの入力画面などはありません。 被積分関数やパラメータを変えたい場合などは、コード内の記述を直接書き換えてください （ または、上のおすすめリンクにある、GUI版の積分ツールをご利用ください ）。

コードはVCSSLで記述されていますが、 適当に別の言語などに書き換えて使って頂いても構いません（ C言語版のコード と C++版のコード も下で掲載しています）。 このコードは著作権フリー（パブリックドメイン / CC0）なので、自己責任でご自由にご使用ください。

[ 数値積分シリーズの連載記事 ]

• 第1回： 矩形法（長方形近似）による数値積分
• 第2回： 台形法（台形近似）による数値積分（前回）
• 第3回： シンプソン法による数値積分（今回）
• 総集編 ： 積分値を求めるプログラム
• 応用 ： 入力された数式を積分して値とグラフを表示するツール

スポンサーリンク

---

使用方法

ダウンロードと展開（解凍）

まず、PC（スマホは未対応）で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。

その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開（解凍）してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…

プログラムの起動

Windows をご使用の場合

上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください：

もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。

正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。

» うまく起動できずにエラーになってしまう場合は…

Linux 等をご使用の場合

ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください：

» javaコマンドが使用できない等のエラーが表示される場合は…

起動後

起動すると数値積分の処理が実行され、積分された関数のグラフが表示されます。また、コンソールに積分値も表示されます。

被積分関数やパラメータの変更など

被積分関数や積分区間、精度などを変えたい場合は、 プログラム「 IntegralSimpson.vcssl 」をテキストエディタで開いて、 コードを適当に書き換えてください。

題材解説

前回までのおさらい
―― 矩形法と台形法

今回は、前回までの記事の発展的な内容となります。そこで、前回までの内容を軽くおさらいしてみましょう。 なお、土台をしっかりと固めたいという方は、先に前回までの記事をご参照ください。

• 第1回： 矩形法（長方形近似）による数値積分
• 第2回： 台形法（台形近似）による数値積分

第1回では、積分の値を数値的に求める方法として、 X軸と \( y = f(x) \) との間に囲まれた面積を、長方形で細かく区切って求めました。 イメージとしては下図の通りになります。

このようにX軸と \( y = f(x) \) との間に囲まれた面積が、 理論上は関数 \( f(x) \) の積分値と等しいのでしたね。 なので、積分の数式を手で解かなくても、プログラムで数値的にこの面積を求めてやれば、積分の値が求まるわけです。

第1回で扱った矩形法（長方形近似）では、実際に上の細長い長方形の面積を足し上げる事で、全体の面積（ ＝ 積分値 ）を求めました。 でも、そうすると、どうしても「 誤差 」が生じます。 これは上図からも明らかで、領域の上端の線は本来はなめらかな曲線なのに、その領域を長方形を並べたもので強引に近似しているわけですから、 上端がギザギザになって、本来の曲線からはみ出ていたり欠けていたりする箇所が生じています。 これが大きな誤差の元になっていたわけです。

そこで第2回で扱った台形法（台形近似）では、上図の細長い長方形の部分を台形に置き換えて改良する事で、 はみ出たり欠けている部分を減らして、より高精度に全体の面積（ ＝ 積分値 ）を求めました。このイメージは下図の通りです。

上の図を見ただけでも、矩形法では上端の曲線からガタガタに外れていた部分が、台形法では大分ならされているので、誤差が一気に減りそうな事が直感的にわかります。

実際、細長い微小領域の幅 \( \Delta x \) を小さくしていった際に、矩形法では積分値の誤差がそれに比例して小さくなっていったのに対し、 台形法では \( \Delta x \) の2乗に比例して（より急激に）小さくなっていきました。 具体的には、\( \Delta x \) の幅を10分の1にすると、矩形法では誤差も10分の1に減るのに対して、台形法は誤差が一気に100分の1まで減る、といった具合です。 \( \Delta x \) を100分の1にするなら、矩形法なら誤差も100分の1ですが、台形法なら1万分の1になります。この違いは大きいですね。

台形からさらに、誤差をどうやって軽減するか ?
―― 微小領域の上端を曲線で近似

さて、少し長いおさらいでしたが、ここから今回の本題です。 台形法から、もう一歩改良して、精度をさらに向上させてみましょう。

と言っても、大々的に複雑な改修をやるわけではないので、安心してください。最終的に変更するのは、プログラムのコードだと、たったの1行です。 具体的には、細長い微小領域の面積を計算していた部分です。 そう、結局は今回も前回同様、「微小領域の面積をいかに精密に求めるか」という話に尽きるわけですね。 それ以外の部分の考え方や処理の流れは、これまでとぜんぶ同じです。

それではまず、台形法の誤差の原因を見つめ直してみましょう。細長く刻んだ微小領域（ ＝ いまは台形で近似 ）を1個抜き出して注目します。 わかりやすくするために、f(x) のカーブの曲がり具合は先ほどの図よりも誇張しています。

微小領域の上端の線を見ると、矩形法よりはよっぽどマシですが、それでもやはり曲線であるはずの \( f(x) \) を直線で近似しているため、 本来の微小領域の面積と比べて、台形の面積はわずかに食い違っています（図中のオレンジ色の部分）。この食い違っている部分は、\( f(x) \) のカーブ具合に応じて、場所によってはみ出ていたり欠けていたりするでしょう。 この過不足が台形法での誤差の根本原因です。 この微小領域の面積を足して、積分区間全体での面積を求め、それを積分値とするわけですが、その際にこの小さな過不足が積み重なって、最終的に積分値の誤差を生み出しているのです。

では、どうするか？ 「 微小領域の上端で、本来は曲線である \( f(x) \) を直線で近似しているから誤差が生じるなら、そこを曲線で近似してやれば精度が上がるのでは？ 」 これが、今回の土台となる考え方です。

直感的にも、より欠けたりはみ出したりする部分が減り、誤差が減って精度が向上しそうですね。

微小領域の上端をどういう曲線で近似するか？
―― 単純で手でも積分可能な二次関数を使用

ここから先は、微小領域の上端で、\( f(x) \) をどういう曲線で近似するか？ という話になるわけですが、ここが今回の肝になるポイントです。

一番誤差が無く完璧なのは \( f(x) \) そのものを使う事ですが、そうすると微小領域の面積を求めるためにまず \( f(x) \) を手計算で積分する必要が出てくるので、意味がなくなってしまいます。 そもそもそれが難しい（ or 不可能な ）場合のために、プログラムを書いてまでコンピューターに積分をやってもらおうとしているのでしたね。

ならどうするかと言うと、二次関数を使います。二次関数は、下図にも示す通りシンプルな曲線で、そして高校などでも習うように、手でも厳密に積分できますね。 なので、二次関数で \( f(x) \) をいい感じに近似できれば、細長い微小領域の面積は手でも解ける ＝ 式として求まるので、 その式で、台形法で使っていた台形の面積の式を置き替えればいいわけです。

さて、二次関数は中学校でも習う通り、係数などのパラメータによって丸みや位置が変わります（上図参照）。 「 y = 2 x2 + 3 x - 4 のグラフを描け 」みたいなのがテストに出て、四苦八苦しながら位置や丸みを求めて手描きした思い出がある方も、結構多いのではないでしょうか。 ここでは \( f(x) \) の近似が目的なので、微小領域の範囲内で \( f(x) \) のカーブになるべくピタッとフィットするように、こちら側で二次関数の係数を設定してやりましょう。

多項式を関数にフィットさせるための係数の選び方は、一般には色々あります。ただ、あまり計算が複雑になるのも避けたいところです。 そこでここでは結構単純に、細長い微小領域の「両端と真ん中」の3つの地点 ※ において、関数 \( f(x) \) と同じ点を通るように、二次関数の係数を設定します。 両端の2点だけなら、そこを通る二次関数はいくらでも存在しますが、その中には上に膨らんでいたり下に凹んでいたりと様々なものがあるので、 「真ん中の点も通る」というルールを加える事で、\( f(x) \) のカーブ具合にそこそこフィットしてもらおう、というわけですね。

このように3点を通る二次関数は1通りに絞られて、ここでは割愛しますが、手計算で求められます。 そうすると、あとはその二次関数を積分すれば、細長い微小領域の面積が求まります。 これはシンプソンの公式として知られていて、以下のようにかなり簡単な形で書けます：
\[ i番目の微小領域の面積 = \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \]
ここで \( x_i \) は微小領域の左端、\( x_i+\Delta x \) は右端、\( x_i + \Delta x/2 \) は真ん中の x 値であり、\( \Delta x \) は微小領域の幅（ \( = x_{i+1} - x_{i} \) ）です。

あとはこれを、全ての微小領域について足しあげれば、領域全体の面積が求まります：
\[ 積分値 = \sum_{i=0}^{N-1} \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \]

式の上では、もう少し重複部分について整理した形にできますが、コードとの対応では上の形が一番ストレートでわかりやすいので、このままにしておきましょう。 要は、矩形法や台形法で、 足し上げる中身が長方形や台形の面積だったのが、シンプソンの公式を使った二次関数近似の面積に置き換わっただけですね。

---

スポンサーリンク

---

コード解説

それでは、実際にプログラミングで、シンプソン法により数値積分するコードを見ていきましょう。 このプログラムのコードはVCSSLで記述されています。 VCSSLのコードは大体C言語っぽい感覚で読めると思います。

また、今回は C言語で記述したコード と C++で記述したコード も一緒に掲載します（※）。

コード全体

それでは、コード全体を見てみます。まずは、VCSSLで記述したコードです。50行に満たない程度の短いコードです。

以上です。流れとしては、先頭あたりで被積分関数やパラメータ変数を設定し、 真ん中あたりで for 文のループを回して数値積分の計算を行い、 最後に積分の過程をグラフにプロットするという具合です。

なお、C言語で記述したコードは以下の通りです。

C++で記述したコードは以下の通りです。

以上です。それでは、もう少し具体的に見てみましょう。

前回との違いは 1 行だけ !

コード各部の細かい解説については、矩形法や台形法からほとんど変わっていないので、 矩形法の数値積分のコード解説をご参照ください。 矩形法や台形法のコードと違っている部分は、以下のたった一か所だけです。

[ ↓ 矩形法のコード ]

[ ↓ 台形法のコード ]

[ ↓ 今回（シンプソン法）のコード ]

これはつまるところ、 長方形や台形の面積を求めていた箇所を、先ほど説明した、シンプソンの公式を使って二次関数近似で面積を求める式：
\[ i番目の微小領域の面積 = \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \] に変更しただけですね。

さて誤差は ?

それでは、このたった1行だけの変更で、どれだけ誤差が減るのか見てみましょう。 このプログラムを実行すると、まず積分過程がグラフで表示されます。 今は被積分関数を \( f(x) = \cos x \) と置いているので、 このグラフはそれを積分した \( \sin x \) の形をしているはずです。

上の図で、赤色の線（重なっていて見づらいですが…）が、今回の計算結果です。 青色の線は、厳密な値と比べるために、\( \sin x \) のグラフを重ね描きしたものです。 前回扱った台形法の時点でも既に両者はかなり重なっていましたが、今回もほぼ重なってしまって区別できません※。 この結果は刻み幅 N = 10、つまり微小区間をかなり粗く刻んでいるのにもかかわらず、です。

グラフをものすご〜く拡大すれば一応は厳密値とのズレが見られますが、この水準になると、もう数値で直接ズレを見た方がわかりやすいでしょう。 という事でコンソールを見ると、以下のような内容が出力されているはずです：

このように、a から b まで積分していっている最中の x 値が左の列、 その地点までの積分値が右の列に出力されています（上で示したグラフでは、このデータを読ませてプロットしています）。 1行ごとに、積分値に微小領域1個の面積が足されています。

最後の行の積分値 0.8414710140343371 が、最終的に求まった \[ \int^{1}_{0} \cos(x) dx \] の値で、厳密な値は sin 1 なので 0.84147098480789650... です。小数点以下 6 桁目から食い違っていますね。 この食い違いが誤差です。

前回扱った台形法では、同じ刻み数 N=10 の条件下で、小数点以下 3 桁目から違っていました。 という事は、シンプソン法で精度が一気に上がった事がわかりますね。 ちなみに、矩形法では小数点以下 2 桁目から違っていました。

さて、微小区間の刻み幅 \( \Delta x \) を細かくするほど、つまり刻み数 N を大きくするほど※、誤差は小さく抑えられるはずです （※ 詳細は矩形法参照ですが、aを積分区間の下端、bを積分区間の上端として、それをN個の微小区間に刻むため、\( \Delta x = (b-a)/N \) です ）。

実際にプログラムにおいて N を増やしていき、先ほどの最終行の値がどうなっていくか見てみましょう：

わかりやすいように、厳密な値（ \(\sin 1\) ）から食い違っている部分を赤色で塗りました。

N = 1 から 1000 までの結果に注目してみると、刻み数Nを1桁増やすごとに、合っている桁数がおおまかに4桁くらいずつ増えている事がわかります（※ラスト2桁の誤差については後述）。 この結果から、シンプソン法は「 4次精度 」の数値積分アルゴリズムである事がわかります。

ところで、前回扱った台形法は2次精度で、Nを1桁増やすごとに、合っている桁数が大体2桁くらいずつ増えるのでしたね。 また、最初に扱った矩形法は1次精度だったので、1桁ずつしか増えませんでした。 一方で、4次精度を誇るシンプソン法を使えば、Nをちょっと増やすだけで一気に「4段飛ばし」で精度が増えていくわけです。 これらの数値積分法のプログラムのコードは、それぞれたった1行違うだけなのに、その実力差は歴然ですね。 このあたりから、「 計算アルゴリズムっておもしろいなー 」と感じてもらえると幸いです。 一見すると少し難解な計算アルゴリズムでも、自分の手でコードを書いて実装して、高精度な値がバッチリ出た瞬間は、不思議さと楽しさが入り混じって最高にテンションが上がりますよ！

[ 数値積分シリーズの連載記事 ]

• 第1回： 矩形法（長方形近似）による数値積分
• 第2回： 台形法（台形近似）による数値積分（前回）
• 第3回： シンプソン法による数値積分（今回）

ライセンス

このVCSSL/Vnanoコード（ 拡張子が「.vcssl」や「.vnano」のファイル ）は実質的な著作権フリー（パブリックドメイン) である CC0 の状態で公開しています※。 記事中にC言語/C++/Java言語などでのサンプルコードが掲載されいてる場合は、それらについても同様です。 そのままでのご利用はもちろん、改造や流用などもご自由に行ってください。

※ ただし、このコードの配布フォルダ内には、ダウンロード後すぐに実行できるように、 VCSSLの実行環境も同梱されており、そのライセンス文書は「 License 」フォルダ内に同梱されています （要約すると、商用・非商用問わず自由に使用できますが、使用の結果に対して開発元は一切の責任を負いません、といった具合の内容です）。 配布フォルダ内の各構成物の一覧やライセンスについては「 ReadMe_使用方法_必ずお読みください.txt 」をご参照ください。

※ Vnano の実行環境については、別途スクリプトエンジンのソースコードも一般公開しており、 何らかのソフトウェア内に組み込んでご利用いただく事も可能です。詳細はこちらをご参照ください。

この記事中の商標などについて

• OracleとJavaは、Oracle Corporation 及びその子会社、関連会社の米国及びその他の国における登録商標です。文中の社名、商品名等は各社の商標または登録商標である場合があります。
• Windows は、米国 Microsoft Corporation の米国およびその他の国における登録商標です。この記事は独立著作物であり、Microsoft Corporation と関連のある、もしくはスポンサーを受けるものではありません。
• Linux は、Linus Torvalds 氏の米国およびその他の国における商標または登録商標です。
• その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。