16 - 2x^3 = 0 の解を求めてください。
簡単ですね、x = 2 が答えです。
x - y = 1, x + y = 3 なんてどうでしょう。
(x, y) = (2, 1)が答えです、これもすぐに求まりました。
では、5x - exp(x) = 0 の解はいくつでしょう。
答えは、x = 0.259171, 2.542641 です、ペンと紙があっても時間がかかりそうです。
このような複雑化された方程式の解を求める数値解法アルゴリズムには、
ニュートン法
二分法
Regula-falsi法
Bairstow法
などがあります。
今回は、そのうちのニュートン法をPythonで実装してみたいと思います。
目次
#1 ニュートン法とは
まず初めに、予想される真の解に近いと思われる値をひとつとる。次に、そこでグラフの接線を考え、その x 切片を計算する。このx切片の値は、予想される真の解により近いものとなるのが一般である。以後、この値に対してそこでグラフの接線を考え、同じ操作を繰り返していく。
https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95
要は適当にxの値を決めて、徐々に真の解へと値を近づけていくというものです。
プログラムの流れは以下のような感じで実装します。
初期値 x = x0 (ここで f(x0) ≠ 0)とします。
真の解を xtrue とし、xtrue = x0 + Δx を仮定すると
f(x0 + Δx) = 0 ‥‥‥①
が成り立ちます。
①をテイラー展開すると
f(xi + Δx) = 0 = f(xi) + f'(xi)Δx + ‥‥ ‥‥‥②
②式のΔxの一次項までを考慮してΔxの近似値を計算すると
(Δxの近似値) ≅ - f(xi) / f'(xi)
次に
xi += (Δxの近似値)
これを繰り返して、 (Δxの近似値) ≅ 0 となった時、
xi を真の解 xtrue としてプログラムを終了します。
#2 実装
from math import exp
def func(x):
return 5 * x - exp(x)
def diff_func(x):
return 5 - exp(x)
def main():
max = 10
dx = 0.0
x = float(input("初期値x0："))
for i in range(max):
dx = - func(x)/diff_func(x)
print("{}：dx = {:.6f}\n".format(i+1, dx))
x += dx
if -1.0e-6 < dx and dx < 1.0e-6:
print("x = {:.6f}\n".format(x))
exit()
print("error\n")
if __name__ == "__main__":
main()
>>出力1
初期値x0：1
1：dx = -1.000000
2：dx = 0.250000
3：dx = 0.009157
4：dx = 0.000015
5：dx = 0.000000
x = 0.259171
>>出力2
初期値x0：2
1：dx = 1.092877
2：dx = -0.385907
3：dx = -0.145130
4：dx = -0.018900
5：dx = -0.000298
6：dx = -0.000000
x = 2.542641
20数行とコンパクトに収めることができました。
特段難しい箇所があるわけでもなく、理解ができていれば実装は容易だと思います。
注意する点は、解が2つ存在するので初期値によって収束する値が異なるといったことくらいでしょうか。