计算并联电阻

Question

下面两个并联电阻的代数等效公式：

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

多个并联阻抗——万无一失

多个impedances in parallel 的总阻抗对应于这些阻抗的harmonic mean 除以并联阻抗的数量。

from scipy import stats

def par(array):
    """Calculate the impedance of an array of parallel impedances"""
    return stats.hmean(array)/len(array)

In [1]: par([2.7, 3.3, 10])
Out[1]: 1.29299

Answer 3

为了证明这不是舍入误差，这里有一些代码使用这两种方法来计算并联电阻。它并不是为了简洁或特别有效，只是为了表明两种计算方式产生完全相同的结果。

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]

您的代码中的某个地方有错误...