【알고리즘】 25강. 수치해석 알고리즘

25강. 수치해석 알고리즘 모음 (출처 : KREYSZIG 공업수학 개정 10판)

 

추천글 : 【알고리즘】 알고리즘 목차 

---

1. Thomas 알고리즘 [본문]

2. Gauss 소거법 알고리즘 [본문]

3. Doolittle 법 [본문]

4. Cholesky 법 [본문]

5. Jacobi 반복법 [본문]

6. Gauss-Seidel 알고리즘 [본문]

7. 개선된 Euler 알고리즘 (Heum의 방법) [본문]

8. 고전적인 RUNGE-KUTTA 알고리즘 [본문]

9. 개선된 RUNGE-KUTTA 알고리즘 [본문]

10. RUNGE-KUTTA-FEHLBERG 알고리즘 [본문]

11. Adams-Bashforth 알고리즘 [본문]

12. Adams-Moulton 알고리즘 [본문]

13. Explicit Method [본문]

14. Fully (Simple) Implicit Method [본문]

15. Crank-Nicolson Method [본문]

16. 적분의 근사식 [본문]

17. 뉴턴-랩슨법 [본문]

18. Optimal Transport Theorem [본문]

19. Nelder-Mead method [본문]

---

 

1. Thomas 알고리즘 (A → A* = [a_jk] = [Ar]) [목차]

이 알고리즘은 Tridiagonal Matrix인 A가  유일해, x = [x_j]를 계산하거나 유일해가 존재하지 않음을 나타낸다. 

(단, A의 대각성분을 f, 대각아래 성분을 e, 대각위 성분을 g라고 하자.) 

 

INPUT: a_{j, n+1} = b_j인 n × (n+1) 첨가행렬 A* = [a_jk] 

OUTPUT: 식 (1)의 해 x = [x_j] 또는 식이 유일해를 갖지 않음을 알리는 메시지

// DECOMPOSITION

DO k = 2, ..., n

e_k = e_k  / f_k-1

f_k = f_k - e_k × g_k-1

END DO

// FORWARD SUBSTITUTION

DO k = 2, ..., n

r_k = r_k - e_k × r_k-1

END DO

// BACKWARD SUBSTITUTION

x_n = r_n / f_n

DO k = n - 1, ..., 1

 x_k = (r_k - g_k × x_k+1) / f_k

END DO

End GAUSS 

 

 

2. Gauss 소거법 알고리즘 (A → A* = [a_jk] = [Ab]) [목차]

이 알고리즘은 유일해, x = [x_j]를 계산하거나 유일해가 존재하지 않음을 나타낸다.

 

INPUT: a_{j, n+1} = b_j인 n × (n+1) 첨가행렬 A* = [a_jk]

OUTPUT: 식 (1)의 해 x = [x_j] 또는 식이 유일해를 갖지 않음을 알리는 메시지

For k = 1, ..., n - 1 do:

m = k

For j = k + 1, ..., n, do:

If |a_mk| < |a_jk|    then m = j    End

End

If a_mk = 0    then 

OUTPUT "유일해 없음"

Stop

End

Else k행을 m행과 교환

For j = k + 1, ..., n, do:

m_jk = a_jk / a_kk

For p = k + 1, ..., n + 1, do:

a_jp = a_jp - m_jka_kp

End

End

End

If a_nn = 0    then 

OUTPUT "유일해가 존재하지 않음"

Stop

End

x_n = a_n,n+1 / a_nn [역대입 시작] 

For i = n - 1, ..., 1, do:

x_i = (1/a_ii)(a_i,n+1 - (a_i,i+1x_i+1 + ... + a_{i, n}x_n))

End

OUTPUT x = [x_j]

Stop

End GAUSS 

 

 

3. Doolittle법 [목차]

Ax = LUx = b

Ux = L^-1b = y

x = U^-1y 

 

 

4. Cholesky법 [목차]

Ax = LL^Tx = b

L^Tx = L^-1b = y

x = (L^T)^-1y

 

 

5. Jacobi 반복법 [목차]

ex.

x₁ - 0.25x₂ - 0.25x₃ = 50            x₁ = 0.25x₂ + 0.25x₃ + 50

-0.25x₁ + x₂ - 0.25x₄ = 50          x₂ = 0.25x₁ + 0.25x₄ + 50

-0.25x₁ + x₃ - 0.25x₄ = 25          x₃ = 0.25x₁ + 0.25x₄ + 25

-0.25x₂ - 0.25x₃ + x₄ = 25          x₄ = 0.25x₂ + 0.25x₃ + 25

cf.

x^(m+1) = b - Lx^(m) - Ux^(m) = b - (L+U)x^(m)

 

 

6. Gauss-Seidel 알고리즘 (A, b, x⁽⁰⁾, ε, N) [목차]

이 알고리즘은 A = [a_jk]가 a_jj ≠ 0, j = 1, ..., n인 n × n 행렬일 때 주어진 초기 근사해 x⁽⁰⁾에 대해 Ax = b의 해 x를 계산한다.

 

INPUT: A, b, 초기 근사값 x⁽⁰⁾, 허용 오차 ε, 최대 반복수 N

OUTPUT: 근사해 x^(m) = [x_j^(m)] 또는 x^(N)이 허용오차 조건을 만족하지 않는다는 오류 메시지

For m = 0, ..., N - 1 do:

For j = 1, ..., n do:

x_j^(m) = (1/a_jj)(b_j - (a_j1x₁^(m+1) + ... + a_{j, j-1}x_j-1^(m+1)) - (a_{j, j+1}x_j+1^(m) + ... + a_jnx_n^(m)))

End

If  ^∀j, |x_j^(m+1) -  x_j^(m)|    then

OUTPUT x^(m+1)

Stop

End

End

OUTPUT

End GAUSS-SEIDEL

 

ex.

x₁ - 0.25x₂ - 0.25x₃ = 50            x₁ = 0.25x₂ + 0.25x₃ + 50

-0.25x₁ + x₂ - 0.25x₄ = 50          x₂ = 0.25x₁ + 0.25x₄ + 50

-0.25x₁ + x₃ - 0.25x₄ = 25          x₃ = 0.25x₁ + 0.25x₄ + 25

-0.25x₂ - 0.25x₃ + x₄ = 25          x₄ = 0.25x₂ + 0.25x₃ + 25

cf. 

x^(m+1) = b - Lx^(m+1) - Ux^(m)  ⇔  x^(m+1) = (I + L)^-1b - (I + L)^-1Ux^(m)

 

 

7. 개선된 Euler 알고리즘 (Heum의 방법) (f, x₀, y₀, h, N) [목차]

이 알고리즘은 초기값 문제 y' = f(x, y), y(x₀) = y₀의 해를 등간격 x₁ = x₀ + h, x₂ = x₀ + 2h, ..., x_N = x₀ + Nh에서 계산한다.

단, f는 [x₀, x_N] 구간에서 이 문제가 유일한 해를 갖는 함수이다.

 

INPUT: 초기값 x₀, y₀, 단계간격 h, 단계수 N

OUTPUT: x_n+1 = x₀ + (n+1)h에서 해 y(x_n+1)에 관한 근사화 y_n+1 (n=0, ..., N-1)

For n = 0, 1, ..., N-1 do:

x_n+1 = x_n + h

k₁ = hf(x_n+y_n)

k₂ = hf(x_n+1, y_n+k₁)

y_n+1 = y_n + 0.5(k₁ + k₂)

OUTPUT x_n+1, y_n+1

End

Stop

End EULER

 

 

8. 고전적인 RUNGE-KUTTA 알고리즘 (f, x₀, y₀, h, N) [목차]

이 알고리즘은 초기값 문제 y' = f(x, y), y(x₀) = y₀의 해를 등간격점

 

x₁ = x₀ + h, x₂ = x₀ + 2h, ..., x_N = x₀ + Nh

 

에서 계산한다. 여기서 f는 주어진 문제가 구간 [x₀, x_N]에서는 유일한 해를 갖도록 정해진다.

 

INPUT: 초기값 x₀, y₀, 단계간격 h, 단계수 N

OUTPUT: 단계점 x_n+1 = x₀ + (n+1)h에서 해 y(x_n+1)의 근사값 y_n+1 (n = 0, ..., N-1)

For n = 0, 1, ..., N-1 do:

k₁ = hf(x_n, y_n)

k₂ = hf(x_n + 0.5h, y_n + 0.5k₁)

k₃ = hf(x_n + 0.5h, y_n + 0.5k₂)

k₄ = hf(x_n + h, y_n + k₃)

x_n+1 = x_n + h

y_n+1 = y_n + (1/6)(k₁ + 2k₂ + 2k₃ + k₄)

OUTPUT x_n+1, y_n+1

End

Stop

End RUNGE-KUTTA

 

 

9. 개선된 RUNGE-KUTTA 알고리즘 (f, x₀, y₀, h, N) [목차]

이 알고리즘은 초기값 문제 y' = f(x, y), y(x₀) = y₀의 해를 등간격점

x₁ = x₀ + h, x₂ = x₀ + 2h, ..., x_N = x₀ + Nh

에서 계산한다. 여기서 f는 주어진 문제가 구간 [x₀, x_N]에서는 유일한 해를 갖도록 정해진다.

 

INPUT: 초기값 x₀, y₀, 단계간격 h, 단계수 N

OUTPUT: 단계점 x_n+1 = x₀ + (n+1)h에서 해 y(x_n+1)의 근사값 y_n+1 (n = 0, ..., N-1)

For n = 0, 1, ..., N-1 do:

k₁ = hf(x_n, y_n)

k₂ = hf(x_n + (1/4)h, y_n + (1/4)k₁)

k₃ = hf(x_n + (3/8)h, y_n + (3/32)k₁ + (9/32)k₂)

k₄ = hf(x_n + (12/13)h, y_n + (1932/2197)k₁ - (7200/2197)k₂ + (7296/2197)k₃)

k₅ = hf(x_n + h, y_n + (439/216)k₁ - 8k₂ + (3680/513)k₃ - (845/4104)k₄)

x_n+1 = x_n + h

y_n+1 = y_n + (25/216)k₁ + (1408/2565)k₃ + (2197/4104)k₄ + (-1/5)k₅

OUTPUT x_n+1, y_n+1

End

Stop

End RUNGE-KUTTA

 

 

10. RUNGE-KUTTA-FEHLBERG 알고리즘 (f, x₀, y₀, h, N) [목차]

이 알고리즘은 초기값 문제 y' = f(x, y), y(x₀) = y₀의 해를 등간격점

 

x₁ = x₀ + h, x₂ = x₀ + 2h, ..., x_N = x₀ + Nh

 

에서 계산한다. 여기서 f는 주어진 문제가 구간 [x₀, x_N]에서는 유일한 해를 갖도록 정해진다.

 

INPUT: 초기값 x₀, y₀, 단계간격 h, 단계수 N

OUTPUT: 단계점 x_n+1 = x₀ + (n+1)h에서 해 y(x_n+1)의 근사값 y_n+1 (n = 0, ..., N-1)

For n = 0, 1, ..., N-1 do:

k₁ = hf(x_n, y_n)

k₂ = hf(x_n + (1/4)h, y_n + (1/4)k₁)

k₃ = hf(x_n + (3/8)h, y_n + (3/32)k₁ + (9/32)k₂)

k₄ = hf(x_n + (12/13)h, y_n + (1932/2197)k₁ - (7200/2197)k₂ + (7296/2197)k₃)

k₅ = hf(x_n + h, y_n + (439/216)k₁ - 8k₂ + (3680/513)k₃ - (845/4104)k₄)

k₆ = hf(x_n + (1/2)h, y_n - (8/27)k₁ + 2k₂ + (3544/2565)k₃ - (1859/4104)k₄ - (11/40)k₅)

x_n+1 = x_n + h

y_n+1 = y_n + (16/135)k₁ + (6656/12825)k₃ + (28561/56430)k₄ + (-9/50)k₅ + (2/55)k₆

OUTPUT x_n+1, y_n+1

End

Stop

End RUNGE-KUTTA

 

 

11. Adams-Bashforth 알고리즘 (x_n, x_n-1, x_n-2, x_n-3, f_n, f_n-1, f_n-2, f_n-3, h, N) [목차]

이 알고리즘은 x₀을 포함한 어떤 열린 구간에서 유일한 해를 가지는 초기값 문제

 

y' = f(x, y),    y(x₀) = y₀,    f_n = (x_n, y_n)

 

를 고려한다. 이 방법은 피적분함수 f(x, y(x))를 Newton의 후진차분공식을 써서 보간다항식으로 대체하고 적분하는 것이다. 

 

INPUT: x_n, x_n-1, x_n-2, x_n-3, f_n, f_n-1, f_n-2, f_n-3, h, N

OUTPUT: 단계점 x_n+m = x₀ + (n+m)h에서 해 y(x_n+m)의 근사값 y_n+m (m = 1, ..., N)

k₄ = hf_n

k₃ = hf_n-1

k₂ = hf_n-2

k₁ = hf_n-3

For m = 1, ..., N do:

x_n+m = x_n+m-1 + h

y_n+m = y_n+(m-1) + (1/24)(55k₄ - 59k₃ + 37k₂ - 9k₁)

k₁ = k₂

k₂ = k₃

k₃ = k₄

k₄ = hf(x_n+m, y_n+m)

OUTPUT x_n+m, y_n+m

End

Stop

End ADAMS-BASHFORTH

 

 

12. Adams-Moulton 알고리즘 (x_n, x_n-1, x_n-2, x_n-3, f_n, f_n-1, f_n-2, f_n-3, h, N) [목차]

이 알고리즘은 x₀을 포함한 어떤 열린 구간에서 유일한 해를 가지는 초기값 문제

 

y' = f(x, y),    y(x₀) = y₀,    f_n = (x_n, y_n)

 

를 고려한다. 이 방법은 예측자-수정자 방법을 추가한 것이다. 오차추정은 다음과 같다.

 

ε_n+1 = (1/15)(y_n+1 - y*_n+1)

 

INPUT: x_n, x_n-1, x_n-2, x_n-3, f_n, f_n-1, f_n-2, f_n-3, h, N

OUTPUT: 단계점 x_n+m = x₀ + (n+m)h에서 해 y(x_n+m)의 근사값 y_n+m (m = 1, ..., N)

k₄ = hf_n

k₃ = hf_n-1

k₂ = hf_n-2

k₁ = hf_n-3

For m = 1, ..., N do:

x_n+m = x_n+m-1 + h

y*_n+m = y_n+(m-1) + (1/24)(55k₄ - 59k₃ + 37k₂ - 9k₁)

k* = f(x_n+m, y*_n+m)

y_n+m = y_n+(m-1) + (1/24)(9k* + 19k₄ - 5k₃ + k₂)

k₁ = k₂

k₂ = k₃

k₃ = k₄

k₄ = hf(x_n+m, y_n+m)

OUTPUT x_n+m, y_n+m

End

Stop

End ADAMS-MOULTON

 

 

13. Explicit Method [목차]

x^(m+1) = Ax^(m) + b

 

ex.

3×3: A = ((1 - 2λ, λ, 0)^T, (λ, 1 - 2λ, λ)^T, (0, 1 - 2λ, λ)^T), b = (λx₀, 0, λx_n+1)^T

 

 

14. Fully (Simple) Implicit Method [목차]

x^(m+1) = A^-1x^(m) - A^-1b

 

ex.

3×3: A = ((1 + 2λ, -λ, 0)^T, (-λ, 1 + 2λ, -λ)^T, (0, 1 + 2λ, -λ)^T), b = (λx₀, 0, λx_n+1)^T

 

 

15. Crank-Nicolson Method [목차]

x^(m+1) = A_im^-1(A_exx^(m) + b_ex) - A^-1b_im

 

ex.

3×3: A = ((1 + 2λ, -λ, 0)^T, (-λ, 1 + 2λ, -λ)^T, (0, 1 + 2λ, -λ)^T), b = (λx₀, 0, λx_n+1)^T

 

 

16. 적분의 근사식 [목차]

⑴ 심슨 공식(Simpson formula) 등

 

 

17. 뉴턴-랩슨법(Newton-Raphson method) [목차]

 

 

18. Optimal Transport Theorem [목차]

 

 

19. Nelder-Mead method (아메바 방법, 넬더-미드 방법, 활강단체법) [목차] 

이 알고리즘은 다차원 공간의 손실함수의 최솟값 또는 최댓값을 찾기 위한 수치 해석 방법이다.

 

import numpy as np
from scipy.optimize import minimize
from scipy.stats import spearmanr

def model(X, coeffs):
    return np.sum(X.T * coeffs, axis=1)

def objective_function(coeffs, X, target):
    prediction = model(X, coeffs)
    return -spearmanr(target, prediction)[0]

def optimize_spearman(X, target, initial_guess=None):
    if initial_guess is None:
        initial_guess = np.zeros(X.shape[0])  # Default to zero coefficients

    result = minimize(objective_function, initial_guess, args=(X, target), method='Nelder-Mead')
    optimal_coeffs = result.x
    maximized_corr = -objective_function(optimal_coeffs, X, target)
    return optimal_coeffs, maximized_corr

# Example usage
# new_bulk = np.array(...) # Replace with your 11x471 matrix
# target_vector = np.array(...) # Replace with your 471-D target vector

# optimal_coeffs, maximized_corr = optimize_spearman(new_bulk, target_vector)
# print("Optimal coefficients:", optimal_coeffs)
# print("Maximized Spearman Correlation:", maximized_corr)

 

⑴ 구별개념 1. SLSQP method

⑵ 구별개념 2. L-BFGS-B method

 

입력: 2016.12.11 02:10

수정: 2024.03.28 22:20