数値計算のプログラム

メニュー

カレンダー

サイトの内容
├授業のページ
├Hyperresearch
├blog:ちょっと聞いてよ
├研究ノート
├非公開ノート
├国語
├フローサイトメトリー支援
├Pythonバイオ
└その他
　├内部
　├免責事項
　├プライバシー
　└転載

サイトマップ

最新の10件

メニュー編集

山内の授業のページ/07春/研究室ゼミ?

数値計算のプログラムで遊ぼう †

課題１　　円の面積を求めて円周率を計算してみよう †

円周率の計算方法はいろいろと効率のよい方法がありますが、ここでは円の面積を（定）積分によって求めて、そこから円周率πを求めてみます。積分の計算をどうやるか、がここの課題です。

対象とする図形は、半径1の1/4円にしましょう。　完全な円でもいいですが、関数の場合わけをすることになるので、最初から1/4円で考えます。

まず、円の形の関数形を確認しましょう。　原点を中心とする半径1の円は、x^2+y^2=1（x^2の^はべき乗を表すものとします。つまりx^2ならxの2乗、x^3ならxの3乗です）で表されますから、この式をyについて解くと、y=sqrt(1-x^2) となります。ここで関数sqrtは、いわゆる「ルート（平方根）」ですが、ホームページ上で書くのは難しいので、sqrt関数として書いておきます。また、この円の関数全体 sqrt(1-x^2) を f(x) と書くことにします。

この関数をｘ軸で区間０から１までを積分すると、1/4円の面積になります。では、その積分の仕方を考えてみましょう。最初に考えるのは長方形（矩形）です。

ｘ軸で区間０から１までをＮ個に分割し、図のように縦に長い長方形を作ります。この長方形の面積をＮ個分寄せ集めて、この1/4円の面積を計算しようというわけです。Ｎが大きくなれば、本物の円と長方形の寄せ集めとの差が小さくなるので、Ｎを無限に大きくすればいずれは1/4円の面積に近づくだろう、というのが積分の原理でしたよね。

さて、長方形の面積を計算しましょう。便利のために、（左から）ｉ番目の長方形の、左下の隅の点のｘ座標をxiとします。 xi = i/N です。　この点を左下の頂点とする長方形をRiと名づけることにしましょう。この長方形Riの高さは、左上の隅の頂点が円と重なっています。だから、円の関数y=f(x)を用いると、高さyi = f(xi) = sqrt(1-xi^2) となります。繰り返しますが、xi^2はxiの2乗のことです。従って、長方形Riの面積Siは、幅が1/Nで、高さがyiだから、Si = (1/N) ×yi で求められます。

Si = (1/N)×yi = (1/N) ×f(xi) = (1/N) ×sqrt(1-xi^2) = (1/N) × sqrt(1-(i/N)^2)

そこで、このSiを、i = 0 から i = N-1 まで、加え集めれば、面積が計算できるというわけです。これが積分の原理ですよね。

加えた結果は1/4円なので、それを4倍すれば半径1の円の面積になります。半径rの円の面積はπr^2で、半径r = 1ですから、π×1^2 = πになるはずですよね。

ここまでで、計算の仕方（つまりアルゴリズム、どうやって計算すればよいか、計算の手順）は分かったと思います。これをプログラムに書き表してみることにしましょう。

注意1）　xやy、1/N、S などの値は、小数になるので、double型で宣言しましょう。　また、i の値は数え上げる都合があるので、仮に整数にしておきます。ということは、x = i/N を計算するときに、 i と N が int型のままで計算すると、右辺を int型で計算してから double型の x に代入しますから、0 になってしまいます。　i と N をあらかじめdouble型にキャストしておいてから、割り算をすること。

注意２）　Ｃ言語のprintfは、フォーマットを指定する文字列で、何桁をプリントするかが決まります。　内部の計算は、double型の変数を使っているので、「仮数部」が48ビットでしたっけ？　だとすると、2の48乗は10の14乗ぐらいですから、13桁とか14桁ぐらいの精度があるはずです。　printfでそのぐらいの桁数を印刷してみましょうか。　　printf("%f16.14\n", S) などとすると、１６桁のスペースの中に少数点以下１４桁を表示するので、上に３と小数点の２文字を追加して、ちょうどになります。

　　　　　　＞＞矩形近似をプログラムにする

プログラムができたら、Nの値（つまり積分の区間0〜1を何等分するか）をいろいろ変えてみましょう。Nを大きくすれば、長方形は細くなって、演習より上にはみ出している部分は小さくなるでしょう。だからNを大きくすれば大きくするほど、長方形の面積Siの和は本当の1/4円の面積πに近づくはずです。Nの値を変えると、どれだけπに近づくでしょうか？試してみてください。

どのような結果が得られましたか？Nをどれだけにすると、結果が我々の知っているπの値にどれだけ近づくでしょうか？

頭を使う問題です。Nをいくつにすると、真のπの値からのずれ（誤差）はどのぐらいになるか、理論的に考えてみてください。つまり、長方形が円よりも上にはみ出している部分の面積を、およそ見積もってみてください。（たとえば三角形として見積もるとどうなりますか？）

課題１の改良　　　長方形が円より上にはみ出しているのを何とかしたい †

いくつかの改良が考えられるでしょう。自分でも考えてみてください。

改良１ †

もっとも単純なのは、長方形の高さを左端に合わせるのではなく、長方形の幅の真中（前の xi の値に長方形の幅の半分を加えた点）　で高さを計算することができるでしょう。いい加減に見積もっても、今までの方法は円より長方形が一方的に上にはみ出していたわけですが、今度の方法は左半分が下にはみ出すのに対して右半分が上にはみ出していますね。お互いに相殺しそうです。（円でなくて直線なら間違いなく相殺しますね）

とにかく、計算してみましょう。同じＮに対して、こちらのほうの値はどうなっているでしょうか？

改良２　台形公式 †

長方形の代わりに台形を使う「台形公式」が使われることがあります。　上の図の長方形で、左側の高さを下底とし、右側の縦の辺が円弧と交わる点のy座標を上底とする台形を考えます。　要するに、長方形にかかっている円弧を、直線で近似するわけです。　左の図のようになりますが、図では円弧と直線がほとんど重なって見えるので、あまりよく分からないかもしれません。

この面積を計算し、全体を足して見ます。単純には、下底が f(xi)、上底が f(xi+1)ですから、これを用いて（（上底）＋（下底））×（高さ）／２　とすれば求まります。

まずはこの方法で計算してみましょう。　前のプログラムをコピーして、修正してみてください。

世の中では、この台形で近似するやり方（台形法）を、式をもう少し整理してから計算します。というのは、i番目の台形の上底（右側の縦の辺）とi+1番目の台形の下底（左側の縦の辺）は同じ長さですから、２回計算する必要は無いわけです。

　　Si
= Σ((f(xi)+f(xi+1))*(1/N)/2)　　--- (上底+下底)×高さ/2
= (1/2N) * (f(x0) + f(xN-1) + 2*Σi=1N-2f(xi) )　　--- 上に述べた理由で２回計算するのでそれをまとめた

式変形は、自分できちんとフォローしておいて下さい。そのまま信じないように。
最後のシグマの所は、i=1からN-2までの和です。最初と最後を除いた和ですが、0から始めているので、1からN-2までと書いています。　この式が「台形公式」と呼ばれています。

では、この台形公式を用いて計算してみてください。　この形になっていると、このままプログラムに直すことができます。Σのところをループに書き直すだけです。やってみてください。

台形で近似すると、長方形（矩形）で近似したときよりは円周に近くなっている感じがしますね。実際にどれだけ違うか、長方形の時に使ったと同じNの値で、面積を計算してみてください。また、台形法の場合もNを大きくすると近似がよくなることも確かめてください。　その近似の度合いが、同じNの増やし方に対して、矩形法より台形法の方が早くπに近づくということが出来ます。

先に述べたように、単純に計算するより公式を用いたほうが、計算の回数が減っています。　単純な計算では、i=1からN-2までの場合について高さ（＝関数fの値）を２回計算しています。公式では、１回計算して２倍しています。公式のほうが計算が半分で済むので、関数の計算が複雑な場合には、かなりの節約になります。