番外編 台形公式でπを計算する [数値解析]
ねこ騙し数学 番外編 台形公式でπを計算する
を、定積分の近似計算法である台形公式を用いて求めようって企画。
手計算では大変なので、
十進BASIC（プログラム言語であるBASICに準拠したもの）なるソフトを使って計算させるプログラムを作ってみた。
プログラムは以下の通り。
REM 台形公式を使って面積を計算する
REM 積分する関数↓の定義
DEF f(x) = 1/(1+x^2)
REM aは積分の下端、ｂは積分の上端
LET a = 0
LET b =1
REM 分割数ｎの入力
REM LET n = 10
INPUT PROMPT "分割数？" :n
REM 分割の幅hの計算
LET h =(b-a)/n
REM メインの計算
LET sum =0
FOR i = 1 TO n-1
LET x = a + i*h
LET sum = sum + f(x)
NEXT i
LET s1 = f(a) + sum
LET s2 = f(b) + sum
LET s1 = s1*h
LET s2 = s2*h
LET s = (s1+s2)/2
REM 計算結果の出力
PRINT "π/4 ="; PI/4
PRINT "台形公式 S = ";s
PRINT "Σf(x_k-1)*h =";s1
PRINT "Σf(x_k)*h = ";s2
END
もっと簡単に書けるんだけれど、分かりやすくするためにこう書いたケロ。
台形公式とは
台形公式というのは、
a ≦ x ≦ bのとき、つまり、x∈[a,b]だね、これをn個に分割して（分割の幅h = (b – a)/n ）、
を計算するもの計算する方法。
として、この区間の面積を
と考え、これをk =1 ～ nまで足し合わせたってもの。
台形の面積は
台形の面積 ＝ （上辺+底辺）×高さ÷２
でしょ。
になっていることが分かるケロ。
①を使ってプログラムを作ってもいいんだけれど、これは計算の無駄が多いので、
だから、
として計算している。
で、Ｓ１とＳ２をみると、
となって、
は共通している（プログラム中のsumはこれ）。
だから、台形公式は
と書いてもいい。
なんで、台形公式によるＳの他にS1とS2を求めているかというと、これも近似だから。
そして、
は減少関数だから、
リーマン積分の上限の寄せ集めＳ（このＳは台形公式による計算値ではないので注意してください）がS1で,下限の寄せ集めｓがS2になる。
この場合、当然のことながら、S1 > S （このSは台形公式の計算値） > S2になる。
そして、π/4は、S1 > π/4 > S2になる。
計算結果ケロ
ちなみに、分割数ｎ = 10で計算した結果は、以下のとおり。
分割数？10
π/4 = .785398163397448
台形公式 S = .78498149722679
Σf(x_k-1)*h = .80998149722679
Σf(x_k)*h = .75998149722679
分割数を100にすると、
分割数？100
π/4 = .785398163397448
台形公式 S = .785393996730786
Σf(x_k-1)*h = .787893996730786
Σf(x_k)*h = .782893996730786
となる。
台形公式、結構、精度がいいでしょっ。
S1（＝Σf(x_k-1)*h）とS2（＝Σf(x_k)*h）は、
ｎ = 10 のとき、かなり誤差が大きいのに、
これを平均する、つまり、この二つを足して２で割るだけで、とたんに精度がよくなる。
何とも不思議な話だとは思いませんか？
でも、
を計算させると、
分割数？10
π/4 = .785398163397448
台形公式 S = .776129581562081
Σf(x_k-1)*h = .826129581562081
Σf(x_k)*h = .726129581562081
分割数？100
π/4 = .785398163397448
台形公式 S = .78510425794476
Σf(x_k-1)*h = .79010425794476
Σf(x_k)*h = .78010425794476
と精度がかなり落ちる。
分割数ｎを大きくしないと、いい値が出ない。
分割数を100万くらいにすると、
分割数？1000000
π/4 = .785398163397448
台形公式 S = .785398163103477
Σf(x_k-1)*h = .785398663103477
Σf(x_k)*h = .785397663103477
となる。
1000万だと、
分割数？10000000
π/4 = .785398163397448
台形公式 S = .785398163387691
Σf(x_k-1)*h = .785398213387691
Σf(x_k)*h = .785398113387691
予期せぬ計算結果にネムネコ、当惑する。
やけに精度がいいな（ポリポリ）。
普通、分割数をあまり大きくしすぎると、逆に精度が落ちるものなんだが・・・。
「これは有限桁で計算せざるを得ないというコンピュータの宿命みたい」なことを書こうと思ったんだけれど、
これでは、このことを書けない（ポリポリ）。
BASICは、精度がいいんですよ（笑）。
与太話
どうです、十進BASICをパソコンにインストールしたくなりましたか？
容量は、ネムネコのパソコンに入っているのが、1.3MBと至ってコンパクトです。
ダウンロードしても、これは損のないソフトだと思いますよ。これをインストールするだけで、かなり世界が広がるのは間違いがない。
この他に、シンプソン法という、もっと精度のいい定積分の近似計算法もあるんだけれど、
それは、また、後で。
ＣやＪａｖａで書いてもいいんだけれど、
これらのプログラミング言語はプログラムを書けるようになるまでべら棒な時間を要する上に、
これらのコンパイラをパソコンにインストールするのも大変だし、
この後に、パソコンの環境変数などを自分で設定しないと駄目なんで、ちょっと敷居が高い。
環境変数の設定をしそこなうと、
最悪、パソコンを使えない、アプリケーションソフトを使えなくなるなどの不都合が生じたりする。
ネムネコのせいでパソコンが使えなくなった、なんていわれると困るし・・・。
その点、十進BASICは安全ですよ。
数学カテで回答者をやっていたとき、回答者の中にも十進BASICを愛用している人が何人かおりました。
ねこ騙し数学で掲載しているグラフの半分くらいは、この十進BASICに書かせたもの。
ちょっとした計算をしたいとき、グラフや関数の概形を知りたいときなど、結構、重宝します、このソフト。
タグ：微分積分
コメント 0