수치해석 알고리즘 : 이분법

카테고리 없음

2019. 4. 22. 23:50

앞서 설계한 다항식 클래스에 이분법 알고리즘인 bisect 메서드를 추가할 것이다.

 

이분법은 하나의 다항식이 주어졌을 때, [start,end] 구간에서 근을 하나 가진다고 가정했을 때 범위를 반으로 쪼개어가며 수치적 근사해를 구하는 기법이다.

 

1. 중간값 mid = (start+end)/2를 구한다. 

2. mid * end의 부호를 판별한다.

3. 곱한 값이 0보다 크면 근은 start와 mid 사이에 있고, 0보다 작으면 mid와 end 사이에 있다.(만약 0이라면 mid값이 근이다.)

4. 전자의 경우 mid를 end로, 후자의 경우 start로 하여 계속 반복한다.

 

이를 적용한 Polynomial 클래스를 보자.

* poly.h

...

class Polynomial
{
protected:
    vector<Term> terms;
public:
    Polynomial();
    double yVal(double x) const;
    double bisect(double start, double end, double tol) const;
    void show() const;
};

...

* poly.cpp

...

double Polynomial::yVal(double x) const
{
    double sum = 0;
    for(auto it=terms.begin(); it!=terms.end(); ++it) {
        double val = 1.0;
        for(int i=0; i<it->getDeg(); ++i) {
            val *= x;
        }
        val *= it->getCoef();
        sum += val;
    }   

    return sum;
}

double Polynomial::bisect(double start, double end, double tol) const
{
    double mid;

    if(yVal(start) * yVal(end) >= 0) {
        cout << "f(a)f(b) >= 0 !" << endl;
        exit(-1);
    }   

    mid = (start + end) / 2.0;
    while(abs(end - start) > tol) {
        double mulVal = yVal(mid) * yVal(end);
        if(mulVal < 0)
            start = mid;
        else if(mulVal > 0)
            end = mid;
        else
            break;
        mid = (start + end) / 2.0;
    }   

    return mid;
}

...

먼저 구간의 양 끝점에서의 함수값을 구할 필요가 있으므로 주어진 x에서의 함수값을 구해주는 yVal(double x) 메서드를 구현하였다.

이는 각 항들을 순회하면서 차수만큼 인자로 전달받은 x를 곱하고, 차수만큼 곱했으면 계수를 곱하고, sum에 누적하는 방식으로 함수값을 구한다.

 

bisect(double start, double end, double tol)은 주어진 두 구간[start, end]에서 허용오차를 tol로 했을 때의 이분법으로 수치해를 구해주는 메서드이다.

먼저 f(start)f(end)를 구해 0보다 작은지 검사한 후에, 두 구간의 차이가 허용오차보다 작아질 때까지 위에서 설명한 이분법의 알고리즘 순서를 반복한다.

 

이를 적용한 예제는 다음과 같다.

 

* bisect.cpp

#include <iostream> 
#include <poly.h>
using namespace std;

int main(void)
{
    Polynomial poly;
    double start, end, tol;
    cout << "type two boundary values." << endl;
    cin >> start >> end;

    cout << "type tolerance." << endl;
    cin >> tol;

    cout << "sol : " << poly.bisect(start, end, tol) << endl;
    return 0;
}

다항식 x^3 - x^2 -1 이 [1,2] 구간에서 하나의 근을 가질 때 허용오차를 0.0001로 하여 프로그램을 실행시켜 보았다.

 

* 실행결과

$ ./a.out
how many terms are there ?
3
type coefficient and degree of 1th term.
1 3
type coefficient and degree of 2th term.
-1 2
type coefficient and degree of 3th term.
-1 0
type two boundary values.
1 2
type tolerance.
0.0001
sol : 1.46555

비교적 정확한 수치적 근사해를 구해내었다.