본문 바로가기
컴퓨터

[파이썬04]방정식의 해 찾기

by Oh 선생 2019. 12. 27.

4장. 방정식의 해 찾기

0. Introduction

방정식의 해를 어떻게 찾을 수 있을까요? 다항방정식으로 한정한다면, 아마도 우린 인수분해를 하거나 근의 공식을 사용하거나 포기하는(...) 방법을 택할 겁니다. 그런데 컴퓨터를 사용하면 조금 다른 방식으로 해를 찾아나갈 수 있습니다. 이번 장에서는 컴퓨터를 이용해 방정식의 해를 찾는 간단한 방법에 대해 알아봅니다.

1. 이분법(Bisection method)

이분법이란 방정식의 해가 있음직한 구간을 하나 찍어서(!) 해가 있는지를 일단 확인하고, 그 구간을 반으로 잘라 그 반으로 자른 구간에도 해가 있는지를 다시 확인해보는 방법입니다. 이걸 계속 반복하다보면, 구간을 계속 반으로 자르게 되므로 구간은 아주 작아질 것이고, 결국 진짜 해와 구간의 양 끝값의 차가 매우 작아질 겁니다. 그럴 때 양 구간의 끝값을 해로 채택하는 방법이지요. 뭔가 이상하죠? 이게 어떻게 해가 되나 싶죠? 사실 컴퓨터로 무언가 방정식을 푼다고 할 때에는 정확한 해를 찾는 경우는 드뭅니다. 그보다는 문제가 없을 정도로 오차가 작은 근사해를 찾는 걸 주된 목적으로 하지요. 만약 진짜 해와 근사해의 차이가 만큼만 난다면, 실생활 영역에서 큰 문제는 없을 것이고, 필요하다면 오차를 더 줄여버릴 수 있으니까요.

이분법 프로그램을 만들기 전에, 우선 이분법의 알고리즘을 살펴봅시다. 기본적인 이분법 알고리즘은 다음과 같습니다.

알고리즘의 1번에는 여러분이 잘 아는 사잇값정리가 사용되고 있습니다. 구간의 양 끝값의 곱이 0보다 작단 말은 부호가 다르단 얘기고, 그러면 ( f가 연속함수라는 가정 하에) 해가 하나 그 사이에 존재하겠죠. 그 다음은 구간을 계속 줄여나가는 겁니다. 간단하죠? 그러면 지금부턴 코드를 작성해봅시다.

'''
Bisection method
'''
from sympy import Symbol, sympify
x=Symbol('x')
y=Symbol('y')

maxiter=100
tol=0.001
expr=sympify(input('해를 찾으려는 방정식을 입력하세요:'))
a=float(input('해를 찾을 왼쪽 범위를 입력해주세요:'))
b=float(input('해를 찾을 오른쪽 범위를 입력해주세요:'))

value_a=expr.subs(x,a)
value_b=expr.subs(x,b)
if value_a==0:
    print('방정식',expr,'=0의 해는 ',value_a,'입니다.')
if value_b==0:
    print('방정식',expr,'=0의 해는 ',value_b,'입니다.')    

if value_a*value_b>0 :
    print('입력한 구간에는 해가 존재하지 않습니다.')
else:
    for i in range(1, maxiter):
        xm=(a+b)/2
        value_xm=expr.subs(x,xm)
        err=(b-a)/2
        if abs(value_xm)<tol and abs(err)<tol:
            sol=xm
            print(sol, i)
            break
        else:
            if value_xm*value_a>0:
                a=xm
                value_a=value_xm
            else:
                b=xm

실행 결과는 다음과 같습니다.

파이썬에서 제곱은 **으로 표시합니다. , 여기에서 x**2 2는 "x^2-2"를 입력한 겁니다. 그리고 구간을 [0,3]으로 준거죠. 그랬더니 해가 1.414... 가 나온 겁니다. 몇 가지 방정식과 구간을 넣어보며 해를 찾아봅시다.

이제 코드 설명을 합시다.

'''
Bisection method
'''
from sympy import Symbol, sympify
x=Symbol('x')
y=Symbol('y')

처음 세 줄은 각주입니다. 각주는 '''로 첫 줄과 마지막 줄을 나타내줘야 합니다. 각주에 포함된 내용은 프로그램 실행에 영향을 미치지 않지만, 나중에 다른 이들이 코드를 읽을 때 도움을 주는 역할을 합니다. '''는 여러 줄을 각주로 묶을 때 주로 사용하고, 한 줄만 각주로 묶을 땐 줄의 맨 앞에 #을 추가합니다.

from sympy import Symbol, sympifysympy라는 라이브러리에서 Symbolsympify라는 모듈을 가져오라는 명령입니다. sympy는 파이썬에서 대수계산, 즉 문자를 이용해 수학을 할 수 있도록 만들어주는 라이브러리입니다. Symbolx, y같은 문자를 문자로 사용하겠다는 모듈이고, sympify는 입력된 문자열을 수식으로 인식하도록 만드는 모듈입니다.

x=Symbol(‘x’)란 앞으로 x라는 문자(왼쪽)에는 대수적 기호로 쓰이는 x(오른쪽)를 할당하겠다는 말입니다. 원래 x에는 수, 문자열 같은 개체만 할당할 수 있었잖아요? 여기에선 대수적 기호로서 x를 할당한단 말입니다. 이제 지금부터 x를 그냥 문자 x로 사용할 수 있어요.

maxiter=100
tol=0.001
expr=sympify(input('해를 찾으려는 방정식을 입력하세요:'))
a=float(input('해를 찾을 왼쪽 범위를 입력해주세요:'))
b=float(input('해를 찾을 오른쪽 범위를 입력해주세요:'))
value_a=expr.subs(x,a)
value_b=expr.subs(x,b)

maxiter=100은 해를 찾기 위한 최대 반복횟수를 설정해주는 겁니다. 반복횟수가 너무 적으면 해가 오차범위 내로 들어가지 못할 겁니다. tol은 허용가능한 오차범위입니다. 여기에선 0.001로 설정해줬습니다. expr=sympify(input('해를 찾으려는 방정식을 입력하세요:')) 부분은 input문을 통해 입력받은 식(=문자열)sympify 문을 통해 파이썬이 수식으로 인식하게끔 하는 과정입니다. 그리고 그 식을 expr이라는 변수에 할당해준거죠. a, b는 해를 찾을 범위를 실수로 바꾸어 할당하는 과정입니다.

value_a=expr.subs(x,a)expr이라는 식에 subs메소드를 사용해 x라는 문자에 a의 값을 대입하라는 명령어입니다. 그리고 그 계산된 값을 value_a라는 변수에 할당하라는 말이죠. sympy에서 식에 특정한 수나 문자를 대입할 때에는 .subs 메소드를 사용합니다. 그 아래 4줄은 읽어보면 알거라서 따로 설명하지 않겠습니다.

if value_a*value_b>0 :
    print('입력한 구간에는 해가 존재하지 않습니다.')
else:
    for i in range(1, maxiter):
        xm=(a+b)/2
        value_xm=expr.subs(x,xm)
        err=(b-a)/2
        if abs(value_xm)<tol and abs(err)<tol:
            sol=xm
            print(sol, i)
            break
        else:
            if value_xm*value_a>0:
                a=xm
                value_a=value_xm
            else:
                b=xm

이 부분이 이분법의 하이라이트입니다. 우선 value_a*value_b0보다 큰지 확인합니다. 크다면 부호가 같으므로 그 구간엔 해가 없는거죠. 그래서 두 번째 줄과 같은 내용을 인쇄합니다. 그게 아니라면(else) for문을 이용해 최대반복횟수만큼 for문 내부의 내용을 반복합니다.

우선 xm=(a+b)/2는 입력된 두 구간의 중점을 찾는 겁니다. 그리고 value_xm=expr.subs(x,xm)에선 식에 중점 xm의 값을 대입하죠. err=(b-a)/2는 구간의 크기를 반으로 나눈 겁니다. 잠정적인 오차, 즉 해와 근사해의 차이는 최대 이 크기 만큼이겠지요?

if abs(value_xm)<tol and abs(err)<tol: 함수식에 xm을 넣은 값이 오차보다 작고(어차피 0에 가까워야 하겠죠?), err도 오차보다 작을 때, 즉 근사해가 충분히 실제 해에 가까워졌을 때를 나타내는 조건문입니다. 만약 이 조건을 만족한다면, sol이란 변수에 그때까지의 근사해, xm을 넣고, sol과 그때까지의 계산횟수 i를 인쇄하게 되어있습니다. for문 내부에 사용된 break는 반복을 멈추고 반복문 바깥으로 나가라는 말입니다.

만약 조건을 만족하지 않는다면(else), xm을 넣은 값과 a를 넣은 값을 비교해서, 0보다 크면 구간 [a,xm]에 해가 없다는 말이므로 다음 구간을 [xm,b]로 택해주기 위해 a=xm으로 바꿔넣고, 식에 a를 넣은 값을 식에 xm을 넣은 값으로 바꿔넣는 겁니다. 그리고 만약 그게 아니라면 구간을 [a,xm]으로 바꿔주기 위해 b=xm으로 바꾸는 거지요.

코드가 이해되나요? 다음 연습문제를 풀어봅시다.

연습문제

(1) 방정식 x^3 -3x+1=0의 해를 구해보자.

(2) 이분법의 한계에 대해 생각해보자.

 

2. 시컨트법(Secant method)

secant는 우리에게 삼각함수로 익숙하지요. 그런데 기하적 의미로는 곡선을 끊는 할선이라는 의미도 갖고 있습니다. 시컨트법은 할선을 이용해 해를 찾아가는 방법입니다.

이러한 방법은 이분법에 비해 어떤 점이 더 좋을까요? 이 부분은 스스로 추측해보길 바랍니다. 알고리즘은 다음과 같습니다.

, 이제 코딩해봅시다.

'''
Secant method
'''
from sympy import Symbol, sympify
x=Symbol('x')
y=Symbol('y')

maxiter=100
tol=0.0001
expr=sympify(input('해를 찾으려는 방정식을 입력하세요:'))
x0=float(input('해를 찾을 임의의 첫 번째 x좌표를 입력해주세요:'))
x1=float(input('해를 찾을 임의의 두 번째 x좌표를 입력해주세요:'))

for i in range(1, maxiter):
    f0=expr.subs(x,x0)
    f1=expr.subs(x,x1)
    temp_x = x1-((x1-x0)/(f1-f0))*f1
    err=abs(x1-x0)
    if err < tol:
        print(temp_x, i)
        break
    else:
        x0=x1
        x1=temp_x

실행 결과는 다음과 같습니다.

코드를 설명합시다. 기본적인 틀은 이분법 때와 같으므로, 핵심적인 부분만 설명하겠습니다.

for i in range(1, maxiter):
    f0=expr.subs(x,x0)
    f1=expr.subs(x,x1)
    temp_x = x1-((x1-x0)/(f1-f0))*f1
    err=abs(x1-x0)
    if err < tol:
        print(temp_x, i)
        break
    else:
        x0=x1
        x1=temp_x

 

3. 뉴턴법(Newton method)

일반적으로 Secant법은 이분법보다 수렴속도가 빠르고, 시작점을 택하는데 자유롭다는 이점이 있지만 분모값이 작아져 계산이 안되는 경우가 가끔 있습니다. 그럴 땐 어떻게 하면 좋을까요? 편의상 f가 미분가능하다하면,

앞선 시컨트법에서 할선을 사용했다면, 뉴턴법에서는 접선을 사용한다는 겁니다. 그 외에 논리는 똑같아요. 하지만 접선을 사용하는 것만으로, 잡아야 할 점이 하나로 줄어들지요. , 코딩해봅시다.

'''
Newton method
'''
from sympy import Symbol, sympify, diff
x=Symbol('x')
y=Symbol('y')

maxiter=100
tol=0.0001
expr=sympify(input('해를 찾으려는 방정식을 입력하세요:'))
dexpr=diff(expr)
x0=float(input('초기 추측값 x를 입력해주세요:'))

for i in range(1, maxiter):
    f=expr.subs(x,x0)
    df=dexpr.subs(x,x0)
    temp_x = x0 - f/df
    err=abs(temp_x - x0)
    if err < tol:
        print(temp_x, i)
        break
    else:
        x0=temp_x

결과는 다음과 같습니다.

초깃값을 4로 대충 넣었더니, 5번만에 1.414...라는 근사해를 찾아줬지요. 코드를 이해해봅시다.

from sympy import Symbol, sympify, diff
x=Symbol('x')
y=Symbol('y')
 
maxiter=100
tol=0.0001
expr=sympify(input('해를 찾으려는 방정식을 입력하세요:'))
dexpr=diff(expr)

앞과 다르게 이번에는 diff 모듈을 불러옵니다. diff는 이름에서 알 수 있듯, 주어진 식의 미분을 구해주는 모듈입니다. 마지막 줄을 보면, dexpr=diif(expr)이라 되어 있지요? expr에 입력된 다항식을 미분하여 dexpr이라는 변수에 할당하라는 명령입니다.

for i in range(1, maxiter):
    f=expr.subs(x,x0)
    df=dexpr.subs(x,x0)
    temp_x = x0 - f/df
    err=abs(temp_x - x0)
    if err < tol:
        print(temp_x, i)
        break
    else:
        x0=temp_x

연습문제

x^3-2x-1=0의 해를 시컨트법과 뉴턴법을 이용해서 찾아보자.

 

이 밖에도 고정점 반복법 등이 있습니다만, 그런 방법은 후에 다른 기회를 통해 배우도록 하세요. 수고했습니다.

 

참고문헌

정상권, 김영희, 최홍원 (2012). 수치해석: 과학계산의 수학. 교우사. 

 

뉴턴. 뉴턴도 주식만은 폭망했다...

댓글0