計算並聯連接電阻

Question

下面2個代數等價公式並聯電阻：

par1(r1, r2) = (r1 * r2)/(r1 + r2), or 

par2(r1, r2) = 1/(1/r1 + 1/r2) 

繼兩個Python函數，其中的每一個計算parallel_resistors式：

def par1(r1, r2): 
    return div_interval(mul_interval(r1, r2), add_interval(r1, r2)) 

def par2(r1, r2): 
    one = interval(1, 1) 
    rep_r1 = div_interval(one, r1) 
    rep_r2 = div_interval(one, r2) 
    return div_interval(one, add_interval(rep_r1, rep_r2)) 

下面是一個區間運算抽象由前述功能par1和par2使用。

def interval(a, b): 
    """Construct an interval from a to b. """ 
    return (a, b) 

def lower_bound(x): 
    """Return the lower bound of interval x. """ 
    return x[0] 

def upper_bound(x): 
    """Return the upper bound of interval x. """ 
    return x[1] 

def div_interval(x, y): 
    """Return the interval that contains the quotient of any value in x divided 
    by any value in y. 

    Division is implemented as the multiplication of x by the reciprocal of y. 
    >>> str_interval(div_interval(interval(-1, 2), interval(4, 8))) 
    '-0.25 to 0.5' 
    """ 
    assert (lower_bound(y) > 0 or upper_bound(y) < 0), "what it means to divide by an interval that spans zero" 
    reciprocal_y = interval(1/upper_bound(y), 1/lower_bound(y)) 
    return mul_interval(x, reciprocal_y) 


def str_interval(x): 
    """Return a string representation of interval x. 
    >>> str_interval(interval(-1, 2)) 
    '-1 to 2' 
    """ 
    return '{0} to {1}'.format(lower_bound(x), upper_bound(x)) 

def add_interval(x, y): 
    """Return an interval that contains the sum of any value in interval x and 
    any value in interval y. 
    >>> str_interval(add_interval(interval(-1, 2), interval(4, 8))) 
    '3 to 10' 
    """ 
    lower = lower_bound(x) + lower_bound(y) 
    upper = upper_bound(x) + upper_bound(y) 
    return interval(lower, upper) 

def mul_interval(x, y): 
    """Return the interval that contains the product of any value in x and any 
    value in y. 

    >>> str_interval(mul_interval(interval(-1, 2), interval(4, 8))) 
    '-8 to 16' 
    """ 
    p1 = lower_bound(x) * lower_bound(y) 
    p2 = lower_bound(x) * upper_bound(y) 
    p3 = upper_bound(x) * lower_bound(y) 
    p4 = upper_bound(x) * upper_bound(y) 
    return interval(min(p1, p2, p3, p4), max(p1, p2, p3, p4)) 

測試結果：

>>> r1 = interval(1, 2) 
>>> r2 = interval(3, 4) 
>>> par1(r1, r2) 
(0.5, 2.0) 
>>> par2(r1, r2) 
(0.75, 1.3333333333333333) 

我們注意到從par1和par2了不同的結果，它通過計算不同，但代數等價表達式。 對於上面給出的輸入r1和r2，以下是計算。

par1 --> return mul_interval((3, 8), (1/6, 1/4)) = (1/2, 2) 
======= 
rep_r1 = div_interval((1, 1), (1, 2)) = (1/2, 1) 
rep_r2 = div_interval((1, 1), (3, 4)) = (1/4, 1/3) 
par2 --> return div_interval((1, 1), (3/4, 4/3)) = (3/4, 4/3) 

其原因不同間隔是由於IEEE浮點格式，其中每div_interval失去精度。

我的理解是否正確？

Answer 1

恐怕你推理爲什麼兩個函數給出不同的結果是不正確的。浮點舍入可能是一個真正的問題，但對於大多數計算（包括這裏的那些），錯誤非常非常小。只有在您測試完全相等或者結合大量浮點計算時，錯誤可能會變得很重要。

這裏真正的問題是，在par1中，您的r1和r2間隔都出現在除法操作的兩側。當你的除法函數計算其第二個參數的倒數時，它會顛倒區間邊界的順序，所以輸入的上界影響輸出的下界。

當相同的時間間隔是除數和分紅的一部分時，您會得到比您應該更寬的時間間隔。這是因爲代碼不知道對於間隔的副本之一的一個極端的結果必定也意味着另一個相同的極端。

考慮計算r/r，區間r = (a, b)。你的div_interval函數將計算結果爲(a/b, b/a)，當通過基本代數時，我們知道結果應該是(1, 1)（因爲任何數字除以它自己都是1）。然而，分割代碼必須假定它的每個參數都是獨立的。如果除數結果爲其間隔的最大界限（b），則可能會將股利除以其最小界限（a）。它不知道這兩件事情不可能同時發生。

第二個公式避免了這個問題，因爲每個輸入間隔只出現在一個地方。你實際上可以允許一點點重複。只要沒有減法，除以包含負邊界的區間的乘法，應該沒有問題。

Answer 2

爲了證明它不是舍入誤差，這裏有一些代碼使用兩種方法來計算並聯電阻。它並不打算簡潔或者特別有效，只是爲了說明兩種計算方法的結果完全相同。

def par1(r1,r2): 
    r1lower = float(r1[0]) 
    r1upper = float(r1[1]) 
    r2lower = float(r2[0]) 
    r2upper = float(r2[1]) 
    low = 1/((1/r1lower)+(1/r2lower)) 
    high = 1/((1/r1upper)+(1/r2upper)) 
    return(low, high) 

def par2(r1,r2): 
    r1lower = float(r1[0]) 
    r1upper = float(r1[1]) 
    r2lower = float(r2[0]) 
    r2upper = float(r2[1]) 
    low = (r1lower * r2lower)/(r1lower + r2lower) 
    high = (r1upper * r2upper)/(r1upper + r2upper) 
    return(low, high) 

if __name__ == '__main__': 
    print par1((1,2),(3,4)) 
    print par2((1,2),(3,4)) 

輸出：

(0.75, 1.3333333333333333) 
(0.75, 1.3333333333333333) 
[Finished in 0.1s] 

某處在你的代碼中有一個錯誤...