Pythonで非線形方程式を解く (二分法、Newton法)
Pythonでざっくりと数値解析の基礎をしてみるシリーズの2つ目です。 snova301.hatenablog.com
方程式の解法
非線形方程式を解く基本的な考え方は、反復計算を高速に、正確に、計算することで、近似解を求めるということです。
代表的な方法として、二分法(Bisection method)とNewton法があります。
二分法
方程式が連続しており、かつ、ある区間において一つだけ解をもつことがわかっている場合に使用できます。
手順ですが、区間両端から加重平均を取り、その平均値がどちらの端に寄っているかを調べていきます。
そして、両端の範囲を小さくしていくことで、近似解を求めます。
こちらの記事に詳しい説明があります。
具体的に、2x2-3x-4=0を解いてみます。
プログラムは以下
# -*- coding: utf-8 -*- # import numpy as np import matplotlib.pyplot as plt def main(): min_x = 2 max_x = 5 x = np.linspace(min_x, max_x, 100) y = func_y(x) # もとの関数 # ここから二分法 for i in range(50): mid_x = (min_x + max_x)/2 # 平均値 mid_y = func_y(mid_x) # そのときのyの値 if min_x*mid_y < 0: # min_xとmid_xの区間に解が存在していなかった場合 min_x = mid_x else: # min_xとmid_xの区間に解が存在していた場合 max_x = mid_x ans_x = (min_x + max_x)/2 # 近似解 ans_y = func_y(ans_x) # そのときのyの値(確認用) print('x : '+str(ans_x)) print('y : '+str(ans_y)) # 以下プロット plt.plot(x, y, color='b', label='function') plt.plot(ans_x, ans_y, 'o', color='r', label='solution') plt.legend() plt.show() def func_y(x_dat): func = 2*x_dat**2 - 3*x_dat - 4 return func if __name__ == '__main__': main()
結果
x : 2.3507810593582117 y : -1.7763568394002505e-15
また、グラフは
ちなみに、Scipyに二分法のモジュールがありました。
以下は、それを使ったプログラム
# -*- coding: utf-8 -*- # import numpy as np import matplotlib.pyplot as plt from scipy.optimize import bisect def main(): min_x = 2 max_x = 5 x = np.linspace(min_x, max_x, 100) y = func_y(x) # ここから二分法 ans_x = bisect(func_y, min_x, max_x) ans_y = func_y(ans_x) print('x : '+str(ans_x)) print('y : '+str(ans_y)) # 以下プロット plt.plot(x, y, color='b', label='function') plt.plot(ans_x, ans_y, 'o', color='r', label='solution') plt.legend() plt.show() def func_y(x_dat): func = 2*x_dat**2 - 3*x_dat - 4 return func if __name__ == '__main__': main()
結果
x : 2.3507810593569047 y : -8.37196978409338e-12
もちろん、ほとんど同じです。
1つ目の方はfor文をどうにかできないか、悩みどころです。
時間があれば、詳しく考察したいですね。
Newton法
テイラー展開を利用した方法。
微小区間の傾きを利用するので、イメージとしては、高校のときに学んだ微分の考え方と同じです。
詳しい説明はこちら
必要なパラメータは初期値(解の近傍がよい)、微小区間の幅です。
ちなみに、Newton-Raphson法とも言うらしいです。
今回も、2x2-3x-4=0を解いてみます。
以下は、プログラム
# -*- coding: utf-8 -*- # import numpy as np import matplotlib.pyplot as plt def main(): init_val = 4.0 # 初期値 delta_val = 0.0001 # 微小区間の幅 min_x = 2 max_x = 5 x = np.linspace(min_x, max_x, 100) y = func_y(x) # 以下Newton法 for i in range(50): # 回数制限をかける f = func_y(init_val) df = (func_y(init_val + delta_val) - func_y(init_val)) / delta_val # 微分値の取得 init_val = init_val - f/df # 解の更新 ans_x = init_val # 近似解 ans_y = func_y(ans_x) # そのときのyの値(確認用) print('x : '+str(ans_x)) print('y : '+str(ans_y)) # 以下プロット plt.plot(x, y, color='b', label='function') plt.plot(ans_x, ans_y, 'o', color='r', label='solution') plt.legend() plt.show() def func_y(x_dat): func = 2*(x_dat**2) - 3*x_dat - 4 return func if __name__ == '__main__': main()
初期値が-2.0のときの結果
x : -0.8507810593582121 y : 0.0
初期値が4.0のときの結果
x : 2.350781059358212 y : 0.0
初期値によって、どこの解に収束するか決まります。
1つの初期値では、1つの解しか得られません。
また、今回は回数制限のみですが、一般的には終了条件を入れます。
二分法と同様にScipyにNewton法のモジュールがあります。
以下はプログラム
# -*- coding: utf-8 -*- # import numpy as np import matplotlib.pyplot as plt from scipy.optimize import newton def main(): init_val = 4.0 # 初期値 min_x = -5 max_x = 5 x = np.linspace(min_x, max_x, 100) y = func_y(x) ans_x = newton(func_y, init_val) # Newton法 ans_y = func_y(ans_x) print('x : '+str(ans_x)) print('y : '+str(ans_y)) plt.plot(x, y, color='b', label='function') plt.plot(ans_x, ans_y, 'o', color='r', label='solution') plt.legend() plt.show() def func_y(x_dat): func = 2*(x_dat**2) - 3*x_dat - 4 return func if __name__ == '__main__': main()
初期値が-2.0のときの結果
x : -0.8507810593582122 y : 8.881784197001252e-16
初期値が4.0のときの結果
x : 2.3507810593582126 y : 1.7763568394002505e-15
以上2つの手法以外にも非線形方程式を解く方法がありますが、気が向いたらします。
まだ、記事を書くことに慣れていないので、数式の入力やコードの貼り方がうまく出来ていないですね。
次回は、連立方程式の解法についてを予定しています。