王国利信息科学与技术学院中山大学
http://human-robot.sysu.edu.cn
Matlab计算与仿真技术第十六讲,数值计算 -I
http://human-robot.sysu.edu.cn/course
Matlab计算与仿真
数值计算 -I
- 误差第十六讲提纲
Matlab计算与仿真数值分析 -I
误差
- 小误差的自然引入
>>format long e % display lots of digits
>> 2.6 + 0.2
ans =
2.800000000000000e+00
>> ans + 0.2
ans =
3.000000000000000e+00
Matlab计算与仿真数值计算 (续 )
- 小误差的自然引入 (续 )
>> ans + 0.2
ans =
3.200000000000001e+00
(为什么会出现最后一位非零数字呢?)
>> 2.6 + 0.6
ans =
3.200000000000000e+00
Matlab计算与仿真数值计算 (续 )
整数的存储
base 10 conversion base 2 (8bits)
1 1 = 2^0 0000 0001
2 2 = 2^1 0000 0010
4 4 = 2^2 0000 0100
8 8 = 2^3 0000 1000
9 8 + 1 = 2^3 + 2^0 0000 1001
10 8 + 2 = 2^3 + 2^1 0000 1010
(bbbbb...b)_2=b^(n-1)+b^(n-2)+...+b^0
(u)int8/int16/int32 指定存储方式
Matlab计算与仿真数值计算 (续 )
浮点数的存储
12.2792 =(符号 ) 0.123792(尾数 )x10^2(幂次 )
精 度 尾数比特数 幂次比特数
single 23 8
double 53 11
b(sign)bb...bbb(52 bits)bbbbbbbbbbb(11 bits)
尾数的比特数决定精度幂次的比特数决定范围
Matlab计算与仿真数值计算 (续 )
尾数可以表示成 2-k
2^(-0)=1,2^(-1)=0.5,2^(-2)=0.25,2^(-3)=0.125...
尾数的二进制转换算法
r0= x
for k=1,2,...,m
if rk-1 >= 2^(-k)
bk = 1
rk = rk-1-2^(-k)
else
bk = 0
rk = rk-1
endif
endfor
Matlab计算与仿真数值计算 (续 )
转换 算法执行示例输入,x=0.8215
k 2-k bk rk = rk-1- bk 2-k
0 NA NA 0.8125
1 0.5 1 0.3125
2 0.25 1 0.0625
3 0.125 0 0.0625
4 0.0625 1 0.0000
输出,x=(1101)2
Matlab计算与仿真数值计算 (续 )
输入,x=0.1
k 2-k bk rk = rk-1- bk 2-k
0 NA NA 0.1
1 0.5 0 0.1
2 0.25 0 0.1
3 0.125 0 0.1
4 0.0625 1 0.0375
5 0.03125 1 0.00625
6 0.015625 0 0.00625
7 0.0078125 0 0.00625
8 0.00390625 1 0.00234375
9 0.001953125 1 0.000390625
10 0.0009765625 0 0.000390625
....
Matlab计算与仿真数值计算 (续 )
输出,x=(00011 0011,,,)2,但是
0.1无法用有限比特数表示出来
Matlab计算与仿真数值计算 (续 )
Matlab浮点运算
>> format long e
>> u = 29/13
u =
2.230769230769231e+00
>> v = 13*u
v =
29
Matlab计算与仿真数值计算 (续 )
>> v-29
ans =
0
发生了两次误差取舍
>> x = 29/1300
x =
2.230769230769231e-02
>> y = 29 - 1300*x
y =
3.552713678800501e-015
Matlab计算与仿真数值计算 (续 )
- 二次方程的误差取舍解的显式表示考虑方程根分别为
2 0a x b x c
2 4
2
b b a c
x
a
2 5 4,3 2 0,1 0xx
1
2
5 4,3 2 1 8 1 5 8 9 9 5
0,0 0 1 8 4 1 0 0 4 9 5 7 6
x
x
Matlab计算与仿真数值计算 (续 )
注意到四位有效数字取舍计算
2b = 2 9 5 0,7 4 a c = 0,4
22
b - 4 a c = ( 5 4,3 2 ) - 0,4 0 0 0
= 2 9 5 1 - 0,4 0 0 0
= 2 9 5 1
= 5 4,3 2
Matlab计算与仿真数值计算 (续 )
从而正确解为,
2
1,4
4
2
+54,32 + 5 4.32
2.000
108.6
=
2.000
= 5 4.30
b b ac
x
a
1 5 4,3 2 1 8 1 5 8 9 9 5x?
约 0.4%的取舍累及误差
Matlab计算与仿真数值计算 (续 )
又正确解为
2
2,4
4
2
+ 54.3 2 - 54,32
2.000
0.000
=
2.000
= 0.0 00
b b ac
x
a
2 0,0 0 1 8 4 1 0 0 4 9 5 7 6x?
100%的取舍累及误差
Matlab计算与仿真数值计算 (续 )
解决方案,良态化表示求 解表达式
22
2
2
2
44
2 4
2
=
4
b b a c b b a c
x
a b b a c
c
b b a c
Matlab计算与仿真数值计算 (续 )
利用良态的公式计算注意到累积取舍误差限制在 0.05%内
2,4
0,2 0 0 0 0,2 0 0 0 0,0 0 1 8 4 2
5 4,3 2 5 4,3 2 1 0 8,6x
2 0,0 0 1 8 4 1 0 0 4 9 5 7 6x?
Matlab计算与仿真数值计算 (续 )
若同样然而利用其求解
22
1
2
2
44
2 4
2
=
4
b b a c b b a c
x
a b b a c
c
b b a c
1,4
0.20 00 0.20 00
54.3 2 54.3 2 0.00 0x
Matlab计算与仿真数值计算 (续 )
问题,
计算取舍产生灾变差稳健的求解方法,
引入
2 4b a c?
21q b + si gn(b) b 4
2 ac
12,
qcxx
aq
Matlab计算与仿真数值计算 (续 )
- 总结,
有限比特位限制精度产生截取误差截取误差微小,但累积截取误差 …
差异微小的数据相减差异巨大的数据相加都会影响精度截取误差客观存在,不可避免但算法设计可以弱化对精度负面影响
Matlab计算与仿真数值分析 (续 )
- 灾变差的一般概述灾变差产生的一般形式或其中 或考虑
ab
c a b -c a b?
ab?
c a b
0
8
.,.,1 0
.,.,1 0
a x x x x x x x
b y y y y y y y?
Matlab计算与仿真数值计算 (续 )
限 16-比特位运算结果
a x.xxx xxxx xxxx xxxx
b 0.000 0000 yyyy yyyy yyyy yyyy
c x.xxx xxxx zzzz zzzz (将被忽略 )
a 的有效位均被保留
b 的一半有效位丢失不难看出这是由两者巨大差异引起的
Matlab计算与仿真数值计算 (续 )
同样考虑双精度浮点运算结果
a x.xxx xxxx xxxx 1
b x.xxx xxxx xxxx 0
c 1.uuu uuuu uuuu X 10-12
u 为强制清零的比特位
c a b
.1
.0
a x x x x x x x x x x x x sssss s
b x x x x x x x x x x x x tt tt tt
Matlab计算与仿真数值计算 (续 )
误差度量与比较绝对差异,
|x-y|
相对差异,
|x-y|/|x|
Matlab计算与仿真数值计算 (续 )
- 浮点数的比较避免使用等量比较,例如
if x=y
…
尽可能采用差异容忍机制
if |x-y| < tol
tol 是设定的容忍误差上界
Matlab计算与仿真数值计算 (续 )
- 误差分析绝对误差,
Eabs (? ) =|? -A|
相对误差,
Erel (? ) |? -A|/|A|
其中? 是 A 的估计或近似
Matlab计算与仿真数值计算 (续 )
级数逼近分析
sin(x) 的级数展开形式
sin(x)=x-x3/3!+x5/5!+…
一阶逼近形式
sin(x)≈x,|x|<1
逼近误差为
Eabs=|x-sin(x)|=|x3/3!-x5/5!+…|
Erel=|x-sin(x)|/|sin(x)|
Matlab计算与仿真数值计算 (续 )
绝对误差 vs相对误差
Matlab计算与仿真数值计算 (续 )
- 迭代收敛的停机控制
{xk} 是 ξ 的 数值逼近序列对任意给定的误差上界 ε,存在 k
|xn-ξ|<ε,n>k
等价的时态差分形式
if |xn+1-xn| < ε’,n>k
ε’ 是时态误差上界
Matlab计算与仿真数值计算 (续 )
采用绝对误差停机控制
|x-xold| < tol
Matlab 程序模板
x =,.,% initialize
xold =,..
while abs(x-xold) > tol
xold = x;
update x
end
Matlab计算与仿真数值计算 (续 )
采用相对误差停机控制
|x-xold|/|xold| < tol
Matlab 程序模板
x =,.,% initialize
xold =,..
while abs((x-xold)/xold) > tol
xold = x;
update x
end
Matlab计算与仿真数值计算 (续 )
实例分析求解方程
cos(x) = x
算法,
1 猜测初值 x0
2 置 xold=x0
3 x=cos(xold)
4 if x=xold 停机,否则 xold=x 转 3
Matlab计算与仿真数值计算 (续 )
Matlab实现代码
x0 =,.,% initial guess
k = 0;
xnew = x0;
while NOT_CONVERGED & k < maxit
xold = xnew;
xnew = cos(xold);
k = k + 1;
end
Matlab计算与仿真数值计算 (续 )
NOT_CONVERGED 的设计不可取的方案
xnew=xold
绝对误差方案
|xnew-xold|>tol
相对误差方案
|(xnew-xold)/xold|/>tol
Matlab计算与仿真数值计算 (续 )
绝对误差方案
|xnew-xold|>tol
>> xold = 100;
>> xnew = 1;
>> delta = 5e-9;
>> (xnew-xold) > delta
ans =
0
Matlab计算与仿真第十六讲预告：数值计算 -II
（ 2008年 06月 18日）
结束语
http://human-robot.sysu.edu.cn
Matlab计算与仿真技术第十六讲,数值计算 -I
http://human-robot.sysu.edu.cn/course
Matlab计算与仿真
数值计算 -I
- 误差第十六讲提纲
Matlab计算与仿真数值分析 -I
误差
- 小误差的自然引入
>>format long e % display lots of digits
>> 2.6 + 0.2
ans =
2.800000000000000e+00
>> ans + 0.2
ans =
3.000000000000000e+00
Matlab计算与仿真数值计算 (续 )
- 小误差的自然引入 (续 )
>> ans + 0.2
ans =
3.200000000000001e+00
(为什么会出现最后一位非零数字呢?)
>> 2.6 + 0.6
ans =
3.200000000000000e+00
Matlab计算与仿真数值计算 (续 )
整数的存储
base 10 conversion base 2 (8bits)
1 1 = 2^0 0000 0001
2 2 = 2^1 0000 0010
4 4 = 2^2 0000 0100
8 8 = 2^3 0000 1000
9 8 + 1 = 2^3 + 2^0 0000 1001
10 8 + 2 = 2^3 + 2^1 0000 1010
(bbbbb...b)_2=b^(n-1)+b^(n-2)+...+b^0
(u)int8/int16/int32 指定存储方式
Matlab计算与仿真数值计算 (续 )
浮点数的存储
12.2792 =(符号 ) 0.123792(尾数 )x10^2(幂次 )
精 度 尾数比特数 幂次比特数
single 23 8
double 53 11
b(sign)bb...bbb(52 bits)bbbbbbbbbbb(11 bits)
尾数的比特数决定精度幂次的比特数决定范围
Matlab计算与仿真数值计算 (续 )
尾数可以表示成 2-k
2^(-0)=1,2^(-1)=0.5,2^(-2)=0.25,2^(-3)=0.125...
尾数的二进制转换算法
r0= x
for k=1,2,...,m
if rk-1 >= 2^(-k)
bk = 1
rk = rk-1-2^(-k)
else
bk = 0
rk = rk-1
endif
endfor
Matlab计算与仿真数值计算 (续 )
转换 算法执行示例输入,x=0.8215
k 2-k bk rk = rk-1- bk 2-k
0 NA NA 0.8125
1 0.5 1 0.3125
2 0.25 1 0.0625
3 0.125 0 0.0625
4 0.0625 1 0.0000
输出,x=(1101)2
Matlab计算与仿真数值计算 (续 )
输入,x=0.1
k 2-k bk rk = rk-1- bk 2-k
0 NA NA 0.1
1 0.5 0 0.1
2 0.25 0 0.1
3 0.125 0 0.1
4 0.0625 1 0.0375
5 0.03125 1 0.00625
6 0.015625 0 0.00625
7 0.0078125 0 0.00625
8 0.00390625 1 0.00234375
9 0.001953125 1 0.000390625
10 0.0009765625 0 0.000390625
....
Matlab计算与仿真数值计算 (续 )
输出,x=(00011 0011,,,)2,但是
0.1无法用有限比特数表示出来
Matlab计算与仿真数值计算 (续 )
Matlab浮点运算
>> format long e
>> u = 29/13
u =
2.230769230769231e+00
>> v = 13*u
v =
29
Matlab计算与仿真数值计算 (续 )
>> v-29
ans =
0
发生了两次误差取舍
>> x = 29/1300
x =
2.230769230769231e-02
>> y = 29 - 1300*x
y =
3.552713678800501e-015
Matlab计算与仿真数值计算 (续 )
- 二次方程的误差取舍解的显式表示考虑方程根分别为
2 0a x b x c
2 4
2
b b a c
x
a
2 5 4,3 2 0,1 0xx
1
2
5 4,3 2 1 8 1 5 8 9 9 5
0,0 0 1 8 4 1 0 0 4 9 5 7 6
x
x
Matlab计算与仿真数值计算 (续 )
注意到四位有效数字取舍计算
2b = 2 9 5 0,7 4 a c = 0,4
22
b - 4 a c = ( 5 4,3 2 ) - 0,4 0 0 0
= 2 9 5 1 - 0,4 0 0 0
= 2 9 5 1
= 5 4,3 2
Matlab计算与仿真数值计算 (续 )
从而正确解为,
2
1,4
4
2
+54,32 + 5 4.32
2.000
108.6
=
2.000
= 5 4.30
b b ac
x
a
1 5 4,3 2 1 8 1 5 8 9 9 5x?
约 0.4%的取舍累及误差
Matlab计算与仿真数值计算 (续 )
又正确解为
2
2,4
4
2
+ 54.3 2 - 54,32
2.000
0.000
=
2.000
= 0.0 00
b b ac
x
a
2 0,0 0 1 8 4 1 0 0 4 9 5 7 6x?
100%的取舍累及误差
Matlab计算与仿真数值计算 (续 )
解决方案,良态化表示求 解表达式
22
2
2
2
44
2 4
2
=
4
b b a c b b a c
x
a b b a c
c
b b a c
Matlab计算与仿真数值计算 (续 )
利用良态的公式计算注意到累积取舍误差限制在 0.05%内
2,4
0,2 0 0 0 0,2 0 0 0 0,0 0 1 8 4 2
5 4,3 2 5 4,3 2 1 0 8,6x
2 0,0 0 1 8 4 1 0 0 4 9 5 7 6x?
Matlab计算与仿真数值计算 (续 )
若同样然而利用其求解
22
1
2
2
44
2 4
2
=
4
b b a c b b a c
x
a b b a c
c
b b a c
1,4
0.20 00 0.20 00
54.3 2 54.3 2 0.00 0x
Matlab计算与仿真数值计算 (续 )
问题,
计算取舍产生灾变差稳健的求解方法,
引入
2 4b a c?
21q b + si gn(b) b 4
2 ac
12,
qcxx
aq
Matlab计算与仿真数值计算 (续 )
- 总结,
有限比特位限制精度产生截取误差截取误差微小,但累积截取误差 …
差异微小的数据相减差异巨大的数据相加都会影响精度截取误差客观存在,不可避免但算法设计可以弱化对精度负面影响
Matlab计算与仿真数值分析 (续 )
- 灾变差的一般概述灾变差产生的一般形式或其中 或考虑
ab
c a b -c a b?
ab?
c a b
0
8
.,.,1 0
.,.,1 0
a x x x x x x x
b y y y y y y y?
Matlab计算与仿真数值计算 (续 )
限 16-比特位运算结果
a x.xxx xxxx xxxx xxxx
b 0.000 0000 yyyy yyyy yyyy yyyy
c x.xxx xxxx zzzz zzzz (将被忽略 )
a 的有效位均被保留
b 的一半有效位丢失不难看出这是由两者巨大差异引起的
Matlab计算与仿真数值计算 (续 )
同样考虑双精度浮点运算结果
a x.xxx xxxx xxxx 1
b x.xxx xxxx xxxx 0
c 1.uuu uuuu uuuu X 10-12
u 为强制清零的比特位
c a b
.1
.0
a x x x x x x x x x x x x sssss s
b x x x x x x x x x x x x tt tt tt
Matlab计算与仿真数值计算 (续 )
误差度量与比较绝对差异,
|x-y|
相对差异,
|x-y|/|x|
Matlab计算与仿真数值计算 (续 )
- 浮点数的比较避免使用等量比较,例如
if x=y
…
尽可能采用差异容忍机制
if |x-y| < tol
tol 是设定的容忍误差上界
Matlab计算与仿真数值计算 (续 )
- 误差分析绝对误差,
Eabs (? ) =|? -A|
相对误差,
Erel (? ) |? -A|/|A|
其中? 是 A 的估计或近似
Matlab计算与仿真数值计算 (续 )
级数逼近分析
sin(x) 的级数展开形式
sin(x)=x-x3/3!+x5/5!+…
一阶逼近形式
sin(x)≈x,|x|<1
逼近误差为
Eabs=|x-sin(x)|=|x3/3!-x5/5!+…|
Erel=|x-sin(x)|/|sin(x)|
Matlab计算与仿真数值计算 (续 )
绝对误差 vs相对误差
Matlab计算与仿真数值计算 (续 )
- 迭代收敛的停机控制
{xk} 是 ξ 的 数值逼近序列对任意给定的误差上界 ε,存在 k
|xn-ξ|<ε,n>k
等价的时态差分形式
if |xn+1-xn| < ε’,n>k
ε’ 是时态误差上界
Matlab计算与仿真数值计算 (续 )
采用绝对误差停机控制
|x-xold| < tol
Matlab 程序模板
x =,.,% initialize
xold =,..
while abs(x-xold) > tol
xold = x;
update x
end
Matlab计算与仿真数值计算 (续 )
采用相对误差停机控制
|x-xold|/|xold| < tol
Matlab 程序模板
x =,.,% initialize
xold =,..
while abs((x-xold)/xold) > tol
xold = x;
update x
end
Matlab计算与仿真数值计算 (续 )
实例分析求解方程
cos(x) = x
算法,
1 猜测初值 x0
2 置 xold=x0
3 x=cos(xold)
4 if x=xold 停机,否则 xold=x 转 3
Matlab计算与仿真数值计算 (续 )
Matlab实现代码
x0 =,.,% initial guess
k = 0;
xnew = x0;
while NOT_CONVERGED & k < maxit
xold = xnew;
xnew = cos(xold);
k = k + 1;
end
Matlab计算与仿真数值计算 (续 )
NOT_CONVERGED 的设计不可取的方案
xnew=xold
绝对误差方案
|xnew-xold|>tol
相对误差方案
|(xnew-xold)/xold|/>tol
Matlab计算与仿真数值计算 (续 )
绝对误差方案
|xnew-xold|>tol
>> xold = 100;
>> xnew = 1;
>> delta = 5e-9;
>> (xnew-xold) > delta
ans =
0
Matlab计算与仿真第十六讲预告：数值计算 -II
（ 2008年 06月 18日）
结束语
