* 온라인 이차방정식 풀이
* 온라인 이차방정식 풀이 (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
*/