[과학계산] Newton’s method, 뉴턴의 방법 python 코드 구현

이번 글에서는 Newton’s method (뉴턴의 방법)을 구현한 코드를 공개하도록 하겠다. Newton’s method 는 워낙 유명한 방법이니까 알고리즘 설명은 생략한다. 진짜 궁금하다면 위키에 있는 Newton’s method 에 대한 설명을 보도록 하자. 그러면 이제 본격적으로 Newton’s method (뉴턴의 방법)에 대해 잘 알아보도록 하겠습니다요.

Newton’s method (뉴턴의 방법) 구현

inverse matrix를 구해서 업데이트 하는 부분이 필요한데 나는 그것보다 gaussian eliminiation 을 이용해서 해를 구한 후 업데이트 하는 방식을 택했다. 이를 위해 gaussian elimination 방법을 택했다.

Gaussian elimination 구현

이제 이 Gaussian elimination 을 바탕으로 Newton’s method (뉴턴의 방법)을 아래와 같이 구현하였다.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import numpy as np
import matplotlib.pyplot as plt
# partially pivoting
def gaussian_elimination(A, b):
    #Ax=b 의 solution 구하기
    #입력변수 A,b
    n = len(b)
    aug_mat = np.hstack((A, b.reshape(11)))
    for i in range(n):
        pivot_row = np.argmax(np.abs(aug_mat[i:, i])) + i
        aug_mat[[i, pivot_row]] = aug_mat[[pivot_row, i]]
        pivot = aug_mat[i, i]
        aug_mat[i] /= pivot
        for j in range(i + 1, n):
            divided_factor = aug_mat[j, i] / aug_mat[i, i]
            aug_mat[j] = divided_factor * aug_mat[i]
    # back substitution
    x = np.zeros(n)
    for i in range(n  111):
        x[i] = aug_mat[i, 1 np.dot(aug_mat[i, i+1:n], x[i+1:])
    return x
cs

Newton’s method (뉴턴의 방법) 구현 코드

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# Newton’s method
def newton_method(F, J, initial_point, tolerance=1e8, max_iter=100):
    #Newton Method이용하기
    #F: 함수, J: F의 Jacobian matrix, initial_point: 초기 설정값, tolerance: 0이와 가까운지 판단하기 위함 임계치, max_iter: 최대 iteration 횟수
    x = initial_point
    history = [x.copy()]
    for iter in range(1, max_iter + 1):
        F_value = F(x)
        J_value = J(x)
        #partially pivoting
        delta_x = gaussian_elimination(J_value, F_value)
        x = x + delta_x
        history.append(x.copy())
        # convergence criterion
        if np.linalg.norm(F(x)) < tolerance:
            break
    return x, np.array(history)
cs

실제로 Newton’s method (뉴턴의 방법)을 사용하는 방법은 아래와 같이 함수를 정의하고 이함수의 gradient를 직접 입력해주면 된다.

Newton’s method (뉴턴의 방법) 예시

1
2
3
4
5
6
7
8
9
10
11
# 예제 1: F(x1, x2) = (x1^2 + x2^2 – 2, x1^2 – x2^2 – 1)
def func_for_3(x):
    return np.array([x[0]**2 + x[1]**2  2, x[0]**2  x[1]**2  1])
def jacob_3(x):
    return np.array([[2*x[0], 2*x[1]], [2*x[0], 2*x[1]]])
# 초기 추정치
initial_guess1 = np.array([1.50.5])
solution, history = newton_method(F, J, initial_guess)
cs

 

추천글

연립 1계 선형 미분방정식 풀이 가능한 경우, 풀이방법

Leave a Comment