一、问题描述
已知一元二次方程
a
x
2
+
b
x
+
c
=
0
(
a
≠
0
)
ax^2 + bx + c = 0(a\neq0)
ax2+bx+c=0(a=0),求方程的实数根。在数学上,我们很容易可以得到方程的实数根为:
x
1
=
−
b
+
b
2
−
4
a
c
2
a
(1)
x_1=\frac{-b+\sqrt{b^2-4ac}}{2a} \tag 1
x1=2a−b+b2−4ac(1)
x
2
=
−
b
−
b
2
−
4
a
c
2
a
(2)
x_2=\frac{-b-\sqrt{b^2-4ac}}{2a} \tag 2
x2=2a−b−b2−4ac(2)
其中,
b
2
−
4
a
c
≥
0
b^2-4ac\geq0
b2−4ac≥0
然而,在计算机上,由于计算机字长的限制,有些浮点数无法绝对精确地存储。并且,一些计算过程可能会导致精度损失。例如接下来这两个例子。
(
A
)
(A)
(A)设有两个数
p
=
3.1415926536
p=3.1415926536
p=3.1415926536和
q
=
3.1415957341
q=3.1415957341
q=3.1415957341,两者几乎相等,而且都有11位有效数字精度。它们的差为
p
−
q
=
−
0.0000030805
p-q=-0.0000030805
p−q=−0.0000030805,只有5位有效数字精度!
(
B
)
(B)
(B)设
f
(
x
)
=
x
(
x
+
1
−
x
)
,
g
(
x
)
=
x
x
+
1
+
x
f(x)=x(\sqrt{x+1}-\sqrt{x}),g(x)=\frac{x}{\sqrt{x+1}+\sqrt{x}}
f(x)=x(x+1−x),g(x)=x+1+xx。显然,当
x
≥
0
x\geq0
x≥0时,
f
(
x
)
f(x)
f(x)可以通过分母有理化得到
g
(
x
)
g(x)
g(x),在数学上两者是等价的。现在,我们采用6位有效数字精度和舍入法比较
f
(
500
)
f(500)
f(500)与
g
(
500
)
g(500)
g(500)的计算结果。
f
(
500
)
=
500
×
(
501
−
500
)
=
500
×
(
22.3830
−
22.3607
)
=
500
×
0.0223
=
11.1500
f(500)=500\times(\sqrt{501}-\sqrt{500})=500\times(22.3830-22.3607)=500\times0.0223=11.1500
f(500)=500×(501−500)=500×(22.3830−22.3607)=500×0.0223=11.1500。
g
(
500
)
=
500
501
+
500
=
500
22.3830
+
22.3607
=
500
44.7437
=
11.1748
g(500)=\frac{500}{\sqrt{501}+\sqrt{500}}=\frac{500}{22.3830+22.3607}=\frac{500}{44.7437}=11.1748
g(500)=501+500500=22.3830+22.3607500=44.7437500=11.1748。
g
(
500
)
g(500)
g(500)比
f
(
500
)
f(500)
f(500)计算精度更高!(真实值为11.174755300747198…,舍入到6位有效数字,与
g
(
500
)
g(500)
g(500)的计算值相同)
精度损失的现象在我们编程过程中很有可能会出现,有时可能会使计算结果与正确结果大相径庭,产生灾难性的后果,我们必须引起足够重视。精度损失的现象有时候是没办法避免的,有时候则可以通过一些数值计算算法或技巧来避免。显然,式
(
1
)
(1)
(1)和式
(
2
)
(2)
(2)的分子在一定条件下,可能是两个相近的数相减,影响计算精度。那么,如何编写一个计算精度较高的程序来求一元二次方程的实数根呢?
二、推导步骤
编写一个计算精度较高的程序来求一元二次方程
a
x
2
+
b
x
+
c
=
0
(
a
≠
0
)
ax^2 + bx + c = 0(a\neq0)
ax2+bx+c=0(a=0)的实数根,至少需要考虑两种情形:
(
A
)
∣
b
∣
≈
b
2
−
4
a
c
(
B
)
c
=
0
(A)\left|b\right| \approx\sqrt{b^2-4ac}\ \ (B)c=0
(A)∣b∣≈b2−4ac (B)c=0
(
A
)
(A)
(A)若
∣
b
∣
≈
b
2
−
4
a
c
\left|b\right| \approx\sqrt{b^2-4ac}
∣b∣≈b2−4ac,当
b
>
0
b>0
b>0时,式
(
1
)
(1)
(1)中,
−
b
+
b
2
−
4
a
c
-b+\sqrt{b^2-4ac}
−b+b2−4ac会导致计算精度损失,对式
(
1
)
(1)
(1)分子有理化得:
x
1
=
−
2
c
b
+
b
2
−
4
a
c
(3)
x_1=\frac{-2c}{b+\sqrt{b^2-4ac}} \tag 3
x1=b+b2−4ac−2c(3)
当
b
≤
0
b\leq0
b≤0时,式
(
2
)
(2)
(2)中,
−
b
−
b
2
−
4
a
c
-b-\sqrt{b^2-4ac}
−b−b2−4ac会导致计算精度损失,对式
(
2
)
(2)
(2)分子有理化得:
x
2
=
−
2
c
b
−
b
2
−
4
a
c
(4)
x_2=\frac{-2c}{b-\sqrt{b^2-4ac}} \tag 4
x2=b−b2−4ac−2c(4)
(
B
)
(B)
(B)若
c
=
0
c=0
c=0,在数学上,
b
2
−
4
a
c
=
∣
b
∣
\sqrt{b^2-4ac}=\left|b\right|
b2−4ac=∣b∣,但计算机在数值计算过程中,对
b
b
b进行了开方和平方根运算,会引入计算误差(或截断误差)和舍入误差,导致此时
b
2
−
4
a
c
≠
∣
b
∣
\sqrt{b^2-4ac}\neq\left|b\right|
b2−4ac=∣b∣,使方程根的计算误差变大。因而,当
c
=
0
c=0
c=0时,需特殊处理,此时方程的两实数根为:
{ x 1 = 0 x 2 = − b a (5) \left\{ \begin{array}{c} x_1=0 \\ \tag 5 x_2=\frac{-b}{a}\end{array}\right. {x1=0x2=a−b(5)
三、 C C C代码
#include <math.h>
#include <float.h>
#include <stdio.h>
#define UINT32 unsigned int
#define ERR_NO_ERROR 0x00000000
#define ERR_NAN 0x00000001
#define ERR_INF 0x00000002
#define MAX(a, b) ((a) > (b)) ? (a) : (b)
#define MIN(a, b) ((a) < (b)) ? (a) : (b)
/*************************************************
Function: is_number
Description: 判断浮点数是否为nan
Input: 浮点数x
Output: 无
Return: 若浮点数x为nan返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_number(float x)
{
return (x == x);
}
/*************************************************
Function: is_finite_number
Description: 判断浮点数是否为inf
Input: 浮点数x
Output: 无
Return: 若浮点数x为inf返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_finite_number(float x)
{
return (x >= -FLT_MAX && x <= FLT_MAX);
}
/*************************************************
Function: solve_quadratic_equation
Description: 求一元二次方程(a*x^2 + b*x + c = 0)的所有实数根
Input: 方程的系数 p = {c, b, a}
Output: 方程的所有实数根x, 实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_quadratic_equation(float p[], float x[], int *rootCount)
{
int i;
float a, b, c, delta, sqrtDelta;
const float ZERO = FLT_MIN; // min normalized positive value(1.175494351e-38F)
const float EPS = FLT_MIN;
UINT32 errNo = ERR_NO_ERROR;
*rootCount = 0;
for (i = 0; i < 3; i++)
{
if (!is_number(p[i]))
{
errNo = ERR_NAN;
return errNo;
}
if (!is_finite_number(p[i]))
{
errNo = ERR_INF;
return errNo;
}
}
a = p[2];
b = p[1];
c = p[0];
if (fabs(a - 0.0) < EPS)
{
if (fabs(b - 0.0) > EPS)
{
x[0] = -c / b;
*rootCount = 1;
}
}
else
{
b /= a;
c /= a;
a = 1.0;
delta = b * b - 4.0 * a * c;
if (delta > ZERO)
{
if (fabs(c - 0.0) < EPS) //若c = 0,由于计算误差,sqrt(b*b - 4*a*c)不等于|b|
{
x[0] = 0.0;
x[1] = -b / a;
}
else
{
sqrtDelta = sqrt(delta);
if (b > 0.0)
{
x[0] = (-2.0 * c) / (b + sqrtDelta); //避免两个很接近的数相减,导致精度丢失
x[1] = (-b - sqrtDelta) / (2.0 * a);
}
else
{
x[0] = (-b + sqrtDelta) / (2.0 * a);
x[1] = (-2.0 * c) / (b - sqrtDelta); //避免两个很接近的数相减,导致精度丢失
}
}
*rootCount = 2;
}
else if (fabs(delta - 0.0) < EPS)
{
x[0] = x[1] = -b / (2.0 * a);
*rootCount = 2;
}
else
{
*rootCount = 0;
}
}
return errNo;
}
void main(void)
{
float x[2], p[3];
int rootCount;
float x1, x2, a, b, c;
UINT32 errNo = ERR_NO_ERROR;
//一元二次方程测试
//(1)(x - 1000)*(x - 0.001) = x^2 - 1000.001*x + 1 = 0
a = 1;
b = -1000.001;
c = 1;
//(2) 3*x ^ 2 - 1000000*x + 0 = 0
//a = 3;
//b = -1000000;
//c = 0;
//(3) 1.0e-10*x^2 - 2.0e-10*x + 1.0e-10 = 0
//a = 1.0e-10;
//b = -2.0e-10;
//c = 1.0e-10;
p[0] = c;
p[1] = b;
p[2] = a;
errNo = solve_quadratic_equation(p, x, &rootCount);
x1 = (-b - sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
x2 = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
printf("原始方法求得方程的根:x1=%f, x2=%f\n本文方法求得方程的根:x1=%f, x2=%f\n",
MIN(x1, x2), MAX(x1, x2), MIN(x[0], x[1]), MAX(x[0], x[1]));
}
(1)当
a
=
1
,
b
=
−
1000.001
,
c
=
1
a=1,b=-1000.001,c=1
a=1,b=−1000.001,c=1时,方程根的准确值为:0.001和1000,程序运行结果:
(2)当
a
=
3
,
b
=
−
1000000
,
c
=
0
a=3,b=-1000000,c=0
a=3,b=−1000000,c=0时,方程根的准确值为:0和333333.33…,程序运行结果:
(3)当
a
=
1.0
e
−
10
,
b
=
−
2.0
e
−
10
,
c
=
1.0
e
−
10
a=1.0e^{-10},b=-2.0e^{-10},c=1.0e^{-10}
a=1.0e−10,b=−2.0e−10,c=1.0e−10时,方程根的准确值为:1和1,程序运行结果:
四、总结
一元二次方程实数根的c语言实现,在很多人看来非常简单,根本不值一提,然而就是这么简单的问题,要写出计算精度高、简洁优雅、可读性强、健壮的代码也不是一件简单的事。
五、参考文献/资料
Numerical Methods Using MATLAB Fourth Edition. John Mathews, Kurtis D.Fink
数值方法(MATLAB版)(第四版) 周璐,陈渝,钱方等译