黄昏より暗きもの、血の流れより赤きもの

読者です 読者をやめる 読者になる 読者になる

黄昏より暗きもの、血の流れより赤きもの

自分の好きな事を好きなように書いて行きます。

【数値計算】シンプソンの台形公式による定積分の値の計算【C言語】

微分積分 数値計算 C言語

今日はシンプソンの台形公式による定積分の値の計算法について話そう。とその前に定積分がどのようなものかを確認していく。

リーマン積分とは?

この世には様々な積分と言うものが存在する。その内日本の高校生や大学生が普段扱っているものとしておなじみなのがリーマン積分だ。以下、有名な解説本として名高い高木貞治の解析概論 改訂第三版(軽装版 岩波書店)のP91からP95までの内容を端折って説明していく。

ここに関数y=f(x)上の区間[x_i,x_{i+1}]を任意の細区間に分ける。さて各細区間毎の長方形の面積を考えてみる。このとき横の長さはδ_{i+1}-δ_i、高さはf(ε_i)となる。


f:id:program_study:20160204212416j:plain

ここで、各細区間毎の間隔δを限りなく小さくし、δ→0としよう。長方形の幅δと高さεの選択によらず、極限値I = \lim{δ→0} \sum{δ_i f(ε_i)}が存在する事積分可能と言う。さらにこの事を区間[a,b]におけるf(x)の定積分と言い、I =\int_a^b f(x) dxで表す。

シンプソンの台形公式

微分積分微分方程式数値計算法の必要性について

上記のリーマン積分I =\int_a^b f(x) dxの定義からも分かる通り、[a,b]の分割の仕方を問わず原理上定積分の値を求める事が出来る。

しかしながらプログラムのfor文ループで扱う時に不便な他、そもそもfor文ループでは連続した数量を扱う事が出来ない。この為何らかの代替案により計算する事となる。

これは解が存在しない微分方程式微分の場合も同様で、何らかの手法によりその代替案が多数提案されている。例えば画像処理の輪郭(エッジ)検出に微分の概念が使われており、微分数値計算の代替案としてSobelオペレーター*1と言ったアルゴリズムも存在する。

シンプソンの台形公式について

さてシンプソンの台形公式に話を戻そう。シンプソンの台形公式は区間[a,b]をn等分し、n等分された区間(x_1,x_2,….,x_n)の台形の面積の和を定積分の値として代用する方法である。やや語弊はあるものの、丁度数IIIの教科書でおなじみの区分求積法が長方形ならば、それを台形に差し替えて計算する方法と言って差し支えない。


f:id:program_study:20160206083901p:plain

さてk(0 \leq k \leq n)区間目における台形の面積S_kを考える。

まず「台形の面積の公式=1/2*(上底+下底)*高さ」であり、計算のためこれを文字式に置き換えていく。

まず区間[a,b]をn等分するので、高さはh = \frac{b-a}{n}となる。さらに上底の長さがf(x_i)、下底の長さがf(x_{i+1})となる。以上からk区間目における台形の面積はS_k = \frac{1}{2} * (f(x_i)+f(x_{i+1})) * \frac{b-a}{n}

最後に台形の面積S_k0 \leq k \leq nの区間で集計して定積分の代用値Sは以下の式で計算できる。

S=\sum_{k=0}^n {S_k} = \sum_{k=0}^n {\frac{1}{2} * (f(x_k)+f(x_{k+1})) * \frac{b-a}{n}}

シンプソンの台形公式のプログラム(C言語)

以上を踏まえたC言語のプログラムは以下となる。

テスト

以下、作ったプログラムのテスト風景を紹介したい。まず適当な関数y=f(x)を決め、以下の事柄について調べる。

事柄1:プログラムの計算結果と手計算の結果が合っているか?
事柄2:プログラムの計算結果(実験値)と手計算の結果(理論値)とでどれだけ誤差がでるのか?

プログラムの正誤の確認はもちろんの事、今回は代用値(近似計算)を扱っているため、手計算による理論値とどれだけ誤差が出るのか?について確認しておくと良い。

調べる事柄1:プログラムの計算結果と手計算の結果が合っているか?

まずはプログラムの正誤を調べたい。正誤を調べるコツは、できるだけ簡単な事例を用意して調べると言う事だ。実際問題自分でテストケースを作成する必要も有るため、プログラムの概要に合わせてテストケースを作成できるようにしておきたい所だ。

テストケースの作成(手計算による計算結果)

今回はテストケースとして関数f(x) = x^2の区間[1,3]の定積分の値\int_1^3 x^2 dxを考える。シンプソンの台形公式で求めた定積分\int_1^3 x^2 dxの代用値をAとする。Aを求める上で台形を分割する区間をn=2として考えるため、高さh = \frac{3-1}{2} = 1と計算出来る。以上を踏まえると以下の値を得る。

手計算の結果

A =\sum_{k=0}^1 \frac{1}{2} (f(x_{i+1})+f(x_{i})) h
A =\frac{1}{2} *(f(1)+f(2)+f(2)+f(3))
A =\frac{1}{2} *(1+4+4+9)
A = 9

プログラムに依る結果

これに対しプログラムの結果は以下の通り。

XXX$ gcc simpthon.c
XXX$ ./a.out
途中の計算結果=S(1.000000)=2.500000
途中の計算結果=S(2.000000)=9.000000
最終的な計算結果=9.000000

結果

Aの値とプログラムに依る結果一致しているので、期待すべき結果が得られた。

調べる事柄2:プログラムの計算結果(実験値)と手計算の結果(理論値)とでどれだけ誤差がでるのか?

概要

今度は分割数nn = 10n = 100,n = 1000と増やしていき、分割数nを増やすと理論値との誤差がどれだけ変化するかを確認する。ここでnの値を変えて検証する理由は、その値に応じて結果が変わることが予想されるからである。

これは丁度n人の人に「内閣を支持する」「内閣を支持しない」とアンケート調査する上で、n=100人に対し調査するのと、n=1,000人に対し調査するのとでは、最終的な支持率の値に違いが生じると言う現象に似ている。

テストケースの作成(手計算による計算結果)

今度もテストケースとして関数f(x) = x^2の区間[1,3]の定積分の値\int_1^3 x^2 dxを考える。高校の教科書通りの計算で求めた定積分\int_1^3 x^2 dxの理論値をBとして、以下を得る。

手計算の結果

B =\int_1^3 x^2 dx = \frac{27}{3} - \frac{1}{3} = \frac{26}{3} = 8.666

プログラムに依る結果

これに対しプログラムの結果は以下の通り。

n=10のとき=8.680000
n=100のとき=8.666800
n=1000のとき=8.684680

結果や考察

2つの結果を見比べるとn=100のとき、B=8.666に最も近い値が得られる。しかしながらn=100の場合とn=1000の場合を見比べると、n=1000の場合の値の方が、n=100の場合の値よりも高い値が観測された。

これは多くの台形に分割することで、関数y=f(x)に収まりきらない部分の面積が増えたためだと考えられる。以上より分割数nを増やしたからと言って、誤差の少ない代用値を得る事ができるとは限らない事が分かる。

参考書籍

一冊目の本は自分が高校時代に背伸びして買った本で、数学科ってこういう事をやるのねと悟った本であった。まさかこのブログの説明の参考資料として活躍する日が来るとは思いもしなかった。

解析概論 改訂第3版 軽装版

解析概論 改訂第3版 軽装版

Cによる数値計算とシミュレーション

Cによる数値計算とシミュレーション

チャート式基礎からの数学3―新課程

チャート式基礎からの数学3―新課程