python으로 제곱근(root) 계산하는 함수 구현하기

1. Newton-Raphson Method

본문에서는 이 방법으로 제곱근을 구하는 함수를 작성하였고, 이 방법을 이해 해보도록 하자.

먼저 결론부터 말하면 아래와 같은 점화식에 따라 값을 계산하여 제곱근을 구한다.

$\sqrt{k} = x$라고 했을 때 $x_{n+1} = \frac{1}{2} \left(x_n + \frac{k}{x_n} \right)$, $x_0^2>k$

 

이제 위 점화식이 어떻게 유도되었는지를 알아보도록 하자.

 

점화식 유도

1. $\sqrt{k} = x$라고 했을 때, $x^2-k=0$라는 방정식을 세울 수 있다. 여기서 $f(x)=x^2-k$라고 하자.

2. $f(x_n)>0$이 되게 하는 점 $x_n$에 대해 $(x_n, f(x_n))$의 접선의 방정식을 구해보면 아래와 같다.
   - $g(x) = 2x_n (x-x_n) + f(x_n) = 2x_nx - x_n^2 - k$
3. 이 접선에서 $x$ 절편값을 구하면 아래와 같다.
   - $x_{n+1} = \frac{1}{2} \left(x_n + \frac{k}{x_n} \right)$

4. 이 $x_{n+1}$은 $x_n$보다 $\sqrt{k}$에 가까운 값이며 점화식을 반복해갈수록 제곱근값에 가까워진다.

    a) $k < x^2_n \rightarrow \frac{k}{x_n} < x_n  \rightarrow 2x_{n+1}=x_n+\frac{k}{x_n} < 2x_n$
    b) $\left( x_n - \frac{k}{x_n} \right) ^2 > 0 \rightarrow x_n^2 - 2k + \frac{k^2}{x_n^2} > 0 \rightarrow x_n^2 + 2k + \frac{k^2}{x_n^2} > 4k \rightarrow \left( x_n + \frac{k}{x_n} \right) ^2 > 4k \rightarrow 2x_{n+1} = x_n + \frac{k}{x_n} > 2\sqrt{k}$
    $\therefore \sqrt{k} < x_{n+1} < x_n$

 

 

2. Python 구현

 

이차방정식 f(x)

import numpy as np

k = 100

def f(x, k=k):
    return x**2 - k

예시로 $k=100$을 주었다.

 

접선의 방정식 g(x)

x0 = 13

def g(x, xn=x0, k=k):
    return 2*xn*x - xn**2 - k

 

 

X절편

def x_intercept(xn, k):
    return (xn + k/xn)/2

 

사실 이 함수만 있으면 된다. f, g 함수는 필요 없다.

 

 

제곱근 함수

def root(x):
    tol = 1e-10

    xn0 = x // 2
    xn1 = x_intercept(xn0, x)
    while abs(xn0 - xn1) > tol:
        xn0 = xn1
        xn1 = x_intercept(xn0, x)

    return xn1

$x_n$과 $x_{n+1}$의 차이가 `tol`이 되기 전까지 계속 반복하도록 작성하였다. 

 

 

root(k) # 10.0

root(k) == np.sqrt(k) # True
root(5) == np.sqrt(5) # True

`numpy` 내장 함수와도 값이 일치함을 확인할 수 있다.

 

속도도 전혀 느리지 않다.

대체로 10번 아래에서 찾는 것 같다.

10의 제곱근은 6번만에 답을 찾았다.

 

 

3. 전체 코드

github