2017-03-20 59 views
0

我试图在C中创建一个Gauss Eliminator。为此,我不时需要检查一个矩阵是否在数值上是单数的:如果某个数(双)非常非常小。什么是推荐最小epsilon为双?

我的问题是,如果我尝试这样做:

if(0 == matrix->items[from]){ 
     fprintf(stderr,"Matrix is still singular after attempting pivot. Exitig.\n"); 
} 

这不会产生如此。由于double的不准确性,它永远不会完全是0.但是,当试图运行程序时,像这样的情况用inf或NaN填充数字,取决于是乘以还是除以它及其组合。

为了过滤这些,我需要这样的:

#define EPSILON very_small 
// rest of the code 
if(matrix->items[from] < EPSILON){ 
    ...singular 
} 

,这是什么EPSILON推荐值?这是双倍的绝对准确度,还是更大的价值?

顺便说一句,这将是更好的,它定义为一个宏如上,或者用它喜欢:

const double EPSILON = ...; 

很抱歉,如果我不是很清楚,英语不是我的母语。

感谢您的回复。

+3

请不要写'if(0 == matrix-> items [from])',这真的很难看。如果你不小心使用赋值而不是比较,现代编译器会投诉,所以这不再是“*好*”的做法。 –

+0

只有一个正确的值:'DBL_EPSILON'。 – Olaf

+2

@Olaf为什么只有这个唯一正确的值? – immibis

回答

2

我需要检查矩阵是否数值奇异

通常,这是通过防止double溢出检测。

// Check if 1.0/determinant will overflow. 
if (fabs(determinant) <= 1.0/(0.99*DBL_MAX)) { 
    Handle_Singular_Case() 
} else { 
    one_over_det = 1.0/determinant; 
} 

使用DBL_EPSILON(例如:2E-16)通常是溶液。 double数学需要相对比较,以确保从1.0的幅度很远的好计算。

// Rarely the right thing to do. 
#define EPSILON DBL_EPSILON 
if(fabs(matrix->items[from]) < EPSILON){ 

然而,这是非常上下文敏感@Weather Vane


然而OP真实的问题肯定是在这里:“试图运行程序的时候,这样的情况下填写数字与INF或NaN,这取决于是否乘以或与它和它的组合分割。”可以使用各种技术来避免此问题,例如使用partial pivoting来消除此问题。

要解决问题,最好发布代码和示例数据。

+1

关于构造一般数值稳定性算法的好的观察,特别是关于这个目的的部分旋转。将数学转换为优秀的程序并不像将公式转换为所选语言的语法那么简单。 –

+1

@JohnBollinger真。矩阵运算的接受计算隐藏了各种陷阱。我在部分代码旋转代码中遇到了一个错误,因为在查找最大元素时,代码使用了'abs()'而不是正确的'fabs()' - 并且启用了不足的警告。直到所有|元素|都没有真正显示错误小于1.0 - 一个0.0的div发生在后面。 – chux

+0

'(1.0 * DBL_EPSILON)'是什么意思?乘以一个实际上可以做什么? – rici