夏休みの課題4のヒント


おしらせ

課題4のヒントとBASICのプログラム誤りがありました.8月31日PM6:00に修正しましたので、それ以前にこのページをプリントしたものを持っている人は、もう一度このページをプリントアウトしてください.


はじめに

この課題では、基本的に計算方法は各自で考えることになっていますが、それはちょっと難しいかもしれないので、2分法と呼ばれる計算法を説明します.もちろん、自分で計算法を考えた人はそれを使って BASIC のプログラムを作ってください.たとえその計算法が間違っていても、そちらの方が評価は高いです.

2分法

2分法は連続関数の零点(関数の値が 0 となる点)を求める計算法です.今、次のような連続関数 f(x) があり、f(a1)<0, f(b1)>0 (ただし、a1=0.5, b1=0.9) であるとします.

f(x) は連続関数ですから、区間 a1<x<b1 に f(x) の零点が必ず1つはあるのですが、上の図のように、f(x) の零点は1つしかないとします.つまり、 a1 と b1 は f(x) の零点を上の図のように挟んでいます.このとき、次のような ai, bi (i=2,3,4,...,n) を計算していきます.

2番目の条件より、ai と bi は i を大きくするに従い、ある値に近づいて行きます.1番目の条件より、その値は f(x) の零点に他なりません.よって、i を十分大きく取れば ai, bi は f(x) の零点に十分近い値になります.

説明だけでは、わかりにくいかもしれませんが、そういうものかと思ってください.実際に2分法を用いて上のグラフの f(x) の零点を計算してみます.

まず、a1、b1 の中点 (=(a1+b1)/2) での f(x) の値 f((a1+b1)/2) を計算します.このとき、

今の場合、(a1+b1)/2=0.7 ですから、下のグラフより f((a1+b1)/2)<0 です.よって、a2 = 0.7, b2 = 0.9 となり、グラフは次のようななります.

ここで、次の2点に注意してください.

つまり、はじめは f(x) の零点が区間 0.5 < x < 0.9 に1つあるということだけがわかっていましたが、上の操作でその零点は区間 0.7 < x < 0.9 にあることがわかり、零点がある区間を縮めることができたのです.2分法では、このようにして零点がある区間を 1/2 ずつ縮めていく操作を繰り返します.上の操作を一回行うごとに区間の長さが 1/2 になっていくのですから、i 回行えば区間の長さは

となります.つまり、

ですので、上の(条件2)が成り立っているのがわかります.よって i を十分大きくしてやれば、零点が十分な精度で求まることになります.

もう少し、上の操作を続けてみましょう.今、a2 = 0.7, b2 = 0.9 ですから、今度はその中点での f(x) の値 f((a2+b2)/2) = f(0.8)を求め、先ほどと同様に

とすれば、f(x) の零点は区間 a3 < x < b3 に含まれることになります.下のグラフより f((a2+b2)/2)=f(0.8) は正ですから、a3=0.7, b3=0.8 となります.

上のグラフからわかるように、前と同様に

以下、同様に

とする操作を i=1,2,3,...,n まで続けていけば、f(x) の零点は区間 ai < x < bi に含まれ、かつ

であることがわかりますので、上の(条件1),(条件2)を満たす ai, bi が計算されます.

以上で f(a1)<0, f(b1)>0 の場合の計算法を示しました.f(a1)>0, f(b1)<0 の場合には -f(x) の根を計算すれば良いので、次のようにすれば良いのです.まず一番はじめに f(a1) の値を計算し、それが正ならば c=-1、負ならば c=1 とおきます.後は c*f(x) の根を計算すれば、常に c*f(a1)<0 となって上の計算法で f(x) の根が計算されます( f(x) と c*f(x) の根は同じであることに注意してください).

BASICのプログラム

以上で、2分法はわかりましたが、実際の BASIC のプログラムはどうなるのでしょうか?
ここでは、プログラムの大まかな流れだけを示します.実際のプログラムは各自が作ってみてください.

[ステップ1] n に繰り返しの回数を代入し、配列 a(i), b(i) (i=1,2,3,...,n)を宣言する.

[ステップ2] a(1), b(1) を入力する.

[ステップ3] もし f(a(1))>0 ならば c=-1 そうでなければ c=1 とする.

[ステップ4] 以下の処理を i=1,2,...,n-1 に対して行う.

[ステップ 4-1] c*f((a(i)+b(i))/2) の値を計算する.

[ステップ 4-2] c*f((a(i)+b(i))/2) が負ならば、

c*f((a(i)+b(i))/2) が正ならば、

とする.c*f((a(i)+b(i))/2)=0 ならば、求める零点は (a(i)+b(i))/2 であるので、これを表示してプログラムを終了する.

[ステップ5] (a(n)+b(n))/2 を表示して、プログラムを終了する.

上の解説では良くわからない人が多いようなので、もう少し具体的に説明します.

f(x)=x^3+2*x^2-x-1 の 0 < x < 1 にある根を計算してみましょう.

上のプログラムのステップ1に相当する部分は BASIC でかくと

10 input "n=",n

20 dim a(n),b(n)

となりますね.10行目で繰り返しの回数 n を入力し、20行目で配列 a(i), b(i) を宣言しています.これで a(1),a(2),...,a(n), b(1),b(2),...,b(n) が変数として使えるようになります.

ここでは、0 < x < 1 にある根を求めるのですから、上の解説の a(1)=0, b(1)=1 となります.ですから、次に上のプログラムのステップ2に相当する部分は

30 a(1)=0

40 b(1)=1

となります.

ステップ3では f(a(1)) を計算する必要がありますが、ここでは

f(x)=x^3+2*x^2-x-1

ですから、

f(a1)=a(1)^3+2*a(1)^2-a(1)-1

ですね.よって上のプログラムのステップ3に相当する部分は

50 if (a(1)^3+2*a(1)^2-a(1)-1)>0 then c=-1 else c=1

です.

ステップ4では、ステップ4-1, ステップ4-2 を i=1,2,...,n-1 に対して行いますから

60 for i=1 to n-1

ですね.

ステップ4-1 で c*f((a(i)+b(i))/2) の値を y に代入すると、ステップ4-2 は

90 if y<0 then goto 100 else goto 130

100 a(i+1)=(a(i)+b(i))/2

110 b(i+1)=b(i)

120 goto 150

130 a(i+1)=a(i)

140 b(i+1)=(a(i)+b(i))/2

150 next i

となります.ステップ5では

160 print "kon=",(a(n)+b(n))/2

とすれば良いですね.

以上で、プログラムのほとんどの部分が出来ました.あとはステップ 4-1でc*f((a(i)+b(i))/2) の値を y に代入するところを作れば完成です.

ここの部分は各自が作ってみてください.

自分が作ったプログラムがどうして動かないのかどうしてもわからないときは、kitamoto@edu.yamaguchi-u.ac.jp にプログラムを添付して、e-mail をください.できる範囲で答えます.それでは、がんばって課題を行ってください.