矩形法（長方形近似）による数値積分

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

なお、今回のプログラムは、あくまでアルゴリズム解説のためのサンプルコード的なものなので、 リッチな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コマンドが使用できない等のエラーが表示される場合は…

起動後

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

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

被積分関数や積分区間、精度などを変えたい場合は、 プログラム「 IntegralRectangular.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 個で済んでしまいます。

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

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

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

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

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

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

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

でも、方法としては最も基礎的な感じですし（というより区分求積そのもの）、 分かりやすいのは分かりやすいですよね。 積分というもののルーツ自体も、こういう発想が出発点だろうなという感じがしますし、 「そもそも積分とは何をする事なのか」という直感的なイメージを掴むためにも、 こういう単純な数値積分法に触れてみるのはおすすめです。

さて、上の図で積分区間 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++で記述したコード も一緒に掲載します（※）。

コード全体

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

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

なお、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 まで積分した値を求める 」という目的からすれば最後の行だけで十分なのですが、 それだと絵的に地味ですし、 「 短冊を足しあげながら積分していく過程をグラフにしたほうが面白い 」という事で、 計算途中のループ内でも、その時点までの積分値を毎回出力しています。

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

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

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

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

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

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

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

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

誤差と n と精度限界

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

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

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

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

の通りです。最後の行の 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倍増やす必要があるわけです。 これは、色々と種類がある数値積分のアルゴリズムの中でも効率が悪く、高精度な値を出すのに結構な計算量が必要です。

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

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

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

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

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

• 第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 氏の米国およびその他の国における商標または登録商標です。
• その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。