ニュートン法

ニュートン法で平方根を求める。
　
・まずは√2　

def a(x):
    return (x**2 + 2) / (2 * x)

x = 2   # 初期値（0以外の値）
for i in range(6):
    x = a(x)
    print(x)
    
# 1.5
# 1.4166666666666667
# 1.4142156862745099
# 1.4142135623746899
# 1.414213562373095
# 1.414213562373095

　
ニュートン法は接線を利用することで解の近似を行っている。

　
続いて√3

def a(x):
    return (x**2 + 3) / (2 * x)

x = 2   # 初期値（0以外の値）
for i in range(8):
    x = a(x)
    print(x)
    
# 1.75
# 1.7321428571428572
# 1.7320508100147274
# 1.7320508075688772
# 1.7320508075688774
# 1.7320508075688772
# 1.7320508075688774
# 1.7320508075688772

　
誤差のせいか収束しないので、decimalを

from decimal import Decimal

def a(x):
    return Decimal(Decimal(x)**2 + s) / (2 * Decimal(x))

s = 3   # √s
x = 2   # 初期値（0以外の値）
for i in range(8):
    x = a(x)
    print(x)
    
# 1.75
# 1.732142857142857142857142857
# 1.732050810014727540500736377
# 1.732050807568877295254353946
# 1.732050807568877293527446342
# 1.732050807568877293527446342
# 1.732050807568877293527446342
# 1.732050807568877293527446342

収束おっけ
　
一般的に、x(new)=x-f(x)/f'(x)を繰り返すとf(x)=0の解の1つが得られる。

# y=f(x)
def f(x):
    return x**2 - 2

# 導関数
def diff(x):
    return (f(x + 0.001) - f(x)) / 0.001

# ニュートン法の漸化式
def a(x):
    return x - f(x) / diff(x)

x = 2  # 初期値
while f(x) > 10**-10:  # 誤差が10**-10になるまで繰り返す
    x = a(x)
    print(round(x, 10))

# 1.5001249688
# 1.4167014195
# 1.4142166238
# 1.4142135635
# 1.4142135624

　
解が複数あったり、あるいは解がなかったりのときは何か工夫しないと。　
　
matplotlibで描いたグラフのコード

def f(x):
    return x**2 - 2

# 導関数
def diff(x):
    return (f(x + 0.001) - f(x - 0.001)) / 0.002

# x=3 における接線の方程式
def tangent(x):
    return diff(3) * (x - 3) + f(3) 

x = np.linspace(-2, 4, 50)
y = f(x) 

plt.plot(x, y)
ax = plt.gca()
ax.spines['top'].set_position(("data",0))
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_color('none')
ax.xaxis.set_ticks_position("top") 

fn = "Times New Roman"

plt.plot([3, 3], [2, f(3)], "k:")
plt.plot([1, 4], [tangent(1), tangent(4)], color="r")
plt.text(2.9, -1, "$x_0$")
plt.text(1.8, -1, "$x_1$")
plt.title("$y=x^2-2$")

plt.xlim(-2, 4)
plt.show()