このエントリーをはてなブックマークに追加
2015/06/25

台形近似による数値積分

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

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

概要

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

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

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

[ 関連記事 ]

使用方法

ダウンロードと起動

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

なお、ダウンロードボタン横のリンクから、Webページ上でそのまま実行する事も可能です。

起動後

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

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

被積分関数や積分区間、精度などを変えたい場合は、 プログラム「 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 \]
前回からほんの少し変わっただけですね。実際、プログラムに加える変更もわずかで済みます。

コード解説

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

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

前回との違いは 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言語版のコードも同様です。 そのままでのご利用はもちろん、言語の種類を問わず、改造や流用などもご自由に行ってください。



スポンサーリンク


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

台形近似による数値積分
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