python練習問題 | √2などの平方根を計算する

今回は平方根の近似計算します。

2の平方根、1.41421356・・・ヒトヨヒトヨニヒトミゴロですね。

ニュートン法って、計算を繰り返して答えに落ち着くのを待つプログラムを作ってみます。

さらに、解の収束具合・計算回数・初期値選びについても見ていきます。

プログラミング課題

では、この式を計算するプログラムを作ってください。

ただし、

• 計算したxをvに入れなおして、またxを計算するってのを繰り返してください
• 最初にaには2、vには好きな値を入れてください
• 計算を繰り返す回数も好きにしてOKです

この式、あえて分子に()をつけなかったので、計算ミスがでるのもまた勉強です。

どんな値がでましたか？

プログラムは下記の通りです。  

v=3
a=2
for i in range(10):
  x=v+(a-v**2)/(2*v)
  v=x
  print(x)

1.8333333333333333
1.4621212121212122
1.4149984298948028
1.4142137800471977
1.4142135623731118
1.414213562373095
1.4142135623730951
1.414213562373095
1.4142135623730951
1.414213562373095

さて、この値はなんですかなぁ。

「ニュートン法」

「ニュートン法」（1669）という計算方法です。

この式は、コンピュータが登場する300年とか前に考えられているってのがすごいですね。

数式を繰り返し再計算して答えを出していく手法ってのがあるんですね。

私とニュートン法との初めての出会いは、人工衛星の軌道計算でした。

学生時代の授業で取り扱ったので、自分でC言語プログラミングしてみました。人工衛星の軌道は「軌道六要素」という6つの値で表現されています。

作ったプログラムは、軌道六要素と自分の位置（緯度と経度）と時刻を入力したら、人工衛星がどの方向にどの高さで見れるかというのを計算するものでした。方位角と仰角（ぎょうかく）って云います。

ニュートン法による数式計算もすごいなぁと思いましたが、行列による座標変換の大事さが分かりました。

人工衛星（だけではなく軌道上の物体ですが）の軌道要素は、NORAD（ノーラッド）って「北アメリカ航空宇宙防衛司令部」が公開しています。

ひょっとしたら、「クリスマスのサンタクロースの追跡」で知っているかもですね。

当時笑った「大統領が子供の夢を粉砕」したニュースも添えておきますね。

人工衛星の計算は、教科書に載ってた式を使っただけですので、難しいことはしていないです。

解の収束具合をみてみましょう

初期値vは、5でも100でも解は収束します。

# 初期値v=5の場合
2.7
1.7203703703703703
1.44145536817765
1.414470981367771
1.4142135857968836
1.4142135623730951
1.414213562373095
1.4142135623730951
1.414213562373095
1.4142135623730951

# 初期値v=100の場合
50.01
25.024996000799838
12.552458046745901
6.3558946949311395
3.335281609280434
1.967465562231149
1.4920008896897232
1.4162413320389438
1.4142150140500531
1.41421356237384

しかし計算が落ち着く（収束する）までの計算回数に差があります。

初期値は答えに近い値にすると良いですね。  

例えば、2の平方根だったら、1の平方根（つまり1）と4の平方根（つまり2）の間だから、初期値は1.5にしようかな？

例えば、5の平方根だったら、4の平方根（つまり2）と9の平方根（つまり3）の間だから、2.5を初期値にしようかな？

と、いった感じです。

こんな風に「初期値を決めてくれるアルゴリズム」を追加しても良いですね。

私達は「求めたい平方根a」だけを設定して、あとは全部プログラムがやってくれると。

計算回数を調整してみましょう

下記の計算結果は、平方根a=2を初期値v=3000で計算した結果です。

初期値vを3000に設定すると、計算回数10回では答えにたどり着きません。

1500.0003333333334
750.0008333331853
375.0017499984445
187.50354165344456
93.75710405931714
46.8892178881024
23.465935808627755
11.775582867900043
5.972712920318732
3.1537845641618323

計算回数を調整するように工夫したいですね。

「値の更新が0.00001以下になるまで計算をする、に改良してください」

# 解の改善が充分小さくなるまで続ける。
v=30000
a=2
before = v
diff   = 1 # 0.00001より大きければなんでもOK
cnt    = 0
while(diff>0.00001):
    cnt=cnt+1
    v=v+(a-v**2)/(2*v)
    diff = before - v
    print(f"{v} {diff}")
    before=v
print(f"cnt = {cnt}")

5000.000033333334 14999.999966666666
7500.000083333333 7499.99995
3750.0001749999983 3749.999908333335
1875.0003541666533 1874.999820833345
937.5007104165593 937.499643750094
468.751421874138 468.7492885424213
234.37784426393125 234.37357761020672
117.19318874685467 117.18465551707658
58.605127292538874 58.5880614543158
29.319626999702802 29.285500292836073
14.69392034575398 14.625706653948821
7.415015530226805 7.278904815527175
3.8423692493440917 3.5726462808827133
2.181440715408433 1.6609285339356585
1.5491329989172367 0.6323077164911965
1.4200888662915165 0.12904413262572012
1.414225716259008 0.005863150032508591
1.4142135624253205 1.2153833687467142e-05
1.4142135623730951 5.2225335167577214e-11
cnt = 19

計算回数が19回に自動的に調整されました。

ループ処理で、前の値を記録しておくという、よく使う処理を経験することができます。

初期値の話

今回は「初期値はなんでもよいよ」と言いましたが、負の値で試したなら、とってもセンスがありますね！

v=-10 ### 初期値を負の値にしてみる
a=2
for i in range(10):
  x=v+(a-v**2)/(2*v)
  v=x
  print(x)

-5.1
-2.7460784313725486
-1.7371948743795982
-1.444238094866232
-1.4145256551487377
-1.4142135968022693
-1.4142135623730954
-1.4142135623730951
-1.414213562373095
-1.4142135623730951

ということで、初期値って重要なんです。

特に非線形であれば、初期値によって導出される答えが変わる場合があります。

答えに近い初期値が良いですよね。答えがわからんから難しいわけですが。

難しい話ですが、AIのニューラルネットワークでも問題になります。

ニューラルネットワークには、脳みそネットワークの設定値「重み」が数千数万とあります。

良い脳みそになるにはこの重みを調整することで作っていきます。

ですから、重みの初期値をどうするかということが問題になります。