함수의 근 탐색
이 단원에서는 임의의 1 변수 함수의 근을 찾는 함수들을 기술합니다. 라이브러리에서는 다양한 반복 풀이와 수렴 테스트를 위한 저수준의 기능들을 제공합니다. 이 기능들은 반복 단계를 완전히 제어할 수 있어, 잘 이용하면 사용자가 요구하는 적절한 해를 찾을 수 있습니다. 각각의 풀이 방법은 같은 작업 환경을 사용하기 때문에, 프로그램을 재컴파일 할 필요없이 여러 풀이 방법을 구동 중에 바꿀 수 있습니다. 각각의 풀이 인스턴스들은 스스로 상태를 추적하기 때문에, 다중 스레드 프로그램에서도 사용이 가능합니다.
이 단원에서 기술하는 함수들의 원형은 gsl_roots.h 에 정의되어 있습니다.
개요
1차원 근 탐색 알고리즘은 2가지로 나눌수 있습니다. 괄호법 1 과 기울기 연마 방법입니다.
괄호법의 경우 근이 존재하는 닫힌 구간에서 알고리즘이 실행되고, 반복적으로 구간의 크기를 줄여 허용치를 만족할 때까지 반복합니다. 이는 정밀한 오차 추정치를 제공해 줄 수 있습니다.
기울기 연마 방법 경우는 근에 대한 초기 추측을 개선해 나가면 진행합니다. 이러한 알고리즘의 경우는 근에 충분히 가까운 지점에서 시작한 경우에만 적절한 근을 찾을 수 있고, 속도를 위해서 엄격한 오차값의 측정을 희생합니다. 이 알고리즘은 근 근처에서 함수의 동작을 근사해, 더 높은 미분 차수의 개선점을 찾아 초기 조건을 개선합니다. 함수가 알고리즘과 잘 들어맞고, 좋은 초기 조건을 추측할 수 있을 때, 이러한 알고리즘은 빠르게 근으로 수렴합니다.
GSL에서는 이 두 가지 형태의 알고리즘을 비슷한 이름의 작업 공간으로 제공합니다. 사용자는 알고리즘을 위한 고수준의 작업공간을 제공하고, 라이브러리는 각 단계에서 필요한 개별 함수들을 제공합니다. 근 탐색 과정은 3가지로 나뉘어집니다.
알고리즘 \(T\) 대한, 근 풀이 공간 \(s\) 선언.
\(T\) 를 수행해 \(s\) 갱신.
\(s\) 의 수렴성 판별, 필요하면 계속 반복.
존재 구간이 주어진 알고리즘은 \(gsl_root_fsolver\) 작업공간을 줄 수 있고, 갱신 과정은 함수 자체만을 이용합니다(도함수를 사용하지 않습니다). 존재 구간이 주어지지 않은 경우는 \(gsl_root_fdfsolver\) 작업공간을 설정합니다. 이 경우 함수와 함수의 도함수( \(fdf\) 라 부릅니다 )가 갱신에 사용됩니다.
주의 사항
명심해야 할 점은 여기서 제공하는 근 탐색 함수가 1번에 1개의 근만을 찾는다는 사실입니다. 탐색 공간에 여러 근들이 존재한다면 가장 먼저 찾아진 근을 반환합니다. 하지만, 어떤 근이 반환될 지 예측하는 것은 매우 어렵습니다. 근이 여러개 있는 공간에서 한 근을 찾아도 대부분의 경우 오류는 발생하지 않습니다.
다중근을 가지는 함수를 다룰 때는 주의해야 합니다. 예를 들어서 \(f(x) = (x-x_0)^2\) 나 \(f(x) = (x-x_0)^3\) 등의 경우가 있습니다. 괄호법 알고리즘의 경우 다중근 중에서 짝수 중근을 가지는 함수에 사용할 수 없습니다. 알고리즘의 초기 구간은 반드시 0을 지나는 지점을 포함해, 구간의 한쪽은 양수를 다른쪽은 음수를 가져야 합니다. 하지만, 짝수번 중첩된 다중근 지점은 0을 지나지 않고 접하므로 초기 조건을 만족하지 못합니다. 홀수번 중첩된 근에 대해서는 잘 작동합니다(e.g. 3중근, 5중근…). 도함수를 이용하는 기울기 연마 방법은 일반적으로 큰 중첩을 가지는 다중근에서 잘 작동하지만, 수렴 속도가 감소합니다. 이 경우에, 스테퍼슨 방법 ( Steffensen’s method )을 사용해 다중근에 대해 수렴 속도를 늘릴 수 있습니다.
\(f\) 가 탐색 공간에서 반드시 근을 가질 필요는 없습니다만, 수치 근 탐색 함수는 근의 존재 여부를 판별 하기 위해 무분별하게 사용하는 것은 권장하지 않습니다. 근의 존재 여부 판별에는 더 좋은 방법들이 있습니다. 왜냐하면, 수치적 근 탐색은 일반적으로 실패할 가능성이 높기 때문입니다. 따라서, 함수에 대한 정보가 없을 경우, 바로 근 탐색으로 넘어가는 것은 좋은 생각이 아닙니다. 일반적으로 가장 좋은 방법은 근 탐색을 시작하기 전에 그래프를 그려서 시각적으로 확인을 해보는 것입니다.
또한, 함수에 따라 초기 추정 위치를 적절하게 정해주어야 합니다. 만일, 구간 오류가 발생한다면 그래프로 함수의 개형을 판별하고 적절하게 초기 추정치를 수정해 주어야 하는 경우도 있습니다(*).
풀이 시작하기
-
type gsl_root_fsolver
괄호법을 사용하는 근 탐색 알고리즘을 위한 풀이 공간입니다.
-
type gsl_root_fdfsolver
도함수가 필요한 근 탐색 알고리즘(기울기 연마 방법)을 위한 풀이 공간입니다.
-
gsl_root_fsolver *gsl_root_fsolver_alloc(const gsl_root_fsolver_type *T)
\(T\) 풀이 방법을 위한 새 풀이 공간을 할당해 포인터를 반환합니다. 예를 들어서 다음 코드는 이분법 방법을(Bisection method)을 위한 풀이 공간을 생성합니다.
const gsl_root_fsolver_type * T = gsl_root_fsolver_bisection; gsl_root_fsolver * s = gsl_root_fsolver_alloc (T);
만약, 이 공간을 할당하기에 메모리가 충분하지 않다면, 함수는 \(NULL\) 포인터를 반환하고 오류 관리자가 \(GSL_ENOMEM\) 코드의 오류를 보고합니다.
-
gsl_root_fdfsolver *gsl_root_fdfsolver_alloc(const gsl_root_fdfsolver_type *T)
미분 기반 풀이형 \(T\) 위한 새 풀이 공간을 할당해 포인터를 반환합니다. 예를 들어서 다음 코드는 뉴턴-렙슨 방법(Newton-Raphson method)를 위한 풀이 공간을 생성합니다.
const gsl_root_fdfsolver_type * T = gsl_root_fdfsolver_newton; gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc (T);
만약, 이 공간을 할당하기에 메모리가 충분하지 않다면, 함수는 \(NULL\) 포인터를 반환하고 오류 관리자가 \(GSL_ENOMEM\) 코드의 오류를 보고합니다.
-
int gsl_root_fsolver_set(gsl_root_fsolver *s, gsl_function *f, double x_lower, double x_upper)
함수 \(f\) 를 이용해 해를 찾는 풀이 공간 \(s\) 에 대해 초기 탐색 구간을 [ \(x_lower\) \(x_upper\) ] 로 잡고 풀이를 시작/재시작 시킵니다.
-
int gsl_root_fdfsolver_set(gsl_root_fdfsolver *s, gsl_function_fdf *fdf, double root)
함수 \(f\) 와 그 도함수 \(fdf\) 를 이용해 해를 찾는 풀이 공간 \(s\) 에 대해 초기 추정 위치를 \(root\) 로 잡고 풀이를 시작/재시작 시킵니다.
-
void gsl_root_fsolver_free(gsl_root_fsolver *s)
-
void gsl_root_fdfsolver_free(gsl_root_fdfsolver *s)
풀이 공간 \(s\) 할당된 메모리 영역을 해제합니다.
-
const char *gsl_root_fsolver_name(const gsl_root_fsolver *s)
-
const char *gsl_root_fdfsolver_name(const gsl_root_fdfsolver *s)
풀이 방법의 이름을 담은 포인터를 반환합니다. 예를 들어서,
printf("s is a '%s' solver\n", gsl_root_fsolver_name(s);
는 \(s is a 'bisection' solver\) 반환합니다.
근을 찾을 함수
근 탐색을 위해서는 1변수 연속함수를 제공해야하고, 알고리즘에 따라서는 1계 도함수도 제공해야 합니다.
다음과 같이 함수를 정의하면, 일반적인 함수의 계수처럼 함수를 쓸 수 있습니다.
-
type gsl_function
이 자료형은 계수를 갖는 일반적인 함수를 정의합니다.
double (* function) (double x, void * params)이 함수는 인자 \(x\) 계수 \(params\) 대해 \(f(x, params)\) 값을 반환하도록 정의해야합니다.
void * params함수의 계수를 나타내는 포인터입니다.
일반적인 2차 다항 함수를 이용해 예를 들어보도록 하겠습니다.
\(a =3, b=2, c=1\) 이면, 다음 코드와 같이 \(gsl_function\) 로 \(F\) 를 정의할 수 있습니다. 이 \(F\) 는 근 탐색 공간에 쓸 함수로 함수 포인터로 전달할 수 있습니다.
struct my_f_params { double a; double b; double c; };
double
my_f (double x, void * p)
{
struct my_f_params * params = (struct my_f_params *)p;
double a = (params->a);
double b = (params->b);
double c = (params->c);
return (a * x + b) * x + c;
}
gsl_function F;
struct my_f_params params = { 3.0, 2.0, 1.0 };
F.function = &my_f;
F.params = ¶ms;
함수 \(f(x)\) 는 매크로 \(GSL_FN_EVAL(&F, x)\) 를 이용해 평가해 볼 수 있습니다.
이는 gsl_math.h 에 정의되어 있습니다.
-
type gsl_function_fdf
이 자료형은 계수를 갖는 일반적인 함수와 이 함수의 1계 도함수를 정의합니다.
double (* f) (double x, void * params)이 함수는 인자 \(x\) 와 계수 \(params\) 에 대해 \(f(x, params)\) 값을 반환하도록 정의해야합니다.
double (* df) (double x, void * params)이 함수는 인자 \(x\) 와 계수 \(params\) 에 대해 \(f\) 1계 도함수 값; \(f'(x, params)\) 을 반환하도록 정의해야합니다.
void (* fdf) (double x, void * params, double * f, double * df)이 함수는 인자 \(x\) 와 계수 \(params\) 에 대해 함수 \(f\) \(f(x, params)\) 와 그 1계 도함수 \(df\) \(f'(x,params)\) 를 설정합니다. 이러한 방식은 독립된 함수 \(f(x, params)\) 와 \(f'(x,params)\) 를 제공함으로써 최적화를 시키는 방법으로, 도함수를 동시에 계산하는 방식에 비해 항상 빠릅니다.
void * params함수의 계수를 나타내는 포인터입니다.
다음은 \(f(x) = \text{exp}2x)\) 인 경우의 예시입니다.
double my_f (double x, void * params){
return exp (2 * x);
}
double my_df (double x, void * params){
return 2 * exp (2 * x);
}
void my_fdf (double x, void * params, double * f, double * df){
double t = exp (2 * x);
*f = t;
*df = 2 * t; /* uses existing value */
}
gsl_function_fdf FDF;
FDF.f = &my_f;
FDF.df = &my_df;
FDF.fdf = &my_fdf;
FDF.params = 0;
함수 \(f(x)\) 는 매크로 \(GSL_FN_FDF_EVAL_F(&FDF,x)\) 를 이용해서, 도함수 \(f'(x)\) 는 \(GSL_FN_FDF_EVAL_DF(&FDF,x)\) 이용해 평가해 볼 수 있습니다. 함수 \(y =f(x)\) 와 도함수 \(dy = f'(x)\) 는 매크로 \(GSL_FN_FDF_EVAL_F_DF(&FDF,x,y,dy)\) 를 써서 동시에 평가해 볼 수도 있습니다. 이 매크로는 \(f(x)\) 를 \(y\) 인자에 그리고 \(f'(x)\) 를 \(dy\) 인자에 저장합니다. 이 둘은 반드시 \(double\) 포인터여야 합니다.
탐색 경계와 추측 값
라이브러리에 있는 두 종류의 알고리즘 모두 탐색 경계나 추측 값을 필요로 합니다. 이 단락에서는 이러한 탐색 경계와 초기 값의 작동여부 그리고 함수 인자들이 어떻게 이들을 제어하는지 알아봅시다.
추측 값은 단순히 \(x\) 값을 의미합니다. 이 값은 원하는 근 정밀도가 될 때까지 알고리즘 내에서 반복됩니다. 이 값은 \(double\) 자료형을 가집니다.
탐색 경계는 구간의 끝 지점들을 의미합니다. 이는 구간의 길이가 요구되는 정밀도보다 작을 때까지 반복됩니다. 이 구간은 두 개 값, 하한 값과 상한 값으로 정의됩니다. 경계 값들을 간격에 포함할지, 말지 여부는 구간이 사용되는 풀이 흐름에 따라 달라집니다.
반복
다음 함수들은 각각의 알고리즘을 되풀이합니다. 각각의 함수는 상응하는 풀이공간의 상태를 1번의 반복 단계마다 갱신합니다. 동일한 기능이 모든 풀이에 작동하므로 코드를 수정하지 않고 실행 단계에서 다른 방법을 대체할 수 있습니다.
-
int gsl_root_fsolver_iterate(gsl_root_fsolver *s)
-
int gsl_root_fdfsolver_iterate(gsl_root_fdfsolver *s)
풀이공간 \(s\) 대한 단일 반복을 수행합니다. 만약 반복 과정에서 예상치 못한 문제가 생기면, 다음의 오류 값가 반환됩니다.
GSL_EBADFUNC이는 반복과정에서 함수나 도함수가 \(Inf\) 나 \(NaN\) 값이 되는 특이점에 도달했다는 뜻입니다.
GSL_EZERODIV함수의 미분값이 반복 지점에서 소멸해, 알고리즘이 0으로 나누는 것을 막았다는 뜻입니다.
이들은 탐색과정에서 언제나 현재 가장 좋은 근사 값을 유지합니다. 존재 구간이 주어진 경우 이에 더해, 근이 존재하는 가장 좋은 구간 값을 추적합니다. 이 정보는 다음과 같은, 보조 함수들을 이용해 접근할 수 있습니다.
-
double gsl_root_fsolver_root(const gsl_root_fsolver *s)
-
double gsl_root_fdfsolver_root(const gsl_root_fdfsolver *s)
호출시점에 풀이 공간 \(s\) 현재 저장된 근의 근사 값을 반환합니다.
-
double gsl_root_fsolver_x_lower(const gsl_root_fsolver *s)
-
double gsl_root_fsolver_x_upper(const gsl_root_fsolver *s)
호출시점에 풀이 공간 \(s\) 현재 탐색 구간의 하한, 상한 값을 반환합니다.
탐색 정지 인자들
근 탐색 과정은 다음 조건들 중 한가지가 충족되면 정지합니다.
사용자가 정의한 정확도를 만족하는 범주에서 근이 찾아진 경우.
사용자가 정의한 최대 반복 횟수에 도달한 경우.
오류가 발생한 경우.
이 조건들은 사용자가 제어해 볼 수 있습니다. 아래의 함수들은 사용자가 현재 결과를 여러 표준 방법들을 사용해, 정확도를 검증해 볼 수 있게 해줍니다.
-
int gsl_root_test_interval(double x_lower, double x_upper, double epsabs, double epsrel)
구간 [ \(x_lower\) \(x_upper\) ] 의 수렴을 절대 오차 \(epsabs\) 와 상대 오차 \(epsrel\) 를 이용해 검증합니다. 만약, 다음 조건을 만족하면, 수렴한다고 보고 \(GSL_SUCCESS\) 를 반환합니다.
\[|a-b| < epsabs +epsrel \cdot \text{min}(|a|, |b|)\]이 조건은 구간 \(x = [a,b]\) 가 원점을 포함하지 않을 때, 사용됩니다. 만약 구간이 원점을 포함한다면, \(\text{min}(|a|,|b|)\) 는 0으로 대체됩니다. (구간의 \(\|x\|\) 값 중에서 0이 자동으로 최소값이 됩니다. ) 이런 방법은, 원점에 가까운 근에 대해 상대 오차가 정확하게 추정할 수 있습니다.
구간에서 이 조건은, 구간에서 측정된 근 \(r\) 이 실제 근 \(r^*\) 에 대해 같은 조건을 만족시키는 것을 의미합니다.
\[|r - r^{*}| < epsabs + epsrel \cdot r^*\]이때, 실제 근 \(r^*\) 이 구간에 포함되어 있어야 합니다.
-
int gsl_root_test_delta(double x1, double x0, double epsabs, double epsrel)
수열 \(x0\) \(x1\) 수렴을 절대 오차 \(epsabs\) 와 상대 오차 \(epsrel\) 를 이용해 검증합니다. 이 판정은 다음 조건을 충족하면, \(GSL_SUCCESS\) 를 반환합니다.
\[|x_1 - x_0| < epsabs + epsrel \cdot |x_1|\]다른 경우에는 \(GSL_CONTINUE\) 를 반환합니다.
-
int gsl_root_test_residual(double f, double epsabs)
잔존값(residual value) \(f\) 에 대해 검증합니다. 다음 조건이 충족되면 함수는 \(GSL_SUCCESS\) 를 반환합니다.
\[|f| < epsabs``\]그리고 다른 경우에는 \(GSL_CONTINUE\) 를 반환합니다. 이 판정은 \(|f(x)|\) 가 충분히 작고, \(x\) 의 정확한 값이 중요하지 않은 상황에서 사용하는 것이 적절합니다.
괄호법 알고리즘
이 단락에서 설명할 괄호법을 이용한 근 탐색 알고리즘은 근이 포함되어 있음을 보장하는 초기 구간 설정이 필요합니다. \(a\) 와 \(b\) 가 구간의 양 끝 값이라면, \(f(a)\) 의 부호는 반드시 \(f(b)\) 와 달라야합니다. 이는 구간에서 함수 값이 0인 지점을 가로지르는 것을 보장합니다. 만약, 적절한 초기 구간이 제공되고 주어진 함수가 잘 작동한다면, 이 알고리즘들은 적절한 값을 항상 제공할 수 있습니다.
알아두어야 할 점은 짝수번 중첩된 다중근은 괄호법을 통해 찾을 수 없다는 점입니다. 이 경우 근이 \(x\) 축을 지나지 않고 접하기 때문입니다.
-
type gsl_root_fsolver_type
-
gsl_root_fsolver_type *gsl_root_fsolver_bisection
이분법 ( Bisection ) 알고리즘은 가장 간단한 괄호법 근 탐색 알고리즘입니다. 이 라이브러리에서 제공하는 함수중 가장 느린 선형 수렴 알고리즘이기도 합니다.
각각의 반복 단계에서, 주어진 구간이 이분되어 중간점에서의 함수 값이 계산됩니다. 이 값의 부호를 이용해 절반의 구간 중 어느 구간이 근을 포함하고 있지 않은지 판정합니다. 그러한 절반의 구간은 다음 단계에서 사라지고 근을 포함하는 새로운 작은 구간이 넘겨집니다. 이 과정은 간격이 충분히 작아질 때까지 무한정으로 반복될 수 있습니다.
각 단계에서, 근의 추정치는 현재 구간의 중간점으로 간주됩니다.
-
gsl_root_fsolver_type *gsl_root_fsolver_falsepos
할선법 ( false postion algorithm )은 선형 보간법에 기반해 근을 탐색합니다. 똑같이 선형으로 수렴하지만, 이분법보다 빠릅니다.
각 반복 단계에서, 구간의 끝지점의 함수 좌표 \((a, f(a))\) 와 \((b,f(b))\) 를 잇는 선이 그려지고, 그 선이 \(x\) 축과 만나는 지점을 중간점 으로 설정됩니다. 이 지점에서의 함수 값을 계산하고, 그 부호를 이용해 구간의 양 방향중 어느 부분 구간이 근을 가지지 않는지 판정합니다. 그러한 절반의 구간은 다음 단계에서 사라지고 근을 포함하는 새로운 작은 구간이 넘겨집니다. 이 과정은 간격이 충분히 작아질 때까지 무한정으로 반복될 수 있습니다.
근의 최적 추정치는 현재 반복에 대한 구간의 선형 보간값으로 간주됩니다.
-
gsl_root_fsolver_type *gsl_root_fsolver_brent
브렌트-데커 방법 ( Brent-Dekker algorithm )은 이분법에 보간법을 같이 적용한 방법입니다. 이 방법은 강력하면서도 빠른 알고리즘입니다. 이 책에서는 브렌트 방법 이라고 부를 것입니다.
각 반복 단계에서, 브렌트 방법은 보간 곡선을 사용해 함수를 근사합니다. 첫번째, 반복 단계에서 이는 2개의 끝 지점의 선형 보간입니다. 이후의 반복에서는 이전 단계의 3 지점을 역 2차 보간으로 근사해 더 놓은 정확도의 근사를 계산합니다. 보간 곡선과 \(x\) 축과의 교점은 그 단계에서의 근으로 추정되고 이 지점을 이용해 다음 단계에서 더 작은 구간을 생성합니다. 만약 이 교점이 구간을 벗어 났다면, 다시 1차 선형 보간으로 되돌아가 이를 반복합니다.
각 단계의 최적 추정치는 가장 최근 시행된 보간이나 이분법의 결과로 주어집니다.
-
gsl_root_fsolver_type *gsl_root_fsolver_bisection
기울기 연마 알고리즘
이 단락에서 설명할 기울기 연마 알고리즘은 근의 위치에 대한 초기 추정이 필요합니다. 이 방법은 절대적인 수렴 보장을 하지 않습니다. 주어지는 함수는 반드시 기울기 연마 방법으로 풀 수 있는 형태여야 하고, 초기 추정 값은 근에 충분히 가까워야 합니다. 이러한 조건이 갖추어졌을 때, 이 알고리즘은 2 차수의 수렴 속도를 가집니다.
이 알고리즘은 함수와 그 도함수를 필요로 합니다.
-
type gsl_root_fdfsolver_type
-
gsl_root_fdfsolver_type *gsl_root_fdfsolver_newton
이 방법은 뉴턴 방법(Newton’s Method) 으로 불리는 방법으로 표준적인 기울기 연마 알고리즘입니다. 이 알고리즘은 근에 대한 초기 추정으로 부터 시작합니다. 각각의 반복 단계에서, 현재 지점에서의 함수 \(f\) 의 접선이 그려집니다. 이 접선이 \(x\) 축과 만나는 지점이 새로운 추정값이 됩니다. 이러한 반복 과정은 다음과 같은 수열로 표현될 수 있습니다.
\[x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}\]뉴턴 방법은 단일 근에 대해 2차수의 수렴 속도를 가지며, 복수의 근에 대해서는 선형적인 수렴 속도를 가집니다.
-
gsl_root_fdfsolver_type *gsl_root_fdfsolver_secant
할선 방법(secant method) 은 뉴턴 방법과 흡사한 방법으로, 한가지 다른 점은 모든 단계에서 지점의 도함수의 값을 계산할 필요가 없습니다.
먼저, 첫번째 단계에서 뉴턴 방법으로 시작해 초기 추정값에서 다음 추정치를 얻습니다. .. math:
x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)}
이 다음 단계에서는 도함수의 값을 계산하지 않고 현재 지점과 그 이전 지점을 이용해 수치적으로 근사합니다.
\[f'_{est} = \frac{ f(x_i)-f(x_{i-1}) }{(x_i - x_{i-1})}\]따라서, \(i>0\) 일 때, 다음과 같습니다.
\[x_{i+1} =x_i - \frac{f(x_i)}{f'_{est}}\]도함수가 근 근처에서 변화가 크지 않을 경우 이러한 할선법은 속도 측면에서 도함수 값을 계산해야 하는 방법에 비해 빠른 속도를 제공해줍니다. 점근적으로 도함수의 값을 구하는 시간이 원래 함수 값을 구하는 시간의 0.44배 이상일 때, 이 방법이 뉴턴 방법보다 빠릅니다. 다른 수치 미분 방법이 그러하듯이 지점의 차이가 너무 작아지면 추정치가 취소되는 오류가 발생하기도 합니다.
단일근에서 이 방법은 \((1+\sqrt{5})/2\) 의 수렴 속도를 가집니다. 이는 대략 1.62 정도입니다. 다중근에 대해서는 선형적으로 수렴합니다.
-
gsl_root_fdfsolver_type *gsl_root_fdfsolver_steffenson
스테퍼슨 방법(Steffenson Method)**은 라이브러리의 함수 중 가장 빠른 수렴 속도를 제공합니다. 이 방법은 기본적인 뉴턴 알고리즘에 아티켄(Aitken)의 **델타-제곱(delta-suared) 가속을 추가한 방법입니다. 뉴턴 방법의 \(x_i\) 지점에서 가속 과정으로 새로운 수열 \(R_i\) 을 생성합니다.
\[R_i = x_i - \frac{(x_{i+1}-x_i)^2}{(x_{i+2}-2 x_{i+1}+x_i)}\]이 수열은 적절한 조건 하에서 본래 수열보다 빠르게 수렴합니다. 새로운 수열은 1개의 값을 생성하기 위해 3개의 항을 필요로 합니다. 따라서, 첫번째 반복 단계에서는 일반적인 뉴턴 방법의 값을 두번 째 이후의 반복에서 부터 가속된 값을 반환합니다. 만약, 가속 수열의 생성 과정에서 분모가 \(0\) 일 경우, 뉴턴 방법의 값을 이용하게 됩니다.
다른 가속 방법들과 마찬가지로 이 방법은 함수가 적절하지 않으면 불안정한 값을 내놓습니다.
-
gsl_root_fdfsolver_type *gsl_root_fdfsolver_newton
예제
어떤 근 탐색 알고리즘을 사용하건, 풀기 위한 함수를 준비해야 합니다.
예시로 이미 일반적인 2차 다항 함수를 살펴보았습니다.
먼저 함수 계수들을 정의하기 위해 demo_fn.h 헤더 파일이 필요합니다.
struct quadratic_params
{
double a, b, c;
};
double quadratic (double x, void *params);
double quadratic_deriv (double x, void *params);
void quadratic_fdf (double x, void *params,
double *y, double *dy);
함수의 정의는 분리된 다른 소스파일 \(demo_fn.c\) 에 있습니다.
double
quadratic (double x, void *params)
{
struct quadratic_params *p
= (struct quadratic_params *) params;
double a = p->a;
double b = p->b;
double c = p->c;
return (a * x + b) * x + c;
}
double
quadratic_deriv (double x, void *params)
{
struct quadratic_params *p
= (struct quadratic_params *) params;
double a = p->a;
double b = p->b;
return 2.0 * a * x + b;
}
void
quadratic_fdf (double x, void *params,
double *y, double *dy)
{
struct quadratic_params *p
= (struct quadratic_params *) params;
double a = p->a;
double b = p->b;
double c = p->c;
*y = (a * x + b) * x + c;
*dy = 2.0 * a * x + b;
}
첫 번째 프로그램은 \(gsl_root_fsolver_brent\) 를 사용해 브랜트 방법을 써서, 일반화된 2차 다항 함수를 이용해 다음과 같은 방정식을 풀 것입니다.
이 방정식의 해는 \(x= \sqrt{5} = 2.236068\) 입니다.
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include "demo_fn.h"
#include "demo_fn.c"
int
main (void)
{
int status;
int iter = 0, max_iter = 100;
const gsl_root_fsolver_type *T;
gsl_root_fsolver *s;
double r = 0, r_expected = sqrt (5.0);
double x_lo = 0.0, x_hi = 5.0;
gsl_function F;
struct quadratic_params params = {1.0, 0.0, -5.0};
F.function = &quadratic;
F.params = ¶ms;
T = gsl_root_fsolver_brent;
s = gsl_root_fsolver_alloc (T);
gsl_root_fsolver_set (s, &F, x_lo, x_hi);
printf ("using %s method\n",
gsl_root_fsolver_name (s));
printf ("%5s [%9s, %9s] %9s %10s %9s\n",
"iter", "lower", "upper", "root",
"err", "err(est)");
do
{
iter++;
status = gsl_root_fsolver_iterate (s);
r = gsl_root_fsolver_root (s);
x_lo = gsl_root_fsolver_x_lower (s);
x_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (x_lo, x_hi,
0, 0.001);
if (status == GSL_SUCCESS)
printf ("Converged:\n");
printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
iter, x_lo, x_hi,
r, r - r_expected,
x_hi - x_lo);
}
while (status == GSL_CONTINUE && iter < max_iter);
gsl_root_fsolver_free (s);
return status;
}
각각의 반복 시점에서 결과는 다음과 같습니다.
$./a.out
using brent method
iter [ lower, upper] root err err(est)
1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
Converged:
6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666
만약, 프로그램을 수정해서 브랜트 방법이 아닌 이분법 방법을 사용하면 ( \(gsl_root_fsolver_brent\) 를 \(gsl_root_fsolver_bisection\) 로 수정하면 됩니다), 이분법은 브랜트 방법보다 느리게 수렴하는 것을 관찰할 수 있습니다.
$./a.out
using bisection method
iter [ lower, upper] root err err(est)
1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000
2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000
3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000
4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000
5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500
6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250
7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625
8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312
9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656
10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828
11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414
Converged:
12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207
이번 프로그램은 똑같은 함수를 풀지만, 대신에 도함수 알고리즘을 사용해 볼 것입니다.
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include "demo_fn.h"
#include "demo_fn.c"
int
main (void)
{
int status;
int iter = 0, max_iter = 100;
const gsl_root_fdfsolver_type *T;
gsl_root_fdfsolver *s;
double x0, x = 5.0, r_expected = sqrt (5.0);
gsl_function_fdf FDF;
struct quadratic_params params = {1.0, 0.0, -5.0};
FDF.f = &quadratic;
FDF.df = &quadratic_deriv;
FDF.fdf = &quadratic_fdf;
FDF.params = ¶ms;
T = gsl_root_fdfsolver_newton;
s = gsl_root_fdfsolver_alloc (T);
gsl_root_fdfsolver_set (s, &FDF, x);
printf ("using %s method\n",
gsl_root_fdfsolver_name (s));
printf ("%-5s %10s %10s %10s\n",
"iter", "root", "err", "err(est)");
do
{
iter++;
status = gsl_root_fdfsolver_iterate (s);
x0 = x;
x = gsl_root_fdfsolver_root (s);
status = gsl_root_test_delta (x, x0, 0, 1e-3);
if (status == GSL_SUCCESS)
printf ("Converged:\n");
printf ("%5d %10.7f %+10.7f %10.7f\n",
iter, x, x - r_expected, x - x0);
}
while (status == GSL_CONTINUE && iter < max_iter);
gsl_root_fdfsolver_free (s);
return status;
}
다음은 위 프로그램에서 사용한 뉴턴 방법의 결과입니다.
$./a.out
using newton method
iter root err err(est)
1 3.0000000 +0.7639320 -2.0000000
2 2.3333333 +0.0972654 -0.6666667
3 2.2380952 +0.0020273 -0.0952381
Converged:
4 2.2360689 +0.0000009 -0.0020263
알아두면 좋은 점은, 현재 값과 이전 값과의 차이보다, 현재 값과 다음 값과의 차이를 이용해 오차를 더 정확하게 계산할 수 있다는 점입니다. 다른 도함수 풀이 방법은 \(gsl_root_fdfsolver_newton\) 를 \(gsl_root_fdfsolver_secant\) 나 \(gsl_root_fdfsolver_steffenson\) 로 교체해서 사용할 수 있습니다.
참고 문헌과 추가 자료
브렌트-데커 알고리즘 (Brent-Dekker algorithm)에 대한 추가 정보를 얻고 싶다면, 다음을 참고할 수 있습니다.
R. P. Brent, “An algorithm with guaranteed convergence for finding a zero of a function”, Computer Journal, 14 (1971) 422-425
J. C. P. Bus and T. J. Dekker, “Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function”, ACM Transactions of Mathematical Software, Vol.: 1 No.: 4 (1975) 330-345
- 1
근의 존재 구간이 주어진 경우(*)