무리수 √2의 값을 구해보자! : Part 5

안녕하십니까 Mr. K입니다

오랜만의 포스팅이군요 -_-;

지난주에 예고했던 것 처럼 √2에 대한 얘기를 하나 끄적이도록 하겠습니다

아, 그 전에

이번 포스트의 경우 도함수에 대한 얘기가 있으니, √2를 구하는 방법을 하나만 알아도 된다 하시는 분들은 패스하셔도 됩니다

Fixed Point Iteration

Most of the techniques we will discuss are iterative in nature and the first one is called the fixed point iteration scheme. Instead of looking for a zero of the function f(x), it determines a fixed point of an equivalent equation. The equation is rewritten in the form x = g(x), so when x is substituted in the function g(x) it returns the value x, hence the nomenclature fixed point. The converstion of the equation f(x) = 0 into one of the form x = g(x) is not always straightforward and definitely can be done in many ways. …

(위 내용은 AN INTRODUCTION TO Programming and Numerical Methods in MATLAB에서 일부 발췌했음을 알려드립니다)

간단히 얘기하면 다음과 같습니다

f(x) = 0 를 만족하는 x를 찾기 위해

주어진 방정식 f(x) = 0 를 g(x) = x 의 형태로 바꾸어서 고려하겠다는 말입니다

이 때 사용하는 것이 fixed point( 2차원 평면에서 y = x 직선 위에 있는 임의의 점 )입니다

이 iteration의 알고리즘을 간단히 보면 다음과 같습니다

x_0 : initial guess

x_{n+1} = g(x_n)

다시말하면, 초기값을 g(x)의 input으로 집어넣으면 나오는 output이 있는데,

그것을 다시 g(x)의 input으로 집어넣는 일을 반복하는 것입니다

그러면 이제,

우리의 관심사가 g(x)를 어떻게 찾을것인가? 로 바뀌었습니다

그런데 사실 g(x)를 만드는 작업 자체는 그다지 어렵지 않습니다

가장 간단히 만들 수 있는 g(x)는 x - f(x) 입니다

g(x) = x - f(x) = x   ⇔   - f(x) = 0   ⇔   f(x) = 0

이 외에도 f(x)의 모양에 따라 충분히 다양한 모양의 g(x)를 만들 수 있습니다

하나의 g(x)만으로 근을 찾을 수 있는 것이 아니기 때문에 같은 iteration이라도 g(x)는 다를 수 있습니다

이 때 g(x)를 이상한 모양으로 잡게 되면, 근이 존재하는데도 불구하고 이 iteration으로 찾을 수 없는 경우도 있습니다

그래서 g(x)를 조금 더 잘 잡기 위해 다음과 같은 정리를 사용합니다

함수 g와 초기 추정값 p0가 다음을 만족할 때,

만약 [a,b] 안의 모든 x에 대해 |g'(x)| < 1 이면 이 iteration으로 근을 찾을 수 있습니다

만약 [a,b] 안의 모든 x에 대해 |g'(x)| > 1 이면 이 iteration으로는 근을 찾을 수 없습니다

용어에 대해 잠깐 설명하자면,

i)에서

C[a,b]는 구간 [a,b]에서 연속인 함수들을 모아놓은 집합 입니다

iii)에서

A를 뒤집은 것 처럼 생긴 이 ∀ 기호는 for all 또는 arbitrary 의 뜻을 가지고 있습니다

정리가 다소 복잡하게 쓰여있는 느낌입니다만,

간단히 정리하면 다음과 같습니다

g와 g'가 연속함수이고, 초기 추정값을 p0라고 했을 때,

|g'(p0)| < 1 이면 이 iteration으로 근을 찾을 수 있고

|g'(p0)| > 1 이면 이 iteration으로 근을 찾을 수 없습니다

f(x) = x² - 2, 초기 추정값 = 1 이라고 하고 두가지 예를 보시겠습니다

1. g(x) = x - f(x) = x - (x² - 2)

이 경우에 g'(x) = 1 - 2x 이므로 g'(√2) < -1 이지요

(g'(x)에 대한 계산은 보통 손으로 하는 작업이므로 이 때는 √2가 얼마정도인지 짐작하고 있다고 가정합니다)

|g'(√2)| > 1 이므로 위 정리에 의해 근을 찾을 수 없지만, 직접 눈으로 보면 아래와 같습니다

파란색 그래프가 f(x), 붉은색 그래프가 g(x), 검은색 그래프가 직선 y=x 입니다

2. g(x) = 1 + x - x²/2

이 경우에는 g'(x) = 1 - x 이므로 g'(√2) ≒ -0.414 입니다

|g'(√2)| < 1 이므로 위 정리에 의해 근을 찾을 수 있고, 직접 눈으로 보면 아래와 같습니다

파란색 그래프가 f(x), 붉은색 그래프가 g(x), 검은색 그래프가 직선 y=x 입니다

둘의 차이가 보이십니까?

|g'(x)| > 1 인 경우에는 iteration을 수행할 수록 점점 근과 멀어지고,

|g'(x)| < 1 인 경우에는 iteration을 수행할 수록 점점 근과 가까워집니다

그래서 이 방법으로 근을 찾고자 할 경우

주어진 f(x)에 대해 적당히 g(x)를 잡아주는 것이 중요합니다

이제 코드로 구현한 것을 보시겠습니다

#include <stdio.h> double f( double ); double g( double ); double absolute( double ); double squareRoot( double ); const double tolerance = 0.0000000001; double realNumber; int maximum_iteration = 50; void main() { double root; printf("양의 제곱근을 구하고 싶은 실수를 입력해줏메 : "); scanf("%lf", &realNumber); root = squareRoot( realNumber ); if( root == -1 ) printf("\n근을 찾을 수 없었습니다, g(x)를 바꿔주세요\n"); else printf("\n%.4f의 제곱근은 %.8f임\n", realNumber, root); } double f( double x ) { return x * x - realNumber; } double g( double x ) { return x - f(x) / 4; } double absolute( double x ) { if( x >= 0 ) return x; else return -x; } double squareRoot( double r ) { int i; double xbef; // before double xaft; // after xbef = r - 1 / r; xaft = g(xbef); for( i = 0; i < maximum_iteration; i++ ) { if( absolute(xaft-xbef) <= tolerance ) break; xbef = xaft; xaft = g(xbef); } if( i == maximum_iteration ) { if( absolute(xaft-xbef) <= tolerance ) return xaft; else return -1; } else { return xaft; } }

잘 보시면 maximum_iteration이라는 변수가 있습니다

위의 1번 예처럼 iteration을 수행할 수록 근과 멀어지는 경우

iteration을 무한히 반복하면 무한히 멀어지는데, 이걸 막아주기 위한 수단으로 반복 회수를 제한하는 방법을 사용합니다

또, 예전의 √2 포스트에서도 나왔던 tolerance가 있습니다

정해진 회수만큼 iteration을 수행하다가

두 x값의 차이가 허용오차보다도 작아지면 같은 값으로 간주하고 그 중 하나를 반환합니다

어떤 f(x)에 대해 이 코드를 이용하기 위해 g(x)를 잘 잡았는데도 불구하고 원하는 근이 나오지 않는다면

다음과 같은 조치를 취하시면 됩니다

1) maximum_iteration을 조금 더 크게 잡아주세요

1)을 했음에도 불구하고 답이 나오지 않습니다

⇒ 2) f(x)의 일부분이 삼각함수입니까? 그러면 g(x)를 바꿔주세요

삼각함수의 모양때문에 원하는 근 ±(k×pi) 꼴의 근이 나올 수도 있습니다 -_-;

이번 포스트는 여기까지입니다//

(포스팅 소요시간 : 약 80분)