2017-03-20 73 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