このエントリーをはてなブックマークに追加
2014/11/01

長方形近似による数値積分

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

RINEARN CODE - 長方形近似による数値積分
スクリーンショット
※ ダウンロードしたZIPファイルを展開し、中のVCSSL.jarをダブルクリックして下さい。
※ 各種デスクトップOSで動作しますが、動作には Java が必要です。

概要

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

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

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

使用方法

このプログラムは、Webページ上でも、ダウンロードしてでも使用できます。

Webページ上で起動

上の画面の実行ボタンを押すと、実行用のWebページに移動し、プログラムが起動します。

ダウンロードして起動

上の画面のダウンロードボタンを押して ZIP ファイルをダウンロードし、 解凍した中にある「 VCSSL.jar 」をダブルクリックすると、プログラムが起動します。 なお、ZIPファイルの解凍は、右クリックして「 すべて展開 」などで行えます。

起動後

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

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

被積分関数や積分区間、精度などを変えたい場合は、 プログラム「 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 になるから、逆に… 」とかそういう事を一切考えなくても、積分値が求まってしまうわけです。

コード解説

それでは、実際にプログラミングで数値積分するコードを見ていきましょう。 このプログラムのコードは、Cライクな簡易プログラミング言語 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言語版のコードも同様です。 そのままでのご利用はもちろん、言語の種類を問わず、改造や流用などもご自由に行ってください。



スポンサーリンク


このエントリーをはてなブックマークに追加
新着プログラム

台形近似による数値積分
2015年06月25日

長方形近似による数値積分
2014年11月01日

小数(浮動小数点数)から分数へ近似的に変換する
2014年08月10日

ユーザーが入力した数式を2次元グラフにプロットする
2013年11月30日

配列を3次元グラフにプロットする
2013年11月28日

配列を2次元グラフにプロットする
2013年11月28日

ファイルを3次元グラフにプロットする
2013年11月27日

ファイルを2次元グラフにプロットする
2013年11月26日

ウィンドウの生成
2013年08月28日

ボタンが並ぶパネルの配置(ButtonPanel)
2013年08月28日

ボタンの配置
2013年08月28日

2DCGと3DCGの合成
2013年05月04日

凹レンズを通過する波のシミュレーション
2013年03月16日

凸レンズを通過する波のシミュレーション
2013年03月15日

乱雑な密度分布における波のシミュレーション
2013年03月12日

ローレンツアトラクタ(ファイル出力版)
2013年02月28日

波の屈折のシミュレーション
2012年12月05日

力学アルゴリズムによる波のシミュレーション(面上の波)
2012年11月21日

手動で波を発生させるシミュレーション
2012年11月19日

力学アルゴリズムによる波のシミュレーション(線上の波)
2012年11月18日

カラーコードとRGBの相互変換と色表示
2012年11月15日

頂点配列によるモデルの変形アニメーション
2012年11月14日

頂点配列によるモデルの作成(四角形格子メッシュ形式)
2012年11月11日

円周率1万桁の計算(ガウス=ルジャンドル法)
2012年08月22日

2DCG用フレームワークの使用サンプル(アニメーション版)
2012年08月09日
  スポンサー リンク

  おすすめ / 人気のコーナー
フリーソフト
RINEARN では、インストール不要の各種解析ソフトを無償公開しています。Windows 8.1 で動作確認済み !
ピックアップ
【VCSSL実践講座】 3Dグラフ描画ツールを作ろう!
C言語系の簡易プログラミング言語「 VCSSL 」で、色々なものを作っていく連載記事コーナーです。現在は3Dグラフ描画ツールを開発中 !
2013年08月29日
【RINEARN CODE】 凸レンズを通過する波のシミュレーション
凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。まっすぐ入射した波がレンズ通過時に屈折し、焦点へ収束する様子がアニメーションで楽しめます。
2013年03月15日
  新着情報
RSSフィード
RINEARN の新着・更新情報をRSSフィードで配信しています !
新着リスト
台形近似による数値積分

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

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

小数(浮動小数点数)を、適当な誤差の範囲内で、近い分数に変換するプログラムです。
2014年08月10日
ベクターソフトニュース様でリニアングラフ3Dをご紹介頂きました!

オンラインソフトウェア流通サイトVector様の「ベクターソフトニュース」コーナーにて、リニアングラフ3Dを取り上げて頂きました。
2014年07月31日
ユーザーが入力した数式を2次元グラフにプロットする

実行時にユーザーが入力した数式の値を、2次元グラフにプロットするサンプルプログラムです。
2013年11月30日
配列を3次元グラフにプロットする

座標値配列の内容を、3次元グラフにプロットするサンプルプログラムです。
2013年11月28日
配列を2次元グラフにプロットする

座標値配列の内容を、2次元グラフにプロットするサンプルプログラムです。
2013年11月28日
ファイルを3次元グラフにプロットする

座標値ファイルの内容を、3次元グラフにプロットするサンプルプログラムです。
2013年11月27日
ファイルを2次元グラフにプロットする

座標値ファイルの内容を、2次元グラフにプロットするサンプルプログラムです。
2013年11月26日
2013年の不具合修正情報

2013年1月〜12月に実施された不具合修正に関する情報を掲載しています。
2013年11月01日
3Dグラフ描画ツールを作ろう!〜第5回縮尺を自動調整する

VCSSL実践講座、3Dグラフ開発編。第5回では、最大最小値に応じて、縮尺を自動調整する処理を実装します。
2013年09月02日
3Dグラフ描画ツールを作ろう!〜第5回縮尺を自動調整する

VCSSL実践講座、3Dグラフ開発編。第5回では、最大最小値に応じて、縮尺を自動調整する処理を実装します。
2013年09月02日
  Twitter