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

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

積分の値を数値的に求めます。長方形近似の発展として、台形で近似した微小領域を足しあげる方法を使用します。

概要

このVCSSLプログラムでは、 前回に引き続き、積分の値を数値的に求めます。 アルゴリズムとしては、前回からもう少し進歩して、 台形の微小領域の面積を足しあげる方法を使用します。 コード上ではたった一行の改良でしかないのですが、それで計算精度が劇的に向上するのがポイントです。

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

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

[ 関連記事 ]

使用方法

ダウンロードと起動

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

起動後

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

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

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

題材解説

前回のおさらい

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

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

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

このようにX軸と \( y = f(x) \) との間に囲まれた面積が、 理論上は関数 \( f(x) \) の積分値と等しいのでしたね。

なので、積分の数式を手で解かなくても、プログラムで数値的にこの面積を求めてやれば、積分の値が求まるわけです。 でも、そうすると、どうしても「 誤差 」が生じます。 これは上図からも明らかで、本来は曲面である領域を、長方形を並べたもので強引に近似しているわけですから、はみ出ていたり欠けていたりする箇所が生じています。 これが誤差の元になるわけです。

誤差をどうやって軽減するか ?
―― 長方形から台形に

さて、ここからは今回の内容です。

上で挙げた図を見ながら、 誤差を手っ取り早く軽減するにはどうすればよいか ? と考えると、 まず長方形をもっと細かく刻んでやる事を思いつきます。 これを実際にやってみると、確かに長方形の刻み幅の小ささに正比例して、 誤差が減っていく事が分かります(前回の記事参照)。 しかし、ある程度以上に刻み幅を細かくしても、計算回数が増えすぎて別の誤差が蓄積してしまい、 限界があります。計算時間もかかります。

そういうわけで、同じ刻み幅でも、 もっと誤差を小さくできないか ? という話が重要になるわけです。 そこで今回は、X軸と \( y = f(x) \) との間に囲まれた面積を、 長方形ではなく台形で区切って近似してみましょう(下図)。

領域を n 本の細長い台形で刻むんだ図
領域を n 本の細長い台形で刻んだ図
このように、長方形の代わりに台形の面積をすべて足しあげれば、関数 f(x) と x 軸の間で囲まれた面積(つまり積分値)を、より高精度に求められるはずです。

この図を、先に挙げた長方形で刻んだ図と比べてみると、 はみ出たり欠けたりしている部分が少なく、ぱっと見でも明らかに誤差が小さそうですね。 実際、結果を先に言ってしまうと、この「台形法(台形近似)」 は、前回の短冊法(長方形近似)よりも劇的に誤差を軽減できます。

さて、具体的に面積を求める計算についてですが、 まず前回同様に積分区間 a から b までを n 本の台形で刻み、 台形の幅を Δx 、そして i 番目の台形の左下のX値を xi としましょう。つまり:
\[ \Delta x = (b - a) / n \] \[ x_i = a + i \Delta x \]
台形の面積は「( 下底 + 上底 )× 高さ(=幅) ÷ 2 」なので、 i 番目の台形の面積は:
\[ i 番目の台形の面積 = \frac{ ( f(x_i) + f(x_{i+1}) ) }{ 2 } \Delta x \] \[ = \frac{ ( f(x_i) + f(x_i + \Delta x) ) }{ 2 } \Delta x \]
となります。あとはこれを、全ての台形について足しあげれば、領域全体の面積が求まります:
\[ 領域全体の面積(=積分値)= \sum_{i=0}^{n-1} \frac{ ( f(x_i) + f(x_i + \Delta x) ) }{ 2 } \Delta x \]
前回からほんの少し変わっただけですね。実際、プログラムに加える変更もわずかで済みます。

コード解説

それでは、実際にプログラミングで、台形法(台形近似)により数値積分するコードを見ていきましょう。 このプログラムのコードは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++で記述したコードは以下の通りです。

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

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

コード各部の細かい解説については、前回とほとんど同じなので、 前回のコード解説をご参照ください。 前回のコードと違っている部分は、実は以下の一か所だけです。

[ ↓ 前回のコード ]

[ ↓ 今回のコード ]

これはつまるところ、 長方形の面積を求めていた箇所を、先ほどの台形の面積を求める式:
\[ i 番目の台形の面積 = \frac{ ( f(x_i) + f(x_{i+1}) ) }{ 2 } \Delta x \] \[ = \frac{ ( f(x_i) + f(x_i + \Delta x) ) }{ 2 } \Delta x \]
に変更しただけですね。

さて誤差は ?

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

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

上の図で、赤色の線(重なっていて見づらいですが…)が、今回の計算結果です。 青色の線は、厳密な値と比べるために、\( \sin x \) のグラフを重ね描きしたものです。 ほとんど重なっていますね。今は刻み幅 N = 10、つまりたった10個の台形で面積を近似しているだけなのに、 これだけ厳密な値と一致しているというのは、なかなか良い精度ではないでしょうか。

グラフに数式の値を重ね描きするには…
現在のVCSSL処理系で採用している2Dグラフソフトは「 リニアングラフ2D 」で、 データに数式を重ねてプロットする機能が付いています。上では、その機能を使いました。 まず、グラフ画面のメニューバーから 「 Tool 」 > 「 Math 」と選択します。 続いて表示されるウィンドウの「 f(<x>)」の入力欄に「 sin(<x>) 」と入力し、 「 PLOT 」ボタンを押します。 すると、グラフの上に、青い線で \( \sin x \) が重ね描きされます。

さて、前回の短冊法(長方形近似)のプログラムで同様に描いたグラフと、 右上の部分(最も厳密値からずれている部分)を比べてみましょう。

台形法と短冊法の比較( N = 10 )
台形法と短冊法の比較( N = 10 )
今回の台形法と、前回の短冊法(長方形近似)の結果の比較。青色の線が厳密値で、赤色の線が数値積分の値。共に刻み数 N=10 で、台形法のほうが誤差を小さく抑えられている事が見て取れます。

この通り、同じ刻み数 N = 10 で比べても、今回の台形法のほうが、 前回の短冊法(長方形近似)よりも厳密値に近い = 誤差を小さく抑えられている事が分かりますね。 予想通りの結果ではないでしょうか。

刻み幅を細かくしていくと…

誤差が小さく抑えられている事はグラフから大体分かりました。 もう少し厳密に、計算結果の数値から誤差を見てみましょう。 コンソールを見ると、以下のような内容が出力されているはずです:

0.0    0.0
0.1    0.0997502082639013
0.2    0.19850374541986468

0.8    0.7167581945005881
0.9    0.7826740283814796
1.0    0.8407696420884198

このように、a から b まで積分していっている最中の x 値が左の列、 その地点までの積分値が右の列に出力されています。 1行ごとに、積分値に台形1個の面積が足されています。 最後の行の積分値 0.8407696420884198 が、最終的に求まった \[ \int^{1}_{0} \cos(x) dx \] の値で、厳密な値は sin 1 なので 0.8414709848... です。小数点以下 3 桁目から食い違っていますね。 この食い違いが誤差です。

さて、台形の刻み幅を細かくするほど、つまり刻み数 N を大きくするほど、誤差は小さく抑えられるはずです。 実際にプログラムにおいて N を増やしていき、先ほどの最終行の値がどうなっていくか見てみましょう:

N=10の場合の結果      0.8407696420884198
N=100の場合の結果     0.8414639725380026
N=1000の場合の結果    0.8414709146853138
N=10000の場合の結果   0.841470984106668
N=100000の場合の結果  0.8414709848008824

わかりやすいように、厳密な値( \(\sin 1\) )から食い違っている部分を赤色で塗りました。 前回の短冊法(長方形近似)の場合と比べて、やはり誤差が急速に減っていっていますね。

上の結果から、台形法では、刻み数 N を 1 桁上げるごとに、 合っている桁数を大体 2 桁ずつ増やせる事が分かります。 こういうのを一般に「 2次の精度 」と言ったりします。 つまり台形法は 2 次の精度の数値積分アルゴリズムという事になります。

それに対して、前回扱った短冊法(長方形近似)では、刻み数 N を 1 桁上げるごとに、 合っている桁数は大体 1 桁ずつ増えました。 つまり短冊法は 1 次の精度の数値積分アルゴリズムという事になります。

数値積分に限らず、数値計算で使用する色々なアルゴリズムには、このような精度の次数がとても重要です。 一般に、次数が高いほうが、より誤差蓄積を抑えつつ、より高速に計算できるためです。 ただし、あまり次数の高いアルゴリズムは、逆に 1 刻みごとの計算量が増えて計算速度が低下してしまい、 刻み数を増やせなかったりもします。 また、計算の途中過程において、誤差がどのように振る舞うか(単調に増加していくのか、 ある程度の幅で振動するのか)という点も重要です。

今回のような練習ではなく、もっと本格的な数値計算の場においては、そのあたりの事情も総合的に考慮して、 どのアルゴリズムを採用するかを決定します。 個人的に知っている限りでは、 1 次の精度のアルゴリズム(短冊法もそうです)というのは、やはり精度や速度的なデメリットが大きく、 あまり実用されているのは見かけません(※)。 今回のような 2 次精度のアルゴリズムからが、ようやく実用範囲に入ってくるという感じがします。

※ ただ、短冊法のような低次精度のアルゴリズムに全く価値が無いかというと、そうでもありません。 個人的には、数値計算プログラムを書く際、最初はとりあえず簡単な( 間違える余地が少ない ) 低次精度のアルゴリズムを使って書きはじめます。そしてバグが無い状態の結果を確認した上で、 より複雑で高次精度のアルゴリズムに書き換える、という事をよくやります。 そうするとバグった際にすぐに分かるからです。

今回の台形法は、覚えやすく簡単なので、そこそこの精度でちょっと積分の値を求めたい、 くらいの場合には実際に結構使えると思います。 なんといっても、あのシンプルな短冊法(長方形近似)から数秒で書き換えられるのに、 精度は一気に跳ね上がるという、 そのコストパフォーマンスの高さはかなか侮れません。 たった一行の違いでこれだけ圧倒的な差が出てくるというのは、 数値計算プログラミングの醍醐味の一つではないでしょうか。

では、台形法を知っていればもう十分かというと… さらにほんの少しだけ書き換えるだけで、なんと 4 次の精度が出るようになります。 という事で、次回は簡単で高速ながらも 4 次の精度を誇る、シンプソン法の数値積分コードを書いてみましょう !

コードのライセンス

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



スポンサーリンク



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

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

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

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

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

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

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

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

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

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

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

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

お知らせ( Twitter )

スポンサーリンク