Demo entry 6357129

Numerical Analysis midterm

   

Submitted by YSKang on Apr 21, 2017 at 07:40
Language: Python. Code size: 7.5 kB.

# -*- coding: utf-8 -*-
# 수치해석, 김상철 교수님
# 1. Incremental Search, Bisection, False Position 정리
# 20163179 홍승환

import numpy as np
import math
from scipy.optimize import fsolve

"""
    증분법을 계산하는 함수입니다.
    func : 대상 함수
    xmin, xmax : x의 최소, 최대값
    ns : 설정한 x 범위를 몇 개의 실제 값으로 쪼갤 것인지
"""
def incsearch(func, xmin, xmax, ns):
    # x 범위 설정
    x = np.linspace(xmin, xmax, ns)

    # 주어진 함수 실행 배열 세팅
    f = func(x)

    # 구간값의 갯수
    nb = 0

    # 구간값을 보관할 배열
    xb = []

    # 0부터 설정한 최고값까지 반복
    for k in np.arange(np.size(x) - 1):
        # 만약 f(x)랑 f(x+1)이 부호가 다르면 = 그래프가 x축을 뚫고 지나갔으면
        if np.sign(f[k]) != np.sign(f[k+1]):
            # 구간 하나 찾음!
            nb = nb + 1

            # 여기부터 여기까지임
            xb.append(x[k])
            xb.append(x[k + 1])

    # 구간 갯수(Number of Brackets), 구간 배열(Root Interval) 반환
    return nb, xb

"""
    이분법을 계산하는 함수입니다.
    func : 대상 함수
    xl : 계산할 x의 최저값 (lower guesses)
    xu : 계산할 x의 최고값 (upper guesses)
"""
def bisect(func, xl, xu):
    # 최대 반복 횟수, 기본 100으로 설정함
    maxit = 100

    # 예상된 오차의 범위, 기본 0.0001%로 설정
    es = 1.0e-4

    # f(최저값)과 x(최고값)을 곱함
    test = func(xl) * func(xu)

    """
        Bisection은, 계산할 x의 최저값과 최고값이 부호가 달라야 합니다.
        중간값 정의에 의해, 그 사이에 무조건적으로 근이 존재하기 때문입니다.
        그러나 f(최저값)과 f(최고값)을 곱한 값이 양수라면, 부호가 같다는 뜻이 됩니다.
        계산이 불가능하다는 뜻이므로, 함수를 종료합니다.
    """
    if test > 0:
        print("No sign change")
        return [], [], [], []
    
    # 반복 카운터
    iter = 0

    # 중심값 초기화
    xr = xl

    # 계산될 오차값 100%로 초기화
    ea = 100

    # 무한 반복
    while True:
        # 우선 현재 중심값을 xrold라는 변수에 저장해둠
        xrold = xr

        # 중심값을 현재의 최고, 최저값을 더한 후 2로 나눈 값(중앙값)으로 설정함
        xr = np.float((xl + xu) / 2)

        # 한 바퀴 더 돌았음
        iter = iter + 1

        # 중심값이 0이 아니라면
        if xr != 0:
            # 다음 식을 계산해서 오차값을 찾음
            # |{(현재 중심값) - (이전 중심값)} / (현재 중심값) | * 100
            ea = np.float(np.abs((np.float(xr) - np.float(xrold))/np.float(xr))*100)
        
        # f(최저값)과 x(최고값)을 곱함
        test = func(xl) * func(xr)

        # 만약 곱한 값이 양수라면
        if test > 0:
            # 최저값을 현재 중심값으로 설정함
            xl = xr
        # 아니고 만약 음수라면
        elif test < 0:
            # 최고값을 현재 중심값으로 설정함
            xu = xr
        # 만약 양수도, 음수도 아닌 0이라면
        else:
            # 오차 없음, 정확히 찾음
            ea = 0
        
        # 만약 (1) 계산된 오차값이 허용된 오차값보다 작다면 = 대강 이쯤되면 답이다 싶으면
        # 또는
        # 만약 (2) 설정한 최대 반복값보다 더 많이 돌았으면
        if np.int(ea < es) | np.int(iter >= maxit):
            # 계산 그만, 반복 종료
            break
    
    # 구한 중심값을 근으로 침
    root = xr

    # f(x)에 구한 근 넣어서 계산
    fx = func(xr)

    # 근, f(근) 값, 계산된 오차값, 반복값
    return root, fx, ea, iter

"""
    가위치법을 계산하는 함수입니다.
    func : 대상 함수
    xl : 계산할 x의 최저값 (lower guesses)
    xu : 계산할 x의 최고값 (upper guesses)
"""
def false_position(func, xl, xu):
    # 최고 반복 한계값, 에러 최댓값 설정
    maxit = 100
    es = 1.0e-4

    # f(최저값)과 x(최고값)을 곱함
    test = func(xl) * func(xu)

    """
        False Position은 Bisection과 마찬가지로 구간을 기반으로 판단하는 알고리즘입니다.
        그러므로 역시 중간값 정리를 적용하여, 양쪽 끝값의 부호가 달라야 계산이 가능합니다.
        그렇지 않다면, 계산을 하지 않고 함수를 종료합니다.
    """
    if test > 0:
        print("No sign change")
        return [], [], [], []
    
    # 반복 카운터 초기화
    iter = 0

    # 중간값 초기화
    xr = xl

    # 오차 초기화
    ea = 100

    # 무한 반복
    while True:
        # 이전 중간값 저장
        xrold = xr

        # False Position 알고리즘 : 삼각형의 닮은꼴 기반 계산법 (PPT 참조)
        xr = np.float(xu - func(xu) * (xl - xu) / (func(xl) - func(xu)))

        # 반복 카운터 1 증가
        iter = iter + 1

        # 중간값이 0이 아니면
        if xr != 0:
            # 다음 식을 계산해서 오차값을 찾음
            # |{(현재 중심값) - (이전 중심값)} / (현재 중심값) | * 100
            ea = np.float(np.abs((np.float(xr) - np.float(xrold))/np.float(xr)) * 100)
        
        # f(최저값)과 x(최고값)을 곱함
        test = func(xl) * func(xr)

        # 만약 곱한 값이 양수라면
        if test > 0:
            # 최저값을 현재 중심값으로 설정함
            xl = xr
        # 아니고 만약 음수라면
        elif test < 0:
            # 최고값을 현재 중심값으로 설정함
            xu = xr
        # 만약 양수도, 음수도 아닌 0이라면
        else:
            # 오차 없음, 정확히 찾음
            ea = 0
        
        # 만약 (1) 계산된 오차값이 허용된 오차값보다 작다면 = 대강 이쯤되면 답이다 싶으면
        # 또는
        # 만약 (2) 설정한 최대 반복값보다 더 많이 돌았으면
        if np.int(ea < es) | np.int(iter >= maxit):
            # 계산 그만, 반복 종료
            break

    # 찾은 중심값을 근으로 함
    root = xr

    # f(근) 계산
    fx = func(xr)

    # 근, f(근), 오차, 반복 횟수 반환
    return root, fx, ea, iter

# 실제적 계산값
print("[!] 0. Real Root\n")

fm = lambda m: np.sqrt(9.81 * m / 0.25) * np.tanh(np.sqrt(9.81 * 0.25 / m) * 4) - 36
m = fsolve(fm, 1)

print("[+] Real Root:", m)

# 증분법
print("[!] 1. Incremental Search\n")
g = 9.81; m = 68.1; cd = 0.25; v = 36; t = 4

fp = lambda mp: np.sqrt(g * np.asarray(mp) / cd) * np.tanh(np.sqrt(g * cd / np.asarray(mp)) * t) - v
nb, xb = incsearch(fp, 1, 200, 50)

# 구간의 갯수
print("[+] Number of Brackets =", nb)

# 구간들
print("[+] Root Interval =", xb)
print()

# 이분법
print("[!] 2. Bisection\n")
fm = lambda m: np.sqrt(9.81 * m / 0.25) * np.tanh(np.sqrt(9.81 * 0.25 / m) * 4) - 36
root, fx, ea, iter = bisect(fm, 40, 200)

# 근
print("[+] root:", root)

# f(근)
print("[+] f(root):", fx, "(Must Be Zero)")

# 계산된 오차
print("[+] Estimated Error:", ea, "(Must Be Zero Error)")

# 근을 찾기 위해 반복한 계산 횟수
print("[+] Iterated Number to Find Root: ", iter)
print()

# 가위치법
print("[!] 3. False Position\n")
fm = lambda m: np.sqrt(9.81 * m / 0.25) * np.tanh(np.sqrt(9.81 * 0.25 / m) * 4) - 36
root, fx, ea, iter = false_position(fm, 40, 200)

# 근
print("[+] root:", root, "(False Position)")

# f(근)
print("[+] f(root):", fx, "(Must Be Zero, False Position)")

# 계산된 오차
print("[+] Estimated Error:", ea, "(Must Be Zero Error, False Position)")

# 근을 찾기 위해 반복한 계산 횟수
print("[+] Iterated Number to Find Root: ", iter, "(False Position)")

# 기대 출력 결과
"""
[!] 0. Real Root

[+] Real Root: [ 142.73763311]

[!] 1. Incremental Search

[+] Number of Brackets = 1
[+] Root Interval = [139.08163265306123, 143.14285714285717]

[!] 2. Bisection

[+] root: 142.73765563964844
[+] f(root): 4.60891335763e-07 (Must Be Zero)
[+] Estimated Error: 5.3450468252827136e-05 (Must Be Zero Error)
[+] Iterated Number to Find Root:  21

[!] 3. False Position

[+] root: 142.73783844758196 (False Position)
[+] f(root): 4.20034974269e-06 (Must Be Zero, False Position)
[+] Estimated Error: 7.781013805708828e-05 (Must Be Zero Error, False Position)
[+] Iterated Number to Find Root:  29 (False Position)
"""

# 설명
"""
    먼저 실제 근값(0번)을 보면 142.73763311을 볼 수 있습니다. 
    이것이 실제 근이고, 이 값과 1, 2에서 구한 값을 비교합니다.

    1의 Root Interval을 보면, 139.08... 부터 143.14... 사이에 
    근이 있다는 것을 알 수 있습니다.
    증분법의 특성상 저 사이 쯤에 근이 있을 것이라는 추정 정도가 가능합니다.
    만약 linspace의 ns 값을 늘려서 대입값을 더 촘촘하게 한다면 
    더욱 정확한 구간 값을 도출할 수 있을 것입니다.

    2의 root를 보면, 142.73765563... 이 나옵니다. 실제 근과 소숫점 뒤
    4자리까지 일치하는 거의 정확한 값입니다.
    근을 x라고 할 때, 근이 맞다면 f(x)는 0이 되어야 합니다. 
    넣어보니 4.6089...e-07, 즉 0.000000460... 의 경미한 오차를 보입니다.
    또, 계산된 오차값을 보면, 위에서 기초적으로 설정한 0.0001%보다 작은 오차로
    21번만에 계산이 종료된 것을 알 수 있습니다.

    3의 root를 보면, 142.73783844... 이 나옵니다. 실제 근과 소숫점 뒤
    3자리까지 일치하는 것을 볼 수 있습니다.
    그러다 보니 f(근)과 계산된 오차 또한 실제값과 거리가 조금 더 있습니다.
    반복도 조금 더 한 것 같습니다.
    즉, 여기서는 이분법이 가위치법보다 더 정확한 값이 나왔다는 뜻이 됩니다.

    이건 원인이 있습니다.
    가위치법은, 근과 그래프의 설정된 양 끝값을 기준으로 삼각형 2개를 그린 후
    그 두 삼각형이 닮은꼴인가를 tan 함수와 비슷하게 판정하게 됩니다.
    그래서, 가위치법은 "원래 함수가 어떻게 생겼느냐"에 따라서 계산 횟수가 바뀌게 됩니다.
    하지만, 이분법은 삼각형과 상관이 없는 계산 방식을 가집니다.
    따라서 우리가 사용한 번지 점프 공식은 가위치법보다는 이분법이 더 적합한
    함수 모양이었다는 결론을 낼 수 있습니다.
"""

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).