NEW GAME!を観ました。プログラミング始めます。（19）

初めてコメントを頂きました！本当にありがとう！！
f:id:kirintt:20171027040636p:plain

※基本的に、愛知大学の有澤先生がプログラミングの講義で使用しているテキスト「Pythonによるプログラミング入門」に沿って進めています。
※Python 2.7.13を使用しています。

今回のテーマは「グラフの交点」です。

グラフの交点

この節ではグラフの交点を求める簡単なアルゴリズムについて勉強する。

　ここに紹介するのは二分法である。交点の x 座標が x0 と x1 の間に入るような x0 と x1 を最初に適当に選ぶ。但し x0 と x1 の間には交点は1つだけ存在しなくてはならない。また話を簡単にするために x0 < x1 とする。

（引用元：Pythonによるプログラミング入門, p. 46）

ふむふむ…。
とりあえずどんなプログラムか見てみよう。

譜 4.8 ２分法を使った交点の計算

#
# root of a function
# Bisection method
#
from math import *
def f(x):
    return cos(x)-x
x0=0
x1=pi/2
while x1 - x0 > 1.0e-8:
    x = (x0+x1)/2
    y = f(x)
    print x,y,x0,x1
    if y > 0: x0=x
    elif y < 0: x1=x
    else: break

（引用元：同上, p. 46）

順番に何をしているか書いてみる。

5行目：pi(円周率)を使う準備。
6, 7行目：関数f(x)の定義。f(x)が0のとき、y=cos(x)とy=xのグラフは交差する。
8, 9行目：最初に見る範囲の設定。
10行目以降： x1 - x0 > 1.0e-8 の間繰り返す。
11行目：変数xにx0とx1を2分する値を入れる。
12行目：変数yにf(x)を入れる。xは11行目で代入された値。
13行目：いくつかの変数を表示。計算が正しく行われているか確認できる。
14, 15行目：yの正負によって、次に狭める範囲がx0かx1かを判断する。
16行目：y=0なら、交点が見つかったということなので計算終了。

結果

0.785398163397 -0.0782913822109 0 1.57079632679
0.392699081699 0.531180450813 0 0.785398163397
0.589048622548 0.242420989754 0.392699081699 0.785398163397
…（中略）…
0.739085137954 -7.9310493728e-09 0.739085126251 0.739085149657
0.739085132102 1.86237991695e-09 0.739085126251 0.739085137954

問１

次の式で表される 2 つのグラフ
　　　　y = − x ** 2 + 1
　　　　y = x
の交点を求めなさい。(2 つとも求めなさい)
（引用元：Pythonによるプログラミング入門, p. 47）

プログラム

#
# root of a function
# Bisection method
#
# from math import *
def f(x):
    return -1*(x**2)+1-x
x0=0
x1=1.0
x00=x0
x2=-2.0

print "When x>0"
while x1 - x0 > 1.0e-8:
    x = (x0+x1)/2
    y = f(x)
    n1 = x
    print x,y,x0,x1
    if y > 0: x0=x
    elif y < 0: x1=x
    else: break
print   # new line
print "When x<0"
while x00 - x2 > 1.0e-8:
    x = (x00+x2)/2
    y = f(x)
    n2 = x
    print x,y,x2,x00
    if y > 0: x00=x
    elif y < 0: x2=x
    else: break
print   # new line
print "result ("+str(n1)+","+str(n1)+"),("+str(n2)+","+str(n2)+")"

譜 4.8 からの変更点・追加点

全体：交点が2つ存在するため、x = 0を境にx >= 0 と x <= 0 の2つに分けて考える。（交点がx >= 0, x <= 0 それぞれの範囲に1つずつあることが前提）
5行目：ここではpiを使わないのでコメントアウト。
9行目：小数点以下まで計算するために".0"を追加
10, 11行目：x00, x2 は x <= 0 の範囲での交点探索に使用する。
13, 23行目：2つの交点を求める過程が表示されるので、どっちがどっちか分かるようにするための文。
17, 27行目：n1, n2 にそれぞれ交点の x 座標を入れる。（最後の結果表示で使用）
22, 32行目：結果を見やすくするために空白行を設けた。
29, 30行目：x <= 0においては、x > (交点) において y > 0, x < (交点) において y < 0 という具合で、x >= 0 の場合と条件が逆転する。これについては y = − x ** 2 + 1, y = x のグラフを1つの図に書いてみるとわかりやすいはず。

33行目：2つの交点を表示。str(n1)は変数n1を文字として表示します、という意味。こうしないと、「文字と数値は"+"でつなげないよ！」といった感じのエラーが出た。

結果

When x>0
0.5 0.25 0 1.0
0.75 -0.3125 0.5 1.0
（中略）
0.618033990264 -3.38550543155e-09 0.618033975363 0.618034005165
0.618033982813 1.32744992776e-08 0.618033975363 0.618033990264

When x<0
-1.0 1.0 -2.0 0
-1.5 0.25 -2.0 -1.0
（中略）
-1.61803399026 -3.38550520951e-09 -1.61803400517 -1.61803397536
-1.61803398281 1.32744992776e-08 -1.61803399026 -1.61803397536

result (0.618033982813,0.618033982813),(-1.61803398281,-1.61803398281)