区间算数

区间算数

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

文章目录

• 区间算数
• • 引言
  • 数学原理
  • • 加
    • 减
    • 乘
    • 除
  • Scheme 实现
  • 潜在问题
  • 附录
  • • python 草稿

引言

区间算数是指一套定义在区间上的运算。用区间表示不确定量的取值范围，常用来衡量测量值的误差，如 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 ∈ [ l x + l y , u x + u y ] z\in[l_x+l_y,u_x+u_y] z∈[lx​+ly​,ux​+uy​]
想象一个在$z$平面上滑动的矩形，可以想象到，上述结论对所有可能的区间都成立，与区间的相对位置、区间的正负均无关。

减

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

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

乘

乘法比较特殊，如图：

想象滑动的矩形，至少有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 = x y = x 1 y z=\frac{x}{y}=x\frac{1}{y} z=yx​=xy1​
问题转化为，已知区间$y$，求 1 y \frac{1}{y} y1​.

绘制 z = 1 y z=\frac{1}{y} z=y1​，

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

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 = R 1 R 2 R 1 + R 2 = 1 1 R 1 + 1 R 2 R=\frac{R_1R_2}{R_1+R_2}=\frac{1}{\frac{1}{R_1}+\frac{1}{R_2}} R=R1​+R2​R1​R2​​=R1​1​+R2​1​1​
对于同样不确定的 R 1 R_1 R1​, R 2 R_2 R2​右式得出区间比左式窄。原因是左式中不确定量 R 1 R_1 R1​和 R 2 R_2 R2​都有重复出现，让同一份不确定性被计算了两次。一个更极端的例子能更清晰地展示这种现象，对区间$x$计算
x x \frac{x}{x} xx​
从语义上说，分子与分母代表同一不确定量，其商应精确地为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, bdef n_gen():a, b = p_gen()return -b, -adef 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 casesdef 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 = intervalreturn (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)) ** .5errors = []
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个顶点中选最值
• 按百分比度量区间长度时， p e r c e n t z ≈ p e r c e n t x + p e r c e n t y percent_z\approx percent_x+percent_y percentz​≈percentx​+percenty​