이차방정식의 근을 (근의 공식을 써서 수치적으로) 안전하게 구하기

* 온라인 이차방정식 풀이

* 온라인 이차방정식 풀이 (2)

* 온라인 방정식 풀이 (3)



이차방정식   x^2 + 2x + 1 = 0   처럼 간단한 것은 보통의 방법으로 풀어도 된다.

하지만   0.00000001 x^2 + 2x + 0.00000001 = 0   처럼 생긴 것 (즉,   ax^2 + bx + c = 0   의 판별식   b^2 - 4ac   가 근사적으로   b^2   이 될 정도로   4ac   가 매우   0   에 가까운 것)은 보통의 근의 공식을 쓰면 수치적으로 오차가 있는 근을 구하게 된다.

이런 경우 (수치적 오차를 줄일 수 있는) 보다 안전한 해법을  생각해보자. 단, 수치적 계산에서 기계입실론 수(즉 유효숫자 개수의 한계로 인해   1   과의 덧셈에서 무시되는 최대의 양수값) 처럼 무시되는 오류는 어쩔 수 없다.

유효숫자 개수를 16개로 한정한 근사값 계산에서   1 + 10^{-20}   의 게산 결과는   1   이다. 여기서는   10^{-20}   이 기계입실론 수 보다 작기 때문에 수치적 덧셈에서 무시되기 때문이다.
  5 \cdot 10^{20} + 3   의 덧셈에서도 이는   5 \cdot 10^{20}   으로 묶으면   5 \cdot 10^{20} (1 + 6 \cdot 10^{-21})   이므로 역시   6 \cdot 10^{-21}   이 기계입실론 수 보다 작기 때문에 결국   3   은   5 \cdot 10^{20}   과 합하는 덧셈에서 무시된다.


다시 원래의 이차방정식   0.00000001 x^2 + 2x + 0.00000001 = 0   (즉   10^{-8} x^2 + 2x + 10^{-8} = 0  ) 으로 돌아가보자.

판별식   D = 2^2 - 4 \cdot 10^{-8} \cdot 10^{-8} = 3.9999999999999996   는 거의   4   에 가깝다.
근의 공식을 쓰면   x = \frac{-2 \pm \sqrt{D}}{2 \cdot 10^{-8}}   이 되어 분자의 복부호에서   +부호인 것은 그 계산 결과가 거의   0   에 가깝다. 이런 경우는 보통의 근의 공식을 직접 사용하면 수치적 오차가 있는 근이 구해진다.


우선 유효숫자 개수가 16개가 되도록 반올림하여서

     \sqrt{D} = \sqrt{3.9999999999999996} \approx \sqrt{4.000000000000000} \approx 2.0

이므로, 근사직으로 두 근은

     x_1 = \frac{-2 - \sqrt{D}}{2 \cdot 10^{-8}} \approx \frac{-4}{2 \cdot 10^{-8}} \approx -2 \cdot 10^8

과

     x_2 = \frac{-2 + \sqrt{D}}{2 \cdot 10^{-8}} \approx \frac{0}{2 \cdot 10^{-8}} \approx 0

이디.

그런데 수치적으로 계산할 때는 이 두 번 째 근사적 근   0   이 잘 얻어지지 않는다. 그 이유는 컴퓨터 계산에서는 2진법을 사용하므로 정확하게 유효숫자 10진법 16자리 근사적 계신이 아니라 유효숫자 2진법 N자리 근사적 계산이기 때문에

     \sqrt{3.9999999999999996} \approx \sqrt{4.000000000000000}

라는 계산과는 다르게 진행된다.


참고로, 주어진 이차방정식은 계수의 대칭성 때문에 만일   x_1   이 근이면 그 역수   \frac{1}{x_1}   도 근이다.

그러므로   \frac{1}{(-2) \cdot 10^8} = -5 \cdot 10^{-9}   가 근이 되어야 함은 마땅하다. 그런데 x_2 \approx 0 이 근이라고 구해졌으니 여기에는 약간의 문제가 있음을 알 수 있다. 그 이유는 x_2 를 구하는 과정에서 -b + \sqrt{D} (단, D는 판별식) 이   0   에 너무 가깝기 때문이다. 그래서 계산 과정에서 이런 경우가 가급적 나타나지 않게 하여 x_2 를 구하는 방법을 찾아본다.



a, \ b, \ c   가 실수이고   ac \ne 0, \ b^2 - 4ac \approx 0   이라고 가정하자.

이차방정식   a x^2 + b x + c = 0   의 근을 구하기 위해서

       q = \begin{cases} - b - \sqrt{b^2 - 4ac} \qquad \textrm{($b \ge 0$일 때)} \\ - b + \sqrt{b^2 - 4ac} \qquad \textrm{($b < 0$일 때)} \\ \end{cases}

처럼  q  를 정하고,   x_1 = \frac{q}{2a}   와   x_2= \frac{2c}{q}  로 두 근   x_1 , \ x_2   를 구하면 수치적으로 오차가 작은 근을 구할 수 있다.

(참고로, 복소수 계수의 이차방정식의 (오차를 줄이기 위한) 수치적 해법은 이것과 조금 다르다.   q = - \left( b + \textrm{(sign)} \sqrt{b^2 - 4ac} \right)   에서 부호   \textrm{(sign)}   은   \overline{b} \cdot \textrm{(sign)} \sqrt{b^2 - 4ac} > 0   되게 정하여   q   를 정한 다음   x_1 = \frac{q}{2a},   x_2 = \frac{2c}{q}   를 써서 두 근   x_1   과   x_2   를 구한다. 여기서   \overline{b}   는   b   의 켤레복소수를 의미한다.)



* 위의 과정을 확인햐기 위해 Python 2.7 인터프리터로 계산 과정을 추적해 보았다.
  (아래에서 두 번 째 근 x_2 가 0 이 아니라 -1.11022302463e-08 로 출력되었다.)

>>> from math import *
>>> a = 1e-8
>>> b = 2.0
>>> c = 1e-8
>>> print a, b, c
1e-08 2.0 1e-08
>>> d = b*b - 4.0*a*c
>>> print d
4.0
>>> print sqrt(d)
2.0
>>> q1 = -b - sqrt(d)
>>> q2 = -b + sqrt(d)
>>> print q1, q2
-4.0 -2.22044604925e-16
>>> x1 = q1/(2.0*a)
>>> x2 = q2/(2.0*a)
>>> print x1, x2
-200000000.0 -1.11022302463e-08

* 왜 위와 같은 결과가 나왔는지 알아보기 위해 이번에는 다소 다르게 테스트해보았다.

>>> from math import *
>>> t = 3.9999999999999996
>>> t
3.9999999999999996
>>> u = sqrt(t)
>>> u
1.9999999999999998
>>> -2 + u
-2.220446049250313e-16
>>> (-2 + u)/2e-8
-1.1102230246251565e-08

*  다음은 C 언어로 작성해본 것인데 위의 Python 인터프리터로 알아본 결과와 비슷하다.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct ROOTS {
 double re1;
 double im1;
 double re2;
 double im2;
};

void printRoot(const char *msg, struct ROOTS *data, int ndx) {
 printf(msg);
 if (ndx == 0) {
  if (data->im1 == 0.0)
         printf("%e\n", data->re1);
     else if (data->re1 == 0.0)
         printf("%e j\n", data->im1);
     else
         printf("%e + %e j\n", data->re1, data->im1);
 }
 else {
  if (data->im2 == 0.0)
         printf("%e\n", data->re2);
     else if (data->re2 == 0.0)
         printf("%e j\n", data->im2);
     else
            printf("%e + %e j\n", data->re2, data->im2);
    }
}

struct ROOTS *solveQuadraticPlain(double a, double b, double c) {
    struct ROOTS *roots = (struct ROOTS *)calloc(sizeof(struct ROOTS), 1);
    double d = b*b - 4.0*a*c;
    double r, x1, x2;
    if (d == 0.0) {
        x1 = -b/(2.0*a);
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x1;
        roots->im2 = 0.0;
 }
 else if (d > 0.0) {
     r = sqrt(d);
        x1 = (-b - r)/(2.0*a);
        x2 = (-b + r)/(2.0*a);
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x2;
        roots->im2 = 0.0;
 }
 else {
     r = sqrt(-d);
        x1 = -b/(2.0*a);
        x2 = r/(2.0*a);
        roots->re1 = x1;
        roots->im1 = -x2;
        roots->re2 = x1;
        roots->im2 = x2;
 }
    return roots;
}

int main() {
    struct ROOTS *sol = (struct ROOTS *)solveQuadraticPlain(1.0, 2.0, 1.0);
    printf("Use the plain quadratic root fomula:\n");
    printf("  The solutions of the equation x^2 + 2x + 1 = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    sol = (struct ROOTS *)solveQuadraticPlain(1.0e-8, 2.0, 1.0e-8);
    printf("\n");
    printf("  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    return 0;
}

/*
Execute Result:
Use the plain quadratic root fomula:
  The solutions of the equation x^2 + 2x + 1 = 0 are
    x1 = -1.000000e+000
    x2 = -1.000000e+000

  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are
    x1 = -2.000000e+008
    x2 = -1.110223e-008
*/

* Java 언어로 작성해 보아도 마찬가지이다.

import java.util.*;

public class SolveQuadraticEqPlain_01 {

    public static void printRoot(String msg, double[] data, int ndx) {
     System.out.printf(msg);
     if (ndx == 0) {
         if (data[1] == 0.0)
             System.out.printf("%e\n", data[0]);
         else if (data[0] == 0.0)
             System.out.printf("%e j\n", data[1]);
         else
             System.out.printf("%e + %e j\n", data[0], data[1]);
     }
     else {
         if (data[3] == 0.0)
             System.out.printf("%e\n", data[2]);
         else if (data[2] == 0.0)
             System.out.printf("%e j\n", data[3]);
         else
             System.out.printf("%e + %e j\n", data[2], data[3]);
     }
    }

    public static double[] solveQuadraticPlain(double a, double b, double c) {
  double[] roots = new double[4];
        double d = b*b - 4.0*a*c;
        double r, x1, x2;
        if (d == 0.0) {
            x1 = -b/(2.0*a);
            roots[0] = x1;
            roots[1] = 0.0;
            roots[2] = x1;
            roots[3] = 0.0;
     }
     else if (d > 0.0) {
         r = Math.sqrt(d);
            x1 = (-b - r)/(2.0*a);
            x2 = (-b + r)/(2.0*a);
        roots[0] = x1;
        roots[1] = 0.0;
        roots[2] = x2;
        roots[3] = 0.0;
     }
     else {
         r = Math.sqrt(-d);
            x1 = -b/(2.0*a);
            x2 = r/(2.0*a);
            roots[0] = x1;
            roots[1] = -x2;
            roots[2] = x1;
            roots[3] = x2;
     }

        return roots;
    }

    public static void main(String[] args) {
        double[]  sol = solveQuadraticPlain(1.0, 2.0, 1.0);
        System.out.printf("Use the plain quadratic root fomula:\n");
        System.out.printf("  The solutions of the equation x^2 + 2x + 1 = 0 are\n");
        printRoot("    x1 = ", sol, 0);
        printRoot("    x2 = ", sol, 1);

        sol = solveQuadraticPlain(1.0e-8, 2.0, 1.0e-8);
        System.out.printf("\n");
        System.out.printf("  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are\n");
        printRoot("    x1 = ", sol, 0);
        printRoot("    x2 = ", sol, 1);
    }
}

/*
Execute Result:
Use the plain quadratic root fomula:
  The solutions of the equation x^2 + 2x + 1 = 0 are
    x1 = -1.000000e+00
    x2 = -1.000000e+00

  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are
    x1 = -2.000000e+08
    x2 = -1.110223e-08
*/

* 참고로 아래는 Octave 로 이차방정식   x^2 + 2x + 1 = 0   과   10^{-8} x^2 + 2x + 10^{-8} = 0   의 근을 구해본 것이다.
   ( -2.0000e+008 의 역수  -5.0000e-009 도 근으로 출력되었다. )

octave.exe:1> c = [1 2 1];
octave.exe:2> roots(c)
ans =

  -1
  -1

octave.exe:3> c = [1e-8 2 1e-8];
octave.exe:4> roots(c)
ans =

  -2.0000e+008
  -5.0000e-009

* 또 아래는 Python 의 numpy 를 이용하여 이차방정식   10^{-8} x^2 + 2x + 10^{-8} = 0   의 근을 구해본 것이다.
   ( -200000000.0 의 역수 -5.0000000000000001e-09 도 근으로 출력되었다. )

from numpy import poly1d

P = poly1d([1e-8, 2.0, 1e-8], variable = 'x')

print("Polynomial p(x) =")
print(P)
print("\nThe complex roots of p(x) = 0 are as a numpy array:")
print(P.r)
print("\nThe complex roots of p(x) = 0 are as a python list:")
print(list(P.r))

"""
Output:
Polynomial p(x) =
       2
1e-08 x + 2 x + 1e-08

The complex roots of p(x) = 0 are as a numpy array:
[ -2.00000000e+08  -5.00000000e-09]

The complex roots of p(x) = 0 are as a python list:
[-200000000.0, -5.0000000000000001e-09]
"""

아래의 두 그림은 각각 Mathematica 와 Maxima 를 써서 이차방정식을 풀어본 결과이다.

* Mathematica 로 이차방정식 풀이
  ( -5.×10^{-8} 이 근으로 출력되었다. )

* Maxima 로 이차방정식 풀이
   ( Maxima 에서는 0 이 근으로 출력되었다. )

* 또 아래 그림은 온라인 이차방정식 풀이 계산기  math.com 을 사용한 결과이다.

* 이차방정식의 근을 (수치적으로 안전하게) 구하기 위해 다시 작성된 C 언어 소스

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct ROOTS {
 double re1;
 double im1;
 double re2;
 double im2;
};

void printRoot(const char *msg, struct ROOTS *data, int ndx) {
 printf(msg);
 if (ndx == 0) {
  if (data->im1 == 0.0)
         printf("%e\n", data->re1);
     else if (data->re1 == 0.0)
         printf("%e j\n", data->im1);
     else
         printf("%e + %e j\n", data->re1, data->im1);
 }
 else {
  if (data->im2 == 0.0)
         printf("%e\n", data->re2);
     else if (data->re2 == 0.0)
         printf("%e j\n", data->im2);
     else
            printf("%e + %e j\n", data->re2, data->im2);
    }
}

struct ROOTS *solveQuadraticPlain(double a, double b, double c) {
    struct ROOTS *roots = (struct ROOTS *)calloc(sizeof(struct ROOTS), 1);
    double d = b*b - 4.0*a*c;
    double r, x1, x2;
    if (d == 0.0) {
        x1 = -b/(2.0*a);
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x1;
        roots->im2 = 0.0;
 }
 else if (d > 0.0) {
     r = sqrt(d);
        x1 = (-b - r)/(2.0*a);
        x2 = (-b + r)/(2.0*a);
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x2;
        roots->im2 = 0.0;
 }
 else {
     r = sqrt(-d);
        x1 = -b/(2.0*a);
        x2 = r/(2.0*a);
        roots->re1 = x1;
        roots->im1 = -x2;
        roots->re2 = x1;
        roots->im2 = x2;
 }
    return roots;
}

struct ROOTS *solveQuadraticSafe(double a, double b, double c) {
    struct ROOTS *roots = (struct ROOTS *)calloc(sizeof(struct ROOTS), 1);
    double d = b*b - 4.0*a*c;
    double r = sqrt(b*b - 4.0*a*c);
    double q, x1, x2;
    if (d == 0.0) {
        x1 = -b/(2.0*a);
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x1;
        roots->im2 = 0.0;
 }
 else if (d > 0.0) {
     r = sqrt(d);
     if (b > 0.0) {
           q = -b - r;
  }
  else {
           q = -b + r;
  }
        x1 = q/(2.0*a);
        x2 = (2.0*c)/q;
        roots->re1 = x1;
        roots->im1 = 0.0;
        roots->re2 = x2;
        roots->im2 = 0.0;
 }
 else {
     r = sqrt(-d);
        x1 = -b/(2.0*a);
        x2 = r/(2.0*a);
        roots->re1 = x1;
        roots->im1 = -x2;
        roots->re2 = x1;
        roots->im2 = x2;
 }

    return roots;
}

int main() {
    struct ROOTS *sol= (struct ROOTS *)solveQuadraticPlain(1.0, 2.0, 1.0);
    printf("Use the plain quadratic root fomula:\n");
    printf("  The solutions of the equation x^2 + 2x + 1 = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    sol= (struct ROOTS *)solveQuadraticPlain(1.0e-8, 2.0, 1.0e-8);
    printf("\n");
    printf("  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    sol= (struct ROOTS *)solveQuadraticSafe(1.0, 2.0, 1.0);
    printf("\n");
    printf("Use the safe quadratic root fomula:\n");
    printf("   The solutions of the equation x^2 + 2x + 1 = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    sol= (struct ROOTS *)solveQuadraticSafe(1.0e-8, 2.0, 1.0e-8);
    printf("  \n");
    printf("The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are\n");
    printRoot("    x1 = ", sol, 0);
    printRoot("    x2 = ", sol, 1);
    free(sol);

    return 0;
}

/*
Execute Result:
Use the plain quadratic root fomula:
  The solutions of the equation x^2 + 2x + 1 = 0 are
    x1 = -1.000000e+000
    x2 = -1.000000e+000

  The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are
    x1 = -2.000000e+008
    x2 = -1.110223e-008

Use the safe quadratic root fomula:
   The solutions of the equation x^2 + 2x + 1 = 0 are
    x1 = -1.000000e+000
    x2 = -1.000000e+000

The solutions of the equation 10^(-8) x^2 + 2x + 10^(-8) = 0 are
    x1 = -2.000000e+008
    x2 = -5.000000e-009
*/