どうも！こんにちはのぶ(@tsureblo_nobu)です！
大学三年生になり研究室を決めない時期になってきたのですが、僕は理論物理系の研究室に行くことを選びました。
そこで必要になってくるのがプログラミングでの数値解析！
なのでこれから少しずつ僕が今まで学んできた数値解析のプログラムをコード例を利用して解説していこうと思います！
大学の授業で数値解析やプログラミング演習などがある方の参考になれば嬉しいです。
1.今回扱うコード例
とりあえず今回扱うコードを見ていきましょう。
「"""」で囲まれている部分や「#」がついている部分は実際にプログラムには関係のない注意書きです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
# -*- coding: utf-8 -*- """ bisec.pyプログラム 二分法による方程式の解法プログラム 使い方 c:\>python bisec.py """ # グローバル変数 a = 2 # f(x)=x*x-a LIMIT = 1e-20 #終了条件 #下請け関数の定義 #f()関数 def f(x) : """関数値の計算""" return x * x - a #f()関数の終わり #メイン実行部 #初期設定 xp = float(input("xpを入力してください:")) xn = float(input("xnを入力してください:")) #繰り返し処理 while (xp - xn) * (xp - xn) > LIMIT: #終了条件を満たすまで繰り返す xmid = (xp + xn) / 2 if f(xmid) > 0: xp = xmid else: xn = xmid print("{:.15f} {:.15f}".format(xn, xp)) #bisec.pyの終わり |
一見長く見えますが怖がらずに１つずつ見ていきましょう。
まず、最初の二行に注目
|
1 2 |
a = 2 LIMIT = 1e-20 #終了条件 |
これはx^2-aの計算を小数第20桁まで計算することを表しています。
関数の定義、初期設定に関しては簡単なので割愛します。
最後に注目したいのは繰り返し処理の部分。ここが一番大切！
|
1 2 3 4 5 6 7 |
while (xp - xn) * (xp - xn) > LIMIT: xmid = (xp + xn) / 2 if f(xmid) > 0: xp = xmid else: xn = xmid print("{:.15f} {:.15f}".format(xn, xp)) |
while文の部分で桁数が10^20を超えるまで以下の動作を繰り返すようになっている。
どこかで計算を切ってあげないと無限ループをしてしまうからLIMITを設けている。
では、次の行のxmidについてはどうだろうか？
xmidは適当に置いたある２つの数xn,xpの間の値をとる。
だからこの解の求め方は二分法と呼ばれています。
それではif文の中身を見ていきましょう！
xmidが0より大きければxpにxmidを代入、0より小さければxnに代入するようになっている。
これがどのような動作をすることになるのかわからない人は下の表を見て欲しい。
二分法のグラフ
これを見れば分かる通りif文の操作をすることで段々と解の位置が限定されてきていることが分かると思う。
ちなみに実行結果は以下の通り
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
xpを入力してください:1.5 xnを入力してください:1.3 1.400000000000000 1.500000000000000 1.400000000000000 1.450000000000000 1.400000000000000 1.425000000000000 1.412500000000000 1.425000000000000 1.412500000000000 1.418750000000000 1.412500000000000 1.415625000000000 1.414062500000000 1.415625000000000 1.414062500000000 1.414843750000000 1.414062500000000 1.414453125000000 1.414062500000000 1.414257812500000 1.414160156250000 1.414257812500000 1.414208984375000 1.414257812500000 1.414208984375000 1.414233398437500 1.414208984375000 1.414221191406250 1.414208984375000 1.414215087890625 1.414212036132813 1.414215087890625 1.414213562011719 1.414215087890625 1.414213562011719 1.414214324951172 1.414213562011719 1.414213943481445 1.414213562011719 1.414213752746582 1.414213562011719 1.414213657379150 1.414213562011719 1.414213609695435 1.414213562011719 1.414213585853576 1.414213562011719 1.414213573932648 1.414213562011719 1.414213567972183 1.414213562011719 1.414213564991951 1.414213562011719 1.414213563501835 1.414213562011719 1.414213562756777 1.414213562011719 1.414213562384248 1.414213562197983 1.414213562384248 1.414213562291116 1.414213562384248 |
段々と解が限定されていき1.4142......に近づいていっていることが分かると思う。
二分法ではこのようにして方程式の解を求める法則でした。
2.最後に
今回扱った二分法はpythonに存在するモジュールを利用すれば一瞬で解くことができます。
だからといって学ばなくて良いとないがしろにするのは勿体無いと思うので今回概要をまとめてみました。
参考になれば嬉しいです。