二分法とニュートン法を用いて以下の非線形方程式の数値解を求める。

sin(x)=0

二分法

区間$x=[a,b]$で$f(x)$が連続かつ$f(a)f(b)<0$の時、区間内で$f(x)=0$となる解を一つ求める方法。

アルゴリズムは

$f(a_0)f(b_0)<0$を満たす解を求めたい範囲で$a_0、b_0$に初期値を代入する。
$c=\frac{(a_i+b_i)}{2}$とし、$f(c)f(a_i)<0$ならば$b_{i+1}=c,a_{i+1}=a_i$、$f(c)f(a_i)>0$ならば$b_{i+1}=b_i,a_{i+1}=c$で更新する。
2の逐次計算を行い、$|f(c)|<ε_1$もしくは$|a_i-b_i|<ε_2$を満たせば収束したとして計算を終了する。

・ソースコード

$x=π$の解を求めるため、初期値$a_0=3.0$、$b_0=3.5$とした。

import math

def f(x):
    return math.sin(x)

EPS1 = 0.00001
EPS2 = 0.00001

a = 3.0
b = 3.5

while True:
    c = (a + b)/2
    if f(a)*f(c) < 0:
        b = c
    else:
        a = c
    if abs(f(c)) < EPS1 or abs(a - b) < EPS2:
        break

print("x = %f" % c)

・結果

>>> print("x = %f" % c)
x = 3.141602

ニュートン法

$x_i$における$f(x)$の接線を求め、接線と$x$軸との交点を次の$x_{i+1}$として更新し、逐次計算を行う方法。

$x_i$における$f(x)$の接線$g(x)$は

g(x)-f(x_i)=f^{'}(x_i)(x-x_i)

接線と$x$軸との交点を次の$x_{i+1}$として更新するため、$g(x)=0$として

0-f(x_i) = f^{'}(x_i)(x_{i+1}-x_i)\\
x_{i+1} = x_i-\frac{f(x_i)}{f^{'}(x_i)} ... ①

アルゴリズムは

$x_0$に初期値を与える。
①式に従って次の$x_{i+1}$を求める。
2を繰り返し、$|f(x)|<ε$となったら収束したとし終了する。

・ソースコード

今回は微分値には解析解($cos(x)$)を用いた。

import math

def f(x):
    return math.sin(x)

def df(x):
    return math.cos(x)

EPS1 = 0.0001

x = 3.5

while True:
    x -= f(x)/df(x)
    if abs(f(x)) < EPS1:
        break
        
print("x = %f" % x)

・結果

>>> print("x = %f" % x)
x = 3.141594