Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002年
http://www.stanford.edu/class/cs138/mlmath.html
黄金分割
1
11 ?? ?
?
1
11 ?? ?
? 012 ??? ??
2
51 ???
p = [1 –1 –1];
r = roots(p);
1)( 2 ??? xxxp
r =
-0.61803398874989
1.61803398874989
C(1)*X^N +,.,+ C(N)*X + C(N+1) = 0
r = solve(‘1/x = x-1’); r =
[ 1/2*5^(1/2)+1/2]
[ 1/2-1/2*5^(1/2)]
符号工具箱
Mapple
phi = r(1)
vpa (phi,50)
1/2*5^(1/2)+1/2
1.6180339887498948482045868343656381177203091798058
phi = double(phi)
1.61803398874989
图形的方法
f = inline('1/x-(x-1)')
ezplot(f,0,4)
phi = fzero(f,1)
hold on
plot(phi,0,’o’)
0 0.5 1 1.5 2 2.5 3 3.5 4
-3
-2
-1
0
1
2
3
4
5
6
7
x
1/x-(x-1)
结果图
?
? - 1
1
1
黄金分割图
%GOLDRECT Golden Rectangle
% GOLDRECT plots the golden rectangle
phi = (1+sqrt(5))/2; x = [0 phi phi 0 0];
y = [0 0 1 1 0]; u = [1 1]; v = [0 1];
plot(x,y,'b',u,v,'b--'); text(phi/2,1.05,'\phi')
text((1+phi)/2,-.05,'\phi - 1');
text(-.05,.5,'1'); text(.5,-.05,'1')
axis equal
axis off
set(gcf,'color','white')
连分式
...
1
1
1
3
2
1
0
?
?
?
?
a
a
a
a
.,,1
1
1
1
1
1
1
?
?
?
???
连分式
Goldfract.m
g =
1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))
g =
21/13
g =
1.61538461538462
err =
0.0026
字符串
有效数字 2位
Fibonacci
Leonardo Pisano Fibonacci was born around 1170 and died
around 1250 in Pisa in what is now Italy,He traveled extensively
in Europe and Northern Africa,He wrote several mathematical
texts that,among other things,introduced Europe to the Hindu-
Arabic notation for numbers,Even though his books had to be
transcribed by hand,they were still widely circulated,In his best
known book,Liber Abaci,published in 1202,he posed the
following problem.
A man put a pair of rabbits in a place surrounded on all sides by
a
wall,How many pairs of rabbits can be produced from that pair in
a
year if it is supposed that every month each pair begets a new
兔子的故事
兔子从出生要一个月，从出生到成熟要过一个月
月 1 2 3 4 5 6
繁殖 1 1 1 2 3 5
出生 0 1 1 2 3 5
成熟 0 0 1 1 2 3
总数 1 2 3 5 8 13
Fibonacci序列
21 ?? ?? nnn fff
2,1 20 ?? ff
Fibonacci.m
1
2
3
5
8
13
21
34
55
89
144
233
递归实现：
Fibnum.m
递归实现
递归实现是 elegant but expensive
tic,fibnum(24),toc
tic,fibonacci(24),toc
试一下，比较时间！
黄金分割和 Fibonacci数
比较一下 goldfract(6)和 fibonacci(7)
连分式,g = g + f;
黄金分割,f(k) = f(k-1) +f(k-2);
n
n
n
f
f ??? 1 ???
?? n
n
n f
f 1lim
n = 40;
f = fibonacci(n);
f(2:n)./f(1:n-1)
f(2:n)./f(1:n-1) - phi
Fibonacci的兔子
Fibonacci的兔子以黄金分割的速度增长。
nn cf ??
21 ?? ?? nnn fff
＋
12 ?? ??
nnn ccf )1(21 ?? ???有两个解 phi和 (1-phi)
由初始条件：
12
)1(,
12 21 ?
???
?? ?
?
?
? cc
))1((12 1 11 ?? ???? nnnf ??? 是有理分式
Fibonacci的兔子
注意：没有半只兔子 ?
format long e
n = (1:40)';
f = (phi.^(n+1)-(1-phi).^(n+1))/(2*phi-1)
f = round(f)
魔方阵
A = magic(3)
sum(A)
sum(A’)
sum(diag(A))
sum(diag(flipud(A)))
sum(1:9)/3
8 1 6
3 5 7
4 9 2
魔方阵的八种组合
for k=0:3
rot90(A,k)
rot90(A’,k)
end
逆时阵旋转
k*90度
基本代数运算
X = inv(A)
det(A) -360
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
format rat 53/360 -13/90 23/360
-11/180 1/45 19/180
-7/360 17/90 -37/360
基本代数运算
r = norm(A)
e = eig(A)
s = svd(A)
e =
15.0000
4.8990
-4.8990
s =
15.0000
6.9282
3.4641
r =
15
NORM(X) is the largest singular value of X,max(svd(X)).
符号计算
A = sym(A)
sum(A)
sum(A’)’
det(A)
inv(A)
eig(A)
svd(A)
100 200 300 400 500
100
200
300
400
500
600
图象显示
load durer
whos
image(X)
colormap(map)
axis image
Melancolia,a
Renaissance etching
by Albrect Durer
50 100 150 200 250 300 350
50
100
150
200
250
300
350
4阶魔方阵
load detail
image(X)
colormap(map)
axis image
4阶魔方阵
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
A = magic(4)
A = A(:,[1 3 2 4])
4阶魔方阵,880个
5阶魔方阵,275305224个
6阶魔方阵：？
秩和奇异阵
A = magic(4)
det(A)
inv(A)
矩阵的秩是矩阵中线性无关的行和列的个数，
一个 n阶矩阵是奇异的当且仅当它的秩小于 n
rank(magic(4))
秩和奇异阵
for n = 1:24,r(n) = rank(magic(n)); end
[(1:24)’ r’]
bar(r)
tiltle(‘魔方阵的秩’ )
0 5 10 15 20 250
5
10
15
20
25
n,奇数，满秩
n=4k:秩是 3
n=2k !=4k:秩是 n/2+2
魔方阵的三维显示
surf(magic(n))
axis off
set(gcf,’doublebuffer’,’on’)
cameratoolbar
for n = 8:11
subplot(2,2,n-7),surf(magic(n))
axis off,view(30,45)
%axis tight
end
图形对象
小鸭子
加 密－－ Matlab字符能力
x = reshape(32:127,32,3)’
c = char(x)
char(32)
c =
!"#$%&'()*+,-./0123456789:;<=>?
@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_
`abcdefghijklmnopqrstuvwxyz{|}~
字符串
d = char(48:57)
double(d) –‘0’
d =
0123456789
ans =
0 1 2 3 4 5 6 7 8 9
字符串（高八位）
char(reshape(160:255,32,3)’)
char(160)
char(160:161)
牎 ⅲ ぅ Ηī 氨渤吹斗腹夯冀究
懒旅呐魄壬仕掏蜗醒矣哉肿刭谯茌捱
模运算
x = [37 –37 37 –37]’;
y = [10 10 -10 -10]’;
r = [x y rem(x,y) mod(x,y)]
37 10 7 7
-37 10 -7 3
37 -10 7 -3
-37 -10 -7 -7
字符加密
取素数,P = 97 (要表示 95个字符）
x = double(‘TV’ )’-32
y = Ax mod P
A = [ 71 2
2 26]
‘TV’?’1U’
If y = Ax mod P
then x = Ay mod P mod(A^2,P) = I
Crypto.m ‘ hello world’?,?{p]>~Y&=
分形蕨
Fern.m
数论中未解决的 3n+1问题
如果 n = 1,停止
如果 n是偶数,n = n/2;
如果 n是奇数,n = 3n+1
问题：这个过程是否一定会 终止？
综述文章,American Mathematical Monthly,92(1985),3-23
http://www.cecm.sfu.ca/organics/papers/lagarias
threenplus1.m
function threenplus1(n)
%"Three n plus 1".
% Study the 3n+1 sequence.
% threenplus1(n) plots the sequence starting with n.
% threenplus1 with no arguments starts with n = 1.
% uicontrols decrement or increment the starting n.
% Is it possible for this to run forever?
if ~isequal(get(gcf,'tag'),'3n+1')
shg
clf reset
uicontrol('position',[260 5 25 22],...
'string','<',...
'callback','threenplus1(''<'')');
uicontrol('position',[300 5 25 22],...
'string','>',...
'callback','threenplus1(''>'')');
set(gcf,'tag','3n+1');
end
if nargin == 0
n = 1;
elseif isequal(n,'<')
n = get(gcf,'userdata') - 1;
elseif isequal(n,'>')
n = get(gcf,'userdata') + 1;
end
if n < 1,n = 1; end
set(gcf,'userdata',n)
y = n;
while n > 1
if rem(n,2)==0
n = n/2;
else
n = 3*n+1;
end
y = [y n];
end
semilogy(y,'.-')
axis tight
ymax = max(y);
ytick = [2.^(0:ceil(log2(ymax))-1) ymax];
if length(ytick) > 8,ytick(end-1) = []; end
set(gca,'ytick',ytick)
title(['n = ' num2str(y(1))]);
x = 32768
y = 0
L,lcad y
shift right 5 bits
add x
store in x
change sign
shift right 5 bits
add y
store in y
plot x y
go to L
画 园
h = 1/32;
x = 1;
y = 0;
while 1
x = x + h*y;
y = y - h*x;
plot(x,y,’.’)
drawnow
end
Circlegen.m
画 园
11
1
??
?
??
??
nnn
nnn
hxyy
hyxx
画园中的特征值
11
1
??
?
??
??
nnn
nnn
hxyy
hyxx
nnn
nnn
yhhxy
hyxx
)1( 21
1
????
??
?
?
???
?
???
?
??? 21
1
hh
hA 01 xAAxx nnn ???
[V,E] = eig(A)
1???? V E VAVEAV 1?? VVEA nn
E是对角阵
有界nk A?? 1?
h = 2*rand,A=[1 h; -h 1-h^2],lambda=eig(A),abs(lambda)
特征值的符号运算
syms h
A = [1 h; -h 1-h^2]
lambda = eig(A)
d = det(A)
d = simple(prod(lambda))
??? s inc o s i??
4/1s i n,2/1c o s 22 hhh ???? ?? 4/12/1 22 hhh ??????
设
theta = acos(1-h^2/2） ;
Lambda=[cos(theta)- i*sin(theta);
cos(theta)+i*sin(theta)]
diff = simple(lambda-Lambda)
浮点数
1985年,IEEE标准委员会和美国国家标准局对二进制
浮点数运算建立 ANSI/IEEE 745-1985标准。
Matlab采用 IEEE双精度浮点运算。
10,2)1( ????? ffx e
f最多用 52位表示 是整数ff 525252 2,220 ??
规范数表示
精度与 f有关,表示的数的范围与 e有关
10231022 ??? e指数 e是整数
双精度浮点数
双精度浮点数以 64位存储,f:52位,e:11位,数的符号,1位。
指数以 e+1023存储 (1,211-2)
整个浮点数的分数部分不是 f,而是 1+f,IEEE标准在 64位
中存储了 65位的信息。
Floatgui.m
eps
精度和舍入误差
eps是 1到一个最近浮点数的距离,eps=2^(-t)
对 IEEE双精度,eps = 2^(-52) ≈2.2204*10^(-16)
最大的相对误差是 eps/2,数之间的最大的相对距
离是 eps
t = 0.1不能被被精确存储
??????????? 121110987654 2 12 02 02 1212 02 0212 110 1
第 1项的 1不用存储，所以重复出现 1001，知道第
53项，而 54项以后的级数大于 eps/2,舍入到 53项
4
131232 216
9
16
9
16
9
16
9
16
91 ???
?
??
?
? ??????? ?t
精度和舍入误差
t = 0.1
format hex
t
3fb999999999999a
3fbh = 1019
1019-1023 = - 4试一下,-0.1是如何存的？
看一下,0.3/0.1和 3的存储
0:0.1:1 的问题
format long
x = 4/3 –1
y = 3*x
z = 1 - y
?
浮点数
Binary Decimal Hex
eps 2^(-52) 2.2204e-16 3cb0000000000000
realmin 2^(-1022) 2.2251e-308 0010000000000000
realmax (2-eps)*2^1023 1.7977e+308 7fefffffffffffff
上溢，下溢
inf 无穷大, e = 1024,f = 0; 7ff0000000000000
nan 不是数,e=1024,f != 0; fff8000000000000
subnormal浮点数
许多机器支持 subnormal浮点数,realmin ~ eps*realmin
标志,e = -1023,例如 eps*realmin ＝ 0000000000000001
Matlab也用浮点数表示整数，我们用 flint来指这些整数，
在 flint浮点运算不引入舍入误差，加减乘的结果还是 flint，
当除和平方根的结果是整数时，也产生 flint
舍入误差和方程求解
2.25.07.1
22517
21
21
??
??
xx
xx
A = [17 5; 1.7 0.5]
b = [22;2.2]
x = A\b
x =
-1.0588
8.0000
t = 1.7/17
A(2,:) = A(2,:) – t*A(1,:)
b(2) = b(2) – t*b(1)
x(2) = b(2) /A(2,2)
x(1) = (22 – 5*x(2))/17
0.985 0.99 0.995 1 1.005 1.01 1.015-5
-4
-3
-2
-1
0
1
2
3
4
5 x 10
-14
舍入误差和多项式
x = 0.988:.0001:1.012;
y = x.^7 – 7*x.^6+21*x.^5-35*x.^4+…
35*x.^3-21*x.^2+7*x-1;
plot(x,y)
-6 -4 -2 0 2 4 6
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
x 105
x
(x-1)7
ezplot
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002年
http://www.stanford.edu/class/cs138/mlmath.html
黄金分割
1
11 ?? ?
?
1
11 ?? ?
? 012 ??? ??
2
51 ???
p = [1 –1 –1];
r = roots(p);
1)( 2 ??? xxxp
r =
-0.61803398874989
1.61803398874989
C(1)*X^N +,.,+ C(N)*X + C(N+1) = 0
r = solve(‘1/x = x-1’); r =
[ 1/2*5^(1/2)+1/2]
[ 1/2-1/2*5^(1/2)]
符号工具箱
Mapple
phi = r(1)
vpa (phi,50)
1/2*5^(1/2)+1/2
1.6180339887498948482045868343656381177203091798058
phi = double(phi)
1.61803398874989
图形的方法
f = inline('1/x-(x-1)')
ezplot(f,0,4)
phi = fzero(f,1)
hold on
plot(phi,0,’o’)
0 0.5 1 1.5 2 2.5 3 3.5 4
-3
-2
-1
0
1
2
3
4
5
6
7
x
1/x-(x-1)
结果图
?
? - 1
1
1
黄金分割图
%GOLDRECT Golden Rectangle
% GOLDRECT plots the golden rectangle
phi = (1+sqrt(5))/2; x = [0 phi phi 0 0];
y = [0 0 1 1 0]; u = [1 1]; v = [0 1];
plot(x,y,'b',u,v,'b--'); text(phi/2,1.05,'\phi')
text((1+phi)/2,-.05,'\phi - 1');
text(-.05,.5,'1'); text(.5,-.05,'1')
axis equal
axis off
set(gcf,'color','white')
连分式
...
1
1
1
3
2
1
0
?
?
?
?
a
a
a
a
.,,1
1
1
1
1
1
1
?
?
?
???
连分式
Goldfract.m
g =
1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))
g =
21/13
g =
1.61538461538462
err =
0.0026
字符串
有效数字 2位
Fibonacci
Leonardo Pisano Fibonacci was born around 1170 and died
around 1250 in Pisa in what is now Italy,He traveled extensively
in Europe and Northern Africa,He wrote several mathematical
texts that,among other things,introduced Europe to the Hindu-
Arabic notation for numbers,Even though his books had to be
transcribed by hand,they were still widely circulated,In his best
known book,Liber Abaci,published in 1202,he posed the
following problem.
A man put a pair of rabbits in a place surrounded on all sides by
a
wall,How many pairs of rabbits can be produced from that pair in
a
year if it is supposed that every month each pair begets a new
兔子的故事
兔子从出生要一个月，从出生到成熟要过一个月
月 1 2 3 4 5 6
繁殖 1 1 1 2 3 5
出生 0 1 1 2 3 5
成熟 0 0 1 1 2 3
总数 1 2 3 5 8 13
Fibonacci序列
21 ?? ?? nnn fff
2,1 20 ?? ff
Fibonacci.m
1
2
3
5
8
13
21
34
55
89
144
233
递归实现：
Fibnum.m
递归实现
递归实现是 elegant but expensive
tic,fibnum(24),toc
tic,fibonacci(24),toc
试一下，比较时间！
黄金分割和 Fibonacci数
比较一下 goldfract(6)和 fibonacci(7)
连分式,g = g + f;
黄金分割,f(k) = f(k-1) +f(k-2);
n
n
n
f
f ??? 1 ???
?? n
n
n f
f 1lim
n = 40;
f = fibonacci(n);
f(2:n)./f(1:n-1)
f(2:n)./f(1:n-1) - phi
Fibonacci的兔子
Fibonacci的兔子以黄金分割的速度增长。
nn cf ??
21 ?? ?? nnn fff
＋
12 ?? ??
nnn ccf )1(21 ?? ???有两个解 phi和 (1-phi)
由初始条件：
12
)1(,
12 21 ?
???
?? ?
?
?
? cc
))1((12 1 11 ?? ???? nnnf ??? 是有理分式
Fibonacci的兔子
注意：没有半只兔子 ?
format long e
n = (1:40)';
f = (phi.^(n+1)-(1-phi).^(n+1))/(2*phi-1)
f = round(f)
魔方阵
A = magic(3)
sum(A)
sum(A’)
sum(diag(A))
sum(diag(flipud(A)))
sum(1:9)/3
8 1 6
3 5 7
4 9 2
魔方阵的八种组合
for k=0:3
rot90(A,k)
rot90(A’,k)
end
逆时阵旋转
k*90度
基本代数运算
X = inv(A)
det(A) -360
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
format rat 53/360 -13/90 23/360
-11/180 1/45 19/180
-7/360 17/90 -37/360
基本代数运算
r = norm(A)
e = eig(A)
s = svd(A)
e =
15.0000
4.8990
-4.8990
s =
15.0000
6.9282
3.4641
r =
15
NORM(X) is the largest singular value of X,max(svd(X)).
符号计算
A = sym(A)
sum(A)
sum(A’)’
det(A)
inv(A)
eig(A)
svd(A)
100 200 300 400 500
100
200
300
400
500
600
图象显示
load durer
whos
image(X)
colormap(map)
axis image
Melancolia,a
Renaissance etching
by Albrect Durer
50 100 150 200 250 300 350
50
100
150
200
250
300
350
4阶魔方阵
load detail
image(X)
colormap(map)
axis image
4阶魔方阵
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
A = magic(4)
A = A(:,[1 3 2 4])
4阶魔方阵,880个
5阶魔方阵,275305224个
6阶魔方阵：？
秩和奇异阵
A = magic(4)
det(A)
inv(A)
矩阵的秩是矩阵中线性无关的行和列的个数，
一个 n阶矩阵是奇异的当且仅当它的秩小于 n
rank(magic(4))
秩和奇异阵
for n = 1:24,r(n) = rank(magic(n)); end
[(1:24)’ r’]
bar(r)
tiltle(‘魔方阵的秩’ )
0 5 10 15 20 250
5
10
15
20
25
n,奇数，满秩
n=4k:秩是 3
n=2k !=4k:秩是 n/2+2
魔方阵的三维显示
surf(magic(n))
axis off
set(gcf,’doublebuffer’,’on’)
cameratoolbar
for n = 8:11
subplot(2,2,n-7),surf(magic(n))
axis off,view(30,45)
%axis tight
end
图形对象
小鸭子
加 密－－ Matlab字符能力
x = reshape(32:127,32,3)’
c = char(x)
char(32)
c =
!"#$%&'()*+,-./0123456789:;<=>?
@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_
`abcdefghijklmnopqrstuvwxyz{|}~
字符串
d = char(48:57)
double(d) –‘0’
d =
0123456789
ans =
0 1 2 3 4 5 6 7 8 9
字符串（高八位）
char(reshape(160:255,32,3)’)
char(160)
char(160:161)
牎 ⅲ ぅ Ηī 氨渤吹斗腹夯冀究
懒旅呐魄壬仕掏蜗醒矣哉肿刭谯茌捱
模运算
x = [37 –37 37 –37]’;
y = [10 10 -10 -10]’;
r = [x y rem(x,y) mod(x,y)]
37 10 7 7
-37 10 -7 3
37 -10 7 -3
-37 -10 -7 -7
字符加密
取素数,P = 97 (要表示 95个字符）
x = double(‘TV’ )’-32
y = Ax mod P
A = [ 71 2
2 26]
‘TV’?’1U’
If y = Ax mod P
then x = Ay mod P mod(A^2,P) = I
Crypto.m ‘ hello world’?,?{p]>~Y&=
分形蕨
Fern.m
数论中未解决的 3n+1问题
如果 n = 1,停止
如果 n是偶数,n = n/2;
如果 n是奇数,n = 3n+1
问题：这个过程是否一定会 终止？
综述文章,American Mathematical Monthly,92(1985),3-23
http://www.cecm.sfu.ca/organics/papers/lagarias
threenplus1.m
function threenplus1(n)
%"Three n plus 1".
% Study the 3n+1 sequence.
% threenplus1(n) plots the sequence starting with n.
% threenplus1 with no arguments starts with n = 1.
% uicontrols decrement or increment the starting n.
% Is it possible for this to run forever?
if ~isequal(get(gcf,'tag'),'3n+1')
shg
clf reset
uicontrol('position',[260 5 25 22],...
'string','<',...
'callback','threenplus1(''<'')');
uicontrol('position',[300 5 25 22],...
'string','>',...
'callback','threenplus1(''>'')');
set(gcf,'tag','3n+1');
end
if nargin == 0
n = 1;
elseif isequal(n,'<')
n = get(gcf,'userdata') - 1;
elseif isequal(n,'>')
n = get(gcf,'userdata') + 1;
end
if n < 1,n = 1; end
set(gcf,'userdata',n)
y = n;
while n > 1
if rem(n,2)==0
n = n/2;
else
n = 3*n+1;
end
y = [y n];
end
semilogy(y,'.-')
axis tight
ymax = max(y);
ytick = [2.^(0:ceil(log2(ymax))-1) ymax];
if length(ytick) > 8,ytick(end-1) = []; end
set(gca,'ytick',ytick)
title(['n = ' num2str(y(1))]);
x = 32768
y = 0
L,lcad y
shift right 5 bits
add x
store in x
change sign
shift right 5 bits
add y
store in y
plot x y
go to L
画 园
h = 1/32;
x = 1;
y = 0;
while 1
x = x + h*y;
y = y - h*x;
plot(x,y,’.’)
drawnow
end
Circlegen.m
画 园
11
1
??
?
??
??
nnn
nnn
hxyy
hyxx
画园中的特征值
11
1
??
?
??
??
nnn
nnn
hxyy
hyxx
nnn
nnn
yhhxy
hyxx
)1( 21
1
????
??
?
?
???
?
???
?
??? 21
1
hh
hA 01 xAAxx nnn ???
[V,E] = eig(A)
1???? V E VAVEAV 1?? VVEA nn
E是对角阵
有界nk A?? 1?
h = 2*rand,A=[1 h; -h 1-h^2],lambda=eig(A),abs(lambda)
特征值的符号运算
syms h
A = [1 h; -h 1-h^2]
lambda = eig(A)
d = det(A)
d = simple(prod(lambda))
??? s inc o s i??
4/1s i n,2/1c o s 22 hhh ???? ?? 4/12/1 22 hhh ??????
设
theta = acos(1-h^2/2） ;
Lambda=[cos(theta)- i*sin(theta);
cos(theta)+i*sin(theta)]
diff = simple(lambda-Lambda)
浮点数
1985年,IEEE标准委员会和美国国家标准局对二进制
浮点数运算建立 ANSI/IEEE 745-1985标准。
Matlab采用 IEEE双精度浮点运算。
10,2)1( ????? ffx e
f最多用 52位表示 是整数ff 525252 2,220 ??
规范数表示
精度与 f有关,表示的数的范围与 e有关
10231022 ??? e指数 e是整数
双精度浮点数
双精度浮点数以 64位存储,f:52位,e:11位,数的符号,1位。
指数以 e+1023存储 (1,211-2)
整个浮点数的分数部分不是 f,而是 1+f,IEEE标准在 64位
中存储了 65位的信息。
Floatgui.m
eps
精度和舍入误差
eps是 1到一个最近浮点数的距离,eps=2^(-t)
对 IEEE双精度,eps = 2^(-52) ≈2.2204*10^(-16)
最大的相对误差是 eps/2,数之间的最大的相对距
离是 eps
t = 0.1不能被被精确存储
??????????? 121110987654 2 12 02 02 1212 02 0212 110 1
第 1项的 1不用存储，所以重复出现 1001，知道第
53项，而 54项以后的级数大于 eps/2,舍入到 53项
4
131232 216
9
16
9
16
9
16
9
16
91 ???
?
??
?
? ??????? ?t
精度和舍入误差
t = 0.1
format hex
t
3fb999999999999a
3fbh = 1019
1019-1023 = - 4试一下,-0.1是如何存的？
看一下,0.3/0.1和 3的存储
0:0.1:1 的问题
format long
x = 4/3 –1
y = 3*x
z = 1 - y
?
浮点数
Binary Decimal Hex
eps 2^(-52) 2.2204e-16 3cb0000000000000
realmin 2^(-1022) 2.2251e-308 0010000000000000
realmax (2-eps)*2^1023 1.7977e+308 7fefffffffffffff
上溢，下溢
inf 无穷大, e = 1024,f = 0; 7ff0000000000000
nan 不是数,e=1024,f != 0; fff8000000000000
subnormal浮点数
许多机器支持 subnormal浮点数,realmin ~ eps*realmin
标志,e = -1023,例如 eps*realmin ＝ 0000000000000001
Matlab也用浮点数表示整数，我们用 flint来指这些整数，
在 flint浮点运算不引入舍入误差，加减乘的结果还是 flint，
当除和平方根的结果是整数时，也产生 flint
舍入误差和方程求解
2.25.07.1
22517
21
21
??
??
xx
xx
A = [17 5; 1.7 0.5]
b = [22;2.2]
x = A\b
x =
-1.0588
8.0000
t = 1.7/17
A(2,:) = A(2,:) – t*A(1,:)
b(2) = b(2) – t*b(1)
x(2) = b(2) /A(2,2)
x(1) = (22 – 5*x(2))/17
0.985 0.99 0.995 1 1.005 1.01 1.015-5
-4
-3
-2
-1
0
1
2
3
4
5 x 10
-14
舍入误差和多项式
x = 0.988:.0001:1.012;
y = x.^7 – 7*x.^6+21*x.^5-35*x.^4+…
35*x.^3-21*x.^2+7*x-1;
plot(x,y)
-6 -4 -2 0 2 4 6
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
x 105
x
(x-1)7
ezplot
