snovaのブログ

主にプログラミングやデジタルコンテンツについて書きます。最近はPython, Flutter, VRに興味があります。

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

また、グラフは

f:id:snova301:20180929164503p:plain
fig01.png

ちなみに、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

f:id:snova301:20180929164507p:plain
fig02.png

もちろん、ほとんど同じです。
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

f:id:snova301:20180929164511p:plain
fig03_01.png

初期値が4.0のときの結果

x : 2.350781059358212
y : 0.0

f:id:snova301:20180929164515p:plain
fig03_02.png

初期値によって、どこの解に収束するか決まります。
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つの手法以外にも非線形方程式を解く方法がありますが、気が向いたらします。

まだ、記事を書くことに慣れていないので、数式の入力やコードの貼り方がうまく出来ていないですね。
次回は、連立方程式の解法についてを予定しています。

参考文献

Google Play and the Google Play logo are trademarks of Google LLC.