루트격리(root isolation)
루트 격리는 근이 여러 개인 다항식에서 정해진 개수의 근이 존재하는 구간을 구하는 것이다. (주로 근이 1개인 구간을 구한다.) 이를 이용해서 근 구하기 알고리즘 등을 사용한다.
1. 스트룸 정리
아래에 있는 데카르트 부호 법칙에 기반한 방법이나, 빈센트 정리에 기반한 방법보다는 느리고 덜 효율적이다. 중근을 한 개로 간주한다는 특징이 있어 덜 효율적이지만 꽤 쓰이는 듯.
스트룸 정리 (Strum's theorem)
루트 격리 (root isolation) 위 다항식의 실수 근을 구하는 알고리즘을 짜려면 어떻게 해야할까? 뉴턴 방법이나 이분 매칭 등을 써서 그 값을 구할 수 있지만 위와 같이 실근이 2개 이상인 경우 모든
tsyang.tistory.com
위 글에서 이미 작성했으니 참고.
2. 데카르트 부호 법칙
데카르트의 부호 법칙은 다항식에서 계수의 부호가 바뀐 횟수가 양의 실근의 개수(복소수 근을 포함한)를 나타낸다는 법칙이다. 중근도 다 따로따로 센다. 그니까 3중근이면 근이 3개가 있다고 치는 셈. 만약 어떤 다항식 $p(x)$에서 계수의 부호가 N번 바뀌었다면 해당 다항식의 양의 실근의 개수는 N-C(0 또는 짝수) 개다.
예를 들어,
$p(x) = x^3-7x+7$ 이라면 부호가 (+,-,+)로 두 번 바뀌었으므로 양의 실근이 2개 혹은 0개이다.
마찬가지로 $p(-x)$의 계수를 이용해 음의 실근의 개수도 구할 수 있다.
그렇다면 루트 격리는 어떻게 할까? 다항식을 x축으로 이동시켜서 구해주면 된다.
예를 들어,
$p(x) = x^3-7x+7$는 부호가 (+,-,+)로 두 번 바뀌므로 양의 실근이 2개 혹은 0개이다.
$p(x+2) = x^3 + 6x^2 + 5x + 1$으로 양의 실근이 0개이다. 즉 열린 구간 (0,2)에서 양의 실근이 2개 혹은 0개가 있다는 것을 의미한다.
데카르트의 부호 법칙을 이용한 방법은 효율적이지만 중근을 따로 세고, 허수근을 포함하기 때문에 대신 느린 스트룸 정리를 쓰기도 한다.
3. 빈센트 정리
빈센트 정리는 가장 효율적인 방법 중 하나인데, 중근이 없는 다항식 (square-free polynomial)의 양의 실근의 개수를 구한다. 어려워서 다루지는 않음
다항식의 근 구하는 법
아벨의 증명에 따르면 5차 이상의 다항식에서는 근의 공식 같은 게 없어서 현실적으로 근의 근사치만을 구할 수 있다. (그니까 다항식을 $(x-a)(x-b)(x-c)...$이런식으로 표현할 수 없다는 뜻). 아래의 방법은 고차 다항식의 근의 근사치를 구하는 방법을 나열했다.
만약 다항식의 모든 근을 구해야 한다면 근을 격리시키거나, 해당근을 r이라고 하면 다항식을 다시 $(x-r)$로 나누어 계산하는 방법이 있다. 후자의 경우 정확도가 점점 떨어진다.
기초적인 방법
1. 이분법
이분법의 경우, $f(a)f(b)<0$인 a, b를 알아야 한다. 이때 구간 (a, b)를 bracket이라고 말한다. 앞서 말한 루트 격리와 쓰기 좋다.
이분법은 그냥 이분 탐색 알고리즘이랑 똑같다. 구간을 반씩 쪼개며 정확도를 높여가는데 이런 걸 linear 하게 수렴한다 그러거나 α=1.0의 차수로 수렴한다고 하는 듯.
빠른 방법은 아니지만 어떤 경우에든 균일한 속도를 보이며 근으로의 수렴을 보장한다. (stable)
2. Regula Falsi (False position method, 가위치법)
한글로는 가위치법이라고 부르는 듯.
이분법과 마찬가지로 bracket이 필요하다. 식으로 설명하면 복잡한데 그림으로 보면 알기 쉽다.
가위치법도 항상 근으로 수렴함이 보장된다. (stable)
가위치법도 linear 하게 수렴한다고 하는데 오히려 이분법보다 느릴 때가 있다.
위 그림과 같은 상황에서 만약 이분법이었다면 매 step을 반복할 때마다 구간이 1/2씩 줄어들었을 것이다. 즉, 위와 같은 상황에서 가위치법은 이분법보다 느리다. 이때 오른쪽 위의 점과 같이 계속 선택되는 점을 정체점(stagnation point)라고 한다.
여기서 수정된 가위치법(modified regular falsi)이 등장하는데, 대부분 저 정체점을 적절하게 바꿔주는 방법을 포함한다.
가장 유명한 Illinois algorithm은, 만약 정체점이 생기면 정체점의 값에 1/2를 곱해준다. 이때, 정체점은 두 번 이상 선택된 점 이런 식으로 정의해줄 수 있다.
Illinois algorithm을 적용한 가위치법은 평균적으로 1.442의 차수를 지닌다. 즉 superlinear 하게 수렴한다. 그 외에도 양쪽 점의 값의 차이로 정체점의 값을 바꿔주는 방법도 있다. (Anderson–Björck algorithm)
+ 2021/9/2 추가) 실제로 가위치법을 돌려봤는데 1분동안 루프를 돌아서 강제종료했다. 확인해보니 x축과 만나는 점의 값이 매우 느리게 변환된다. 특히 뒤로갈수록 값의 변화가 점점 줄었다.
그래서 수정된 가위치법인 일리노이 알고리즘을 적용해봤더니 무려 16번의 반복만에 구해졌다.
간단한 테스트 & 코드
테스트
[Test]
public void GetARootInRange_Illinois_Test_ForNonSqaureRoot()
{
double targetAccuracy = 0.00001;
//평범한 근 (bracket)
Polynomial g = new Polynomial(-1, -3, 0, 0, 0, 1); // x^5 - 3x -1 (근 : -1.21465, -0.334734, 1.38879)
double ansg0 = g.GetARootInRange_Illinois(new DoubleRange(-10, -1));
double ansg1 = g.GetARootInRange_Illinois(new DoubleRange(-1, -0.334));
double ansg2 = g.GetARootInRange_Illinois(new DoubleRange(1, 1000));
Debug.Log($"{ansg0}, {ansg1}, {ansg2}");
Assert.AreEqual(g.GetValue(ansg0).DoubleEqual(0.0d, targetAccuracy), true);
Assert.AreEqual(g.GetValue(ansg1).DoubleEqual(0.0d, targetAccuracy), true);
Assert.AreEqual(g.GetValue(ansg2).DoubleEqual(0.0d, targetAccuracy), true);
}
코드
public double GetARootInRange_Illinois(in DoubleRange doubleRange,
double tolerance = 0.000001)
{
double a = doubleRange.min;
double b = doubleRange.max;
double fa = GetValue(a);
double fb = GetValue(b);
if (fa.DoubleEqual(0d, tolerance)) return a;
if (fb.DoubleEqual(0d, tolerance)) return b;
int signA = fa > 0 ? 1 : -1;
int signB = fb > 0 ? 1 : -1;
if (signA * signB > 0)
{
//중근인 경우에는 미분함수의 근을 구한다. (단, 둘 다 루트격리된 영역에서)
var newRange = StrumChain.IsolateAlsoDifferentiateRoot(this, doubleRange.min, doubleRange.max);
return GetDifferentiated().GetARootInRange_Illinois(newRange, tolerance);
}
int aCnt = 0;
int bCnt = 0;
int stagnationCnt = 2; //두번 연속 히트하면 정체점.
while (false == fa.DoubleEqual(fb, tolerance))
{
double slope = (fb - fa) / (b - a);
double constant = fa - slope * a;
double m = - constant / slope;
double fm = GetValue(m);
if (fm.DoubleEqual(0d, tolerance)) return m;
//signA와 부호가 같음.
if (fm * signA > 0)
{
a = m;
fa = fm;
aCnt = 0;
if (++bCnt >= stagnationCnt)
fb /= 2.0d;
}
else
{
b = m;
fb = fm;
bCnt = 0;
if (++aCnt >= stagnationCnt)
fa /= 2.0d;
}
}
return (fb + fa) / 2.0d;
}
3. 뉴턴법 (newton method)
뉴턴법의 경우 bracket이 필요 없다. 대신 함수가 미분 가능해야 한다. 역시 그림으로 보자.
뉴턴법의 경우 α=2 차수로 수렴한다. 즉 이분법보다 제곱의 속도로 빨리 수렴한다는 뜻.
그러나 제약이 굉장히 많은데 우선 함수가 미분가능해야한다. 또한 미분 값이 0에 가까운 경우 오히려 근에서 멀어지게 되며(발산) 사이클이 생겨서 영영 수렴하지 못하기도 한다. 따라서 매우 unstable 하기 때문에 단독으로 쓰기는 힘들다.
나도 이걸 써보려고 했는데 아마 2차 도함수의 근이 없는 구간이라면 적절한 뉴턴법의 초기점을 구할 수 있을 듯하다. 물론 증명은 안 했다.
4. 할선법 (secant method)
할선법은 뉴턴법과 비슷하지만 미분값 대신 두 점을 지나는 직선을 사용한다. 이 경우 뉴턴법과 달리 함수가 미분 가능하지 않거나 값들의 테이블로 주어진 경우에도 적용이 가능하다.
가위치법과 헷갈리지 말아야 한다. 할선법 역시 뉴턴 방법처럼 외부로 발산할 수 있으며 수렴이 보장되지 않는다. (unstable)
할선법의 경우 α=1.618 차수로 수렴한다.(superlinear) 때문에 뉴턴법보다 느리다고 생각할 수 있는데, 뉴턴법은 미분함수를 구해서 미분값을 구해야 하는 과정이 있으므로 만약 뉴턴법 2번할 때 할선법 3번을 할 수 있다면 할선법이 더 빠르게 된다.
위 방법들을 조합한 방법
위에 적힌 방법들을 그대로 사용하는 경우는 드물다. 특히 할선법이나 뉴턴법 같은 경우는 unstable하기 때문에 stable한 이분법이나 가위치법과 섞어 쓰기도 한다.
1. Brent's method (dekker's method)
Brent's method는 dekker's method를 기반으로 만든 방법이다.
dekker's method는 할선법+이분법 이라고 보면 된다. 할선법으로 진행하다가 막히면 이분법을 시도하는 셈. 이 경우 평균적으로 할선법과 같은 속도를 지니지만 최악의 경우에도 이분법의 속도를 가지며 안정성을 보장할 수 있다.
brent's method는 dekker's method에서 아이디어를 얻어 만들었다고 한다. (그래서 brent-dekker method라고도 함.) 이 경우 IQI(iverse quadratic interpolation)와 가위치법을 같이 섞어 쓴다.
참고로 MATLAB의 fzero가 이 방법을 이용한다고 할 정도로 뭔가 상업적으로 보장된? 성능을 가진 듯.
2. ITP 메서드
할선법과 가위치법을 섞어 이것저것 하는 메서드... 이것도 꽤 효율적이라고 한다.
참고 :
https://en.wikipedia.org/wiki/Real-root_isolation
https://en.wikipedia.org/wiki/Root-finding_algorithms
https://en.wikipedia.org/wiki/Regula_falsi
https://en.wikipedia.org/wiki/Secant_method
https://en.wikipedia.org/wiki/ITP_method
'수학 > 이론' 카테고리의 다른 글
| 한 점과 가장 가까운 곡선 위의 점 구하기 (0) | 2021.09.04 |
|---|---|
| 스트룸 정리 (Strum's theorem) (1) | 2021.08.13 |