区间算数

区间算数

内容基于 SICP ，练习 2.7-2.16。

引言

区间算数是指一套定义在区间上的运算。用区间表示不确定量的取值范围，常用来衡量测量值的误差，如 35.6±0.2°C。带有不确定性的量参与运算，结果也是不确定量，区间算数就是用来形式化这样的运算。本文只涉及四则运算，讨论其数学原理、scheme 实现和潜在的问题。

数学原理

区间的存储结构是上下界组成的序对，如果结果的界是运算数的界的函数，编程会比较方便。

下面着重讨论结果的界是不是运算数的函数，是什么函数。记，
x ∈ [ l x , u x ] y ∈ [ l y , u y ] z = f ( x , y ) , f ∈ { + , − , × , ÷ } x\in[l_x,u_x]\\ y\in[l_y,u_y]\\ z=f(x,y),f\in\{+,-,\times,\div\} x∈[lx​,ux​]y∈[ly​,uy​]z=f(x,y),f∈{+,−,×,÷}
以区间$x$为横轴，$y$为纵轴，$z$为竖轴绘图。从$z$轴俯视：

四个顶点处的$z$值是我们容易求出的。下面结合图像分析四则运算。

加

令$z=x+y$，绘制其图像：

可以看出，最值点在立方体的/对角线上。左下最低，右上最高（将+应用到相应序对上）：
$$z\in [l_x+l_y,u_x+u_y]$$
想象一个在$z$平面上滑动的矩形，可以想象到，上述结论对所有可能的区间都成立，与区间的相对位置、区间的正负均无关。

减

令$z=x-y$，绘图：

可以看出，和加法类似，减法的最值点出现在\对角线上。左上最低，右下最高，运用减法在相应的序对上：
$$z\in [l_x-u_y,u_x-l_y]$$
与相对位置、正负无关。

乘

乘法比较特殊，如图：

想象滑动的矩形，至少有4中可能的最值情况（矩形分别在四个象限内）。但由于平面内没有局部极值点，最值一定在4个顶点中，即，
P = { l x l y , l x u y , u x l y , u x u y } z ∈ [ min ⁡ { P } , max ⁡ { P } ] P=\{l_xl_y,l_xu_y,u_xl_y,u_xu_y\}\\ z\in[\min\{P\},\max\{P\}] P={lx​ly​,lx​uy​,ux​ly​,ux​uy​}z∈[min{P},max{P}]

除

3D图：

总体来说和乘法类似，但在除数为零时趋于无穷。这样的除法会产生开区间，突破从闭区间到闭区间封闭运算的语法前提，同时有“不确定性无穷大”这个语义上的困难。

将除法转化为乘法可以让我们看得更清晰一些：
$$z=\frac{x}{y}=x\frac{1}{y}$$
问题转化为，已知区间$y$，求$\frac{1}{y}$.

绘制$z=\frac{1}{y}$，

图中用红、橙、绿表示出三个区间，分别属于正区间、跨零区间、负区间。对于正区间和负区间都有，
$$z\in [\frac{1}{u},\frac{1}{l}]$$
但对于跨零区间这个结果已经不能代表不确定性，程序应在用户除以跨零区间时发出警告。

Scheme 实现

#lang sicp

(define (make-interval a b)
  (cons a b))

(define (upper-bound interval)
  (cdr interval))

(define (lower-bound interval)
  (car interval))


(define (add-interval x y)
  (make-interval (+ (lower-bound x) (lower-bound y))
                 (+ (upper-bound x) (upper-bound y))))

(define (mul-interval x y)
  (let ((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))))
    (make-interval (min p1 p2 p3 p4)
                   (max p1 p2 p3 p4))))

(define (div-interval x y)
  (mul-interval x
                (make-interval (/ 1.0 (upper-bound y))
                               (/ 1.0 (lower-bound y)))))

(define (sub-interval x y)
  (make-interval (- (lower-bound x) (upper-bound y))
                 (- (upper-bound x) (lower-bound y))))


; 对百分比表示的支持
(define (make-center-width c w)
  (make-interval (- c w) (+ c w)))

(define (center i)
  (/ (+ (lower-bound i) (upper-bound i)) 2))

(define (width i)
  (/ (- (upper-bound i) (lower-bound i)) 2))


(define (make-center-percent c p)
  (make-center-width c (* c p)))

(define (percent i)
  (/ (width i) (center i)))

潜在问题

使用本文的算法，在代数上等价的公式可能算出不同的结果！

例如电阻并联公式，
$$R=\frac{R_1{R}_2}{R_1+R_2}=\frac{1}{\frac{1}{R_1}+\frac{1}{R_2}}$$
对于同样不确定的$R_1$, $R_2$右式得出区间比左式窄。原因是左式中不确定量$R_1$和$R_2$都有重复出现，让同一份不确定性被计算了两次。一个更极端的例子能更清晰地展示这种现象，对区间$x$计算
$$\frac{x}{x}$$
从语义上说，分子与分母代表同一不确定量，其商应精确地为1，但本文的算法没有考虑同名变量相等：

> (let ((x (make-interval 5 10)))
    (div-interval x x))
(0.5 . 2.0)

但可以肯定的是，本算法给出的区间一定包含考虑同名变量的算法给出的区间。

附录

下面的代码摘自实验用jupyter notebook，用于探索乘法时最值点出现的情况以及乘法相对误差的传播规律，未经整理，供参考。

python 草稿

#!/usr/bin/env python
# coding: utf-8

# In[10]:


from random import randint
from itertools import *
from collections import Counter


# In[4]:


def p_gen():
    a = randint(1, 100)
    b = randint(a, 100)
    return a, b

def n_gen():
    a, b = p_gen()
    return -b, -a

def z_gen():
    return randint(-100, 0), randint(0, 100)


# In[8]:


p_intervals = [p_gen() for i in range(100)]
n_intervals = [n_gen() for i in range(100)]
z_intervals = [z_gen() for i in range(100)]


# In[29]:


case_pz = []
for p_interval, z_interval in product(p_intervals, z_intervals):
    bounds = [p_bound * z_bound for p_bound, z_bound in product(p_interval, z_interval)]
    maxi = max(range(len(bounds)), key=lambda x: bounds[x])
    mini = min(range(len(bounds)), key=lambda x: bounds[x])
    cases.append(((p_interval, z_interval), mini, maxi))


# In[30]:


[case[0] for  case in cases if case[1:] == (0, 1)]


# In[31]:


cnt = Counter(case[1:] for case in cases)


# In[ ]:


case_pn = []
for p_interval, n_interval in product(p_intervals, z_intervals):
    bounds = [p_bound * z_bound for p_bound, z_bound in product(p_interval, z_interval)]
    maxi = max(range(len(bounds)), key=lambda x: bounds[x])
    mini = min(range(len(bounds)), key=lambda x: bounds[x])
    cases.append(((p_interval, z_interval), mini, maxi))


# In[39]:


def test(left_intervals, right_intervals):
    cases = []
    for interval1, interval2 in product(left_intervals, right_intervals):
        bounds = [bound1 * bound2 for bound1, bound2 in product(interval1, interval2)]
        maxi = max(range(len(bounds)), key=lambda x: bounds[x])
        mini = min(range(len(bounds)), key=lambda x: bounds[x])
        cases.append(((interval1, interval2), mini, maxi))
    return cases

def count(cases):
    return Counter(case[1:] for case in cases)


# In[42]:


count(test(p_intervals, z_intervals))


# In[43]:


count(test(p_intervals, n_intervals))


# In[44]:


case_pn = test(p_intervals, n_intervals)


# In[45]:


cnt_pn = count(case_pn)


# In[49]:


[case[0] for case in case_pn if case[1:] == (0, 0)]


# In[50]:


count(test(z_intervals, p_intervals))


# In[51]:


count(test(n_intervals, p_intervals))


# In[55]:


for intervals1, intervals2 in combinations_with_replacement((p_intervals, z_intervals, n_intervals), 2):
    print('-' * 30)
    print(count(test(intervals1, intervals2)))
    print(count(test(intervals2, intervals1)))
    print('-' * 30)


# In[57]:


len(list(combinations_with_replacement(range(3), 2)))


# In[58]:


def percent(interval):
    a, b = interval
    return (b - a) / (b + a)


# In[59]:


def mul(i1, i2):
    return i1[0] * i2[0], i1[1] * i2[1]


# In[63]:


errors = []
for i1, i2 in product(p_intervals, repeat=2):
    # print('%.4f  %.4f  %.4f --- %.4f' % (p1:=percent(i1), p2:=percent(i2), p1 + p2, percent(mul(i1, i2))))
    errors.append(percent(i1) + percent(i2) - percent(mul(i1, i2)))


# In[66]:


(sum(i ** 2 for i in errors) / len(errors)) ** .5

errors = []
for i1, i2 in product(p_intervals, repeat=2):
    # print('%.4f  %.4f  %.4f --- %.4f' % (p1:=percent(i1), p2:=percent(i2), p1 + p2, percent(mul(i1, i2))))
    errors.append(percent(i1) + percent(i2) - percent(mul(i1, i2)))


# In[66]:


(sum(i ** 2 for i in errors) / len(errors)) ** .5

可见，

• 乘法最值分布情况复杂，我们不如直接在4个顶点中选最值
• 按百分比度量区间长度时，$percen{t}_z\approx percen{t}_x+percen{t}_y$