Now Loading...
※ ZIPファイルを展開して「 VCSSL.jar 」をダブルクリックして下さい。 動作には Java が必要です。 このVCSSLコードのライセンスは パブリックドメイン(CC0)です(詳細: ReadMe.txt)。

矩形法(長方形近似)による数値積分

積分の値を数値的に求めます。長方形の短冊で近似した微小領域を足しあげる、最も単純な方法を使用します。

概要

このVCSSLプログラムでは、積分の値を数値的に求めます。 アルゴリズムとしては、長方形の短冊(たんざく)で近似した微小領域の面積を足しあげる、 もっとも単純な方法を使用します。

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

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

[ 関連記事 ]

使用方法

ダウンロードと起動

上の画面の「 ダウンロードして実行 」ボタンを押して ZIP ファイルをダウンロードし、 右クリックして「すべて展開」してください。 展開した中にある「 VCSSL.jar(JARファイル) 」をダブルクリックすると、プログラムが起動します。

起動後

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

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

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

題材解説

手計算(解析計算)における「 積分 」

積分 」と聞くと、面倒だなあとか、苦手だなあとか、マイナスなイメージを持つ人も多いかもしれません。 街を行く人々に「 あなたは積分が好きですか? それとも嫌いですか? 」とアンケートを取ったとしたら、 恐らく後者の方がかなり多いのではないでしょうか。

適当に複雑そうな積分の式

その理由としては恐らく、手計算では:

  • 積分を微分から逆算して解くのに、訓練と慣れが必要な事
  • いくつかの定石のようなテクニックを適材適所で使いこなす難しさ
  • 途中計算が長くなりがちでミスしやすい事
などなどが挙げられるでしょう。この記事のテーマも積分なので、内心「うへぇ」と思われるかもしれません。 しかし、今回扱う「 数値積分 」では、上記のような事は必要ないので、 上のような理由で積分が大嫌いだ、という方でも安心して使えます。

ところで、手で解くと上のように何かと面倒な積分ですが、まあそれでも、 最近は数式を手計算のように自動計算してくれるソフトとかが普及してきて、 恐ろしく込み入った積分の答えを求めてくれたりもします。 しかし、そもそも 手計算で(解析的に厳密に)解く事ができない積分 も、たくさんあります。 というより、解けるように作られた学校のテストとかではなく、実用上で積分を計算する必要性に出くわした場合、 手計算では解けない事のほうが多いかもしれません。

そこで、今回扱う「 数値積分 」の出番、というわけです。

積分の値を、数式ではなく数値として求める「数値積分」

数値積分とは、積分の答えを数式ではなく数値として求める計算の事です。 と言っても、これだけではピンとこないと思われますので、例を出しましょう:

\[ \int^{1}_{0} \cos x dx \]

この定積分の値を求めたいとします。従来の手計算(解析計算)では、積分が微分の逆演算である事と、 sin を微分すれば cos になる事を我々は知っているので: \[ \int^{1}_{0} \cos x dx \] \[ = [ \sin x ]^{1}_{0} \] \[ = \sin 1 - \sin 0 \] \[ = \sin 1 \] \[ = 0.8414709848... \] と計算するでしょう。最後の所だけ関数電卓で数値を求めますね。これが、高校数学でも習った、 普通に手で積分を解く流れです。

で、いよいよ本題ですが、数値積分というのは、上の途中式の過程をすっ飛ばして、 つまり「 sin を微分すれば cos になる 」とかを一切知らなくても、 いきなり最後の 0.8414709848... という値を求められる方法です。 それも、けっこう簡単に。プログラム的には for 文のループ 1 個で済んでしまいます。

「 微分の逆 」ではなく「 面積を求める 」、数値積分のイメージ

それではまず、数値積分で積分値を求めるイメージについて解説します。

積分というと、手計算では「微分の逆」というイメージが強いですが、 数値計算を行う上でそのイメージは別に必要ありません。 むしろ数値積分で必要なのは「 積分する = 面積を求める 」というイメージです。

区間 a から b において、関数 y = f(x) と x 軸で囲まれた領域の図
区間 a から b において、関数 y = f(x) と x 軸で囲まれた領域の図
この領域の面積は、f(x) を a から b まで積分した値と一致します。

上の図の通り、f (x) を区間 a から b まで積分するというのは、 f (x) と X 軸との間の面積を a から b まで求めるという事と、同じ事です。 実際に高校数学の区分求積法を用いて、適当な関数でこの面積を計算してみると、 微分の逆演算として求めた定積分の値と、確かに一致する事が分かります。

領域を細長い長方形の短冊(たんざく)で刻んで足しあげ、面積を求める「短冊法(長方形近似)」

では、実際にプログラミングで f (x) と X 軸との間に囲まれた面積をどうやって求めるか、 という処理を考えていきましょう。 今回扱うのは、精度を出す効率はあまり良くないけれども、最も単純で、処理内容の意味も分かりやすい、 「 短冊方(長方形近似) 」による数値積分法です。 これは、領域を細長い長方形の短冊(たんざく)で刻んで、 それを足しあげる事で、面積を近似的に求めるという発想のアルゴリズムです。 短冊で刻む数(分割数)を n とします。

領域を、
	n 本の細長い長方形の短冊で刻んだ図
領域を、n 本の細長い長方形の短冊で刻んだ図
この長方形の面積をすべて足しあげれば、 関数 f(x) と x 軸の間で囲まれた面積、つまり積分の値を(近似的に)求められます。

まあ、大体近い面積は求まるでしょうが、上の辺とか、 本来は曲線であるところをカクカクに近似してるので、あからさまにはみ出てたり欠けてたり、 明らかに誤差が大きそうですね。 「 だったら長方形の刻み(分割数)をめっちゃ細かくすれば、誤差もそれなりに抑えられるでしょ 」という、 やや強引な感じのアルゴリズムでもあります。

でも、方法としては最も基礎的な感じですし(というより区分求積そのもの)、 分かりやすいのは分かりやすいですよね。 なんとなく、積分というもののルーツ自体も、こんな発想から始まってるのかな、という気もします。 単純明快で、原理に近そう、というか。

さて、上の図で積分区間 a から b までを n 本の短冊で刻み、長方形の幅を Δx 、 そして i 番目の長方形の左下のX値を xi としましょう。つまり: \[ \Delta x = (b - a) / n \] \[ x_i = a + i \Delta x \] すると、全ての短冊を足しあげて近似した、領域全体の面積は: \[ 領域全体の面積(=積分値)= \sum_{i=0}^{n-1} f (x_i) \Delta x \] となります。「 積分 = 面積 」から、この面積の値が、 求まった積分値(近似値)です。つまり数値積分の結果です。

と、以上の通りですが、つまるところ数値積分で積分の値を求めるのにどういう計算が必要かというと、 n 等分した点での f (x) の値に、刻み幅 Δx をかけて足しあげるだけ。 例えば cos x を積分したいんだったら cos x に Δx をかけて足しあげればいいわけで、非常に簡単です。 冒頭でも述べた通り、「 sin を微分すると cos になるから、逆に… 」とかそういう事を一切考えなくても、 積分値が求まってしまうわけです。

コード解説

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

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

※ C/C++ の結果をグラフ化したい場合は…
VCSSLのコードでは実行後に自動でグラフプロットまでが行われますが、 C/C++ のコードでは標準出力に積分値が出力されるのみです。 C/C++ の計算結果をグラフ化して見たい場合は、実行時にリダイレクトなどで標準出力内容をファイルに記録し、 適当なグラフソフト(例えば リニアングラフ2D など) でそのファイルを開いて、プロットしてください。

コード全体

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

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

double 型について
上のコードで double 型を使っていますが、 VCSSL では double 型は float 型と同じものとして扱われ、 現行の処理系では共に倍精度なので、別に float 型を使っても構いません。 今回は、C/C++ のコードと統一するために double 型を用いました。

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

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

以下では、各部の処理について、もう少し詳しく見てみましょう。 なお、上の通り全体的に C / C++ / VCSSL で似通ったコードであるため、 これ以降の詳細解説では、VCSSL のコードのみを抜き出して掲載していきます。

計算条件となる、被積分関数やパラメータ定数などの下準備

まず先頭付近(import や include の部分)では、使うライブラリ・ヘッダを読み込んでいます。 続いて、被積分関数の f(x) や積分区間、短冊での刻み数(分割数)など、 計算条件となる関数・パラメータ定数をまとめて宣言しています。

今回は冒頭の解説そのままに \[ \int^{1}_{0} \cos x dx \] の値を求める事にするので、 被積分関数 f(x) は \[ f(x) = \cos x \] とし、積分区間もそのまま \[ a = 0 \] \[ b = 1 \] とします:

このコードはVCSSLのものですが、C++も同じ内容で、C言語では定数を const する代わりに define しているだけです。

ここまでは、いわば下準備です。

数値積分の計算部

さて、いよいよ中核、「 短冊方(長方形近似)」による数値積分の計算部分です。

ここでは、題材解説のパートで述べた通り、素直に \[ 領域全体の面積 ( =積分値 ) = \sum_{i=0}^{n-1} f (x_i) \Delta x \] の値を計算しています。高さ f ( xi ) で幅 Δx の細長い長方形を足しあげて、面積を求めるイメージでしたね。 ここで Δx の値は、a から b までを n 等分するので \[ \Delta x = (b - a) / n \] です。xi は i 番目の長方形の左下のX値で、これも長方形の幅が Δx なのでそのまま素直に \[ x_i = a + i \Delta x \] です。

実際にこの計算を行っている部分のコードは:

の通りです。変数 value に領域全体の面積( = 積分値)を求めています。for 文の中で、 長方形の短冊の面積を、1本ずつ value に足しこんでいく流れです。

なお、積分結果の出力(表示)については、 「 a から b まで積分した値を求める 」という目的からすれば最後の行だけで十分なのですが、 それだと絵的に地味ですし、 「 短冊を足しあげながら積分していく過程をグラフにしたほうが面白い 」という事で、 計算途中のループ内でも、その時点までの積分値を毎回出力しています。

実際にこのプログラムを実行すると出力される内容は以下の通りです:

0.0    0.0
0.1    0.1
0.2    0.1995004165278026

0.8    0.7319228590332298
0.9    0.8015935299679464
1.0    0.8637545267950129

a から b まで積分していっている最中の x 値が左の列、 その地点までの積分値が右の列に出力されています。 1行ごとに、積分値に長方形の短冊1個の面積が足されています。処理のイメージがつかめて面白いのではないでしょうか。

最後の行の積分値 0.8637545267950129 が、つまるところ数値積分で求まった \[ \int^{1}_{0} \cos(x) dx \] の値で、厳密な値は sin 1 なので 0.8414709848... です。ちょっと食い違ってますね。 これが数値積分の誤差です。今は n = 10、 つまり領域をたった 10 個の短冊に刻んで足しあげたので、まあこんなもんでしょう。 n を大きくすると誤差も抑えられていきます(後述)。

さて、このコンソールに出力された内容をグラフにプロットする処理は、プログラムの最後:

の部分です。ここは今回あまり詳しくは触れませんが、load 関数でコンソールの内容を文字列として取得し、 それを setGraph2DData 関数でグラフにプロットさせています。 プロットされるグラフは以下の通りです。

計算過程の積分値の変化をプロットしたグラフ( n = 10 )
計算過程の積分値の変化をプロットしたグラフ
コンソールに出力された、計算過程の積分値の内容を、そのままグラフにプロットしたものです。

このグラフは、積分途中の x の地点までの積分値、つまり上端を変えながら積分値をプロットしているのと同じ事なので、 従って理論上の厳密な関数は \[ F(x) = \int^{x}_{0} \cos(x') dx' = \sin x\] になるべきです。確かになんとなく sin x の原点付近っぽい形ですね。多少誤差で歪んでいるはずですが。

ところで、現在のVCSSL処理系で採用している2Dグラフソフトは「 リニアングラフ2D 」で、 データに数式を重ねてプロットする機能が付いています。なので、実際に sin x のグラフを重ねて見てみましょう。 まず、グラフ画面のメニューバーから 「 Tool 」 > 「 Math 」と選択します。 続いて表示されるウィンドウの「 f(<x>)」の入力欄に「 sin(<x>) 」と入力し、 「 PLOT 」ボタンを押します。すると:

数値積分で求めた積分値( n = 10 )と、厳密な積分値との比較
分割数 n = 10 で数値積分により求めた積分値と、厳密な積分値との比較
先ほどの図に、厳密な積分値である sin x (青い線)を重ねてプロットしたものです。2つの曲線の差が、数値積分の誤差です。

青い線が追加で描画されましたね。赤い線が先ほどの数値積分の計算結果、 青い線が厳密値である sin x です。微妙にずれてますね。この両者の差が、数値積分の誤差です。

誤差と n と精度限界

さて、この誤差は基本的に、本来ならば上辺が曲線的である積分領域の面積を、 短冊で刻んでカクカクに近似している事によるものです。

領域を n 本の細長い長方形の短冊で刻むんだ図
領域を n 本の細長い長方形の短冊で刻んだ図
本来ならば上辺が曲線的であるのに、長方形の短冊で刻んでカクカクに近似したたのが、誤差の原因です。

なので、短冊をどんどん細長く、つまり短冊で刻む数 n をどんどん増やしていけば、 上図で曲線から欠けたりはみ出したりしている部分が抑えられていって、 近似している面積も真の面積に近づいていき、 誤差もどんどん抑えられていくはずです。

という事で、n = 10 から刻みを1桁細かくして、n = 100 にして計算してみましょう。 まずその計算結果のグラフに、先ほどと同様に厳密値の sin x を重ねてみると、以下のようになります:

数値積分で求めた積分値( n = 100 )と、厳密な積分値との比較
分割数 n = 100 で数値積分により求めた積分値と、厳密な積分値との比較
n = 10<の場合と比べると、2つの曲線がより接近して重なっている、つまり誤差がより抑えられた事が見て取れます。

一気に重なるようになりましたね。ぱっと見でも誤差が抑えられています。 数値で誤差を見てみましょう。コンソールの出力結果は:

0.0    0.0
0.01    0.01
0.02    0.019999500004166653

1.0    0.8437624610086617

の通りです。最後の行の 0.8437624610086617 が最終的な積分結果で、 厳密な値は sin 1 = 0.8414709848... です。

先ほどの n = 10 の場合は 0.86... で、小数点以下2桁目から違っていたのに対して、 今度の n = 100 では2桁目は合っていて、3 桁目から違っていますね。 つまり刻み数 n を1桁増やす事で、精度が1桁あがりました。 では n をもう1桁増やして、n = 1000 で計算するとどうかと言うと、0.8417007635323798 になります。 今度は3桁目まで厳密な値と一致して、4桁目から違っていますね。 ここまでの流れだと、n を1桁増やすごとに、 誤差が1桁程度ずつ抑えられていっている、 つまり信頼できる精度が1桁程度ずつ上がっていっている事が分かります。

ところで、n を増やすという事は、プログラムの for 文のループ回数をそれだけ増やすという事ですから、 そっくりそのまま計算量(≒計算時間)も増えます。 つまるところ、今回用いた短冊方(長方形近似)による数値積分では、 信頼できる精度を1桁程度増やすために、計算量も1桁つまり10倍増やす必要があるわけです。 これは、色々と種類がある数値積分のアルゴリズムの中でも効率が悪く、高精度な値を出すのに結構な計算量が必要です。

さて、「 たとえ効率が悪くても、いくら時間がかかってもひたすら待てば、 コンピュータが扱える限界まで精度を上げられるのか 」と言うと、 実は今回のように精度の効率が悪いアルゴリズムだと、限界があります。 というのも、コンピュータは、本来は(割り切れなかったりで)無限に続く数値の桁数を、 有限の桁数に丸めて計算しているわけなので、演算のたびに誤差(丸め誤差)が生じています。 演算をする度にこの誤差が生じ、今回のように繰り返し足しこんでいくような計算だと、その誤差も蓄積していきます。 なので、あまりにも計算量を大きくし過ぎると、 この誤差蓄積が大きくなってきて、逆に計算結果の精度が落ちてきたりするわけです。 本末転倒ですね。

※ コンピュータでの数値の扱いについて詳しく知りたい場合は、IEEE 754 の仕様などを参照してください。 IEEE754 - Wikipedia

より高精度なアルゴリズムへ

じゃあどうするのか、というと、もっと精度の効率の良いアルゴリズムを使用します。 つまり、計算量 n を1桁(10倍)増やすだけで、精度が2桁程度上がったり、 もしくは4桁程度一気に向上するようなアルゴリズムを使うわけです。 そうすると、もっとずっと高速に計算できて、誤差蓄積も抑えられるので精度の限界も上がり、 実用性が向上するわけです。

次回は、今回扱った「 短冊方(長方形近似) 」を少し改修するだけで、 精度の効率が大きく上がる、「 台形法 」について解説します。

コードのライセンス

このVCSSLコードは著作権フリー(パブリックドメイン)で公開しています。C言語版のコードも同様です。 そのままでのご利用はもちろん、言語の種類を問わず、改造や流用などもご自由に行ってください。



スポンサーリンク



台形法(台形近似)による数値積分

積分の値を数値的に求めます。長方形近似よりも高精度な方法として、台形で近似した微小領域を足しあげる方法を使用します。
短冊法(長方形近似)による数値積分

積分の値を数値的に求めます。長方形の短冊で近似した微小領域を足しあげる、最も単純な方法を使用します。
小数(浮動小数点数)から分数へ近似的に変換する

小数(浮動小数点数)を、適当な誤差の範囲内で、近い分数に変換するプログラムです。
円周率1万桁の計算(ガウス=ルジャンドル法)

ガウス=ルジャンドル法により、円周率を1万桁まで計算するプログラムです。
試し割り法による素数判定

試し割り法を用いた、素数判定のプログラムです。
スポンサーリンク

インデックス
台形法(台形近似)による数値積分

積分の値を数値的に求めます。長方形近似よりも高精度な方法として、台形で近似した微小領域を足しあげる方法を使用します。
短冊法(長方形近似)による数値積分

積分の値を数値的に求めます。長方形の短冊で近似した微小領域を足しあげる、最も単純な方法を使用します。
小数(浮動小数点数)から分数へ近似的に変換する

小数(浮動小数点数)を、適当な誤差の範囲内で、近い分数に変換するプログラムです。
円周率1万桁の計算(ガウス=ルジャンドル法)

ガウス=ルジャンドル法により、円周率を1万桁まで計算するプログラムです。
試し割り法による素数判定

試し割り法を用いた、素数判定のプログラムです。

お知らせ( Twitter )

スポンサーリンク