Matlab Math Cleve Morler wbchenfudan edu cn 2002
Matlab Math Cleve Morler著 陈文斌(wbchen@fudan. edu. cn) 复旦大学 2002年
http: //www. stanford. edu/class/cs 138/mlmath. html
p = [1 – 1]; r = roots(p); r= -0. 61803398874989 1. 61803398874989 C(1)*X^N +. . . + C(N)*X + C(N+1) = 0
符号 具箱 r = solve(‘ 1/x = x-1’); Mapple r= [ 1/2*5^(1/2)+1/2] phi = r(1) [ 1/2 -1/2*5^(1/2)] 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’)
结果图 1/x-(x-1) 7 6 5 4 3 2 1 0 -1 -2 -3 0 0. 5 1 1. 5 2 x 2. 5 3 3. 5 4
%GOLDRECT Golden Rectangle % GOLDRECT plots the golden rectangle phi = (1+sqrt(5))/2; x = [0 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')
连分式 Goldfract. m g= 1+1/(1+1/(1+1/(1)))))) g= 字符串 21/13 g= 1. 61538462 有效数字 2位 err = 0. 0026
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
Fibonacci序列 1 2 3 5 8 13 21 Fibonacci. m 递归实现: Fibnum. m 34 55 89 144 233
递归实现是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 = 40; f = fibonacci(n); f(2: n). /f(1: n-1) - phi
Fibonacci的兔子 注意:没有半只兔子 format long e n = (1: 40)'; f = (phi. ^(n+1)-(1 -phi). ^(n+1))/(2*phi-1) f = round(f)
魔方阵 8 1 6 A = magic(3) 3 5 7 sum(A) 4 9 2 sum(A’) sum(diag(A)) sum(diag(flipud(A))) sum(1: 9)/3
基本代数运算 det(A) X = inv(A) format rat -360 0. 1472 -0. 1444 0. 0639 -0. 0611 0. 0222 0. 1056 -0. 0194 0. 1889 -0. 1028 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) r= 15 e= 15. 0000 4. 8990 -4. 8990 s= 15. 0000 6. 9282 3. 4641 NORM(X) is the largest singular value of X, max(svd(X)).
符号计算 A = sym(A) sum(A’)’ det(A) inv(A) eig(A) svd(A)
图象显示 load durer whos image(X) 100 200 colormap(map) 300 axis image 400 Melancolia, a Renaissance etching by Albrect Durer 500 600 100 200 300 400 500
4阶魔方阵 load detail 50 image(X) colormap(map) axis image 100 150 200 250 300 350 50 100 150 200 250 300 350
4阶魔方阵 A = magic(4) 16 2 3 13 A = A(: , [1 3 2 4]) 5 11 10 8 9 7 6 4 14 15 4阶魔方阵: 880个 5阶魔方阵: 275305224个 6阶魔方阵:? 12 1
秩和奇异阵 for n = 1: 24, r(n) = rank(magic(n)); end [(1: 24)’ r’] bar(r) tiltle(‘魔方阵的秩’)
25 20 15 10 5 0 0 5 10 n: 奇数,满秩 n=4 k: 秩是 3 n=2 k !=4 k: 秩是n/2+2 15 20 25
魔方阵的三维显示 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
模运算 x = [37 – 37]’; y = [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’ ’ 1 U’ If y = Ax mod P then x = Ay mod P Crypto. m mod(A^2, P) = I ‘hello world’ , ? {p]>~Y&=
数论中未解决的3 n+1问题 如果n = 1,停止 如果n是偶数,n = n/2; 如果n是奇数,n = 3 n+1 问题:这个过程是否一定会终止? 综述文章:American Mathematical Monthly, 92(1985), 3 -23 http: //www. cecm. sfu. ca/organics/papers/lagarias
threenplus 1. m
function threenplus 1(n) %"Three n plus 1". % Study the 3 n+1 sequence. % threenplus 1(n) plots the sequence starting with n. % threenplus 1 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'), '3 n+1') shg clf reset uicontrol('position', [260 5 25 22], . . . 'string', '<', . . . 'callback', 'threenplus 1(''<'')'); uicontrol('position', [300 5 25 22], . . . 'string', '>', . . . 'callback', 'threenplus 1(''>'')'); set(gcf, 'tag', '3 n+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(log 2(ymax))-1) ymax]; if length(ytick) > 8, ytick(end-1) = []; end set(gca, 'ytick', ytick) title(['n = ' num 2 str(y(1))]);
x = 32768 y=0 L: 画园 Circlegen. m lcad y shift right 5 bits h = 1/32; add x x = 1; store in x y = 0; change sign while 1 shift right 5 bits x = x + h*y; add y y = y - h*x; store in y plot(x, y, ’. ’) plot x y drawnow go to L end
画园中的特征值 [V, E] = eig(A) E是对角阵 h = 2*rand, A=[1 h; -h 1 -h^2], lambda=eig(A), abs(lambda)
特征值的符号运算 syms h A = [1 h; -h 1 -h^2] d = det(A) lambda = eig(A) d= simple(prod(lambda)) 设 theta = acos(1 -h^2/2); Lambda=[cos(theta)- i*sin(theta); cos(theta)+i*sin(theta)] diff = simple(lambda-Lambda)
精度和舍入误差 t = 0. 1 format hex t 3 fb 999999 a 3 fbh = 1019 试一下:-0. 1是如何存的? 看一下: 0. 3/0. 1和3的存储 0: 0. 1: 1 的问题 format long x = 4/3 – 1 y = 3*x z=1 -y 1019 -1023 = - 4
浮点数 Binary Decimal Hex eps 2^(-52) 2. 2204 e-16 3 cb 0000000 realmin 2^(-1022) 2. 2251 e-308 0010000000 realmax (2 -eps)*2^1023 1. 7977 e+308 7 fefffffff 上溢,下溢 inf 无穷大: e = 1024, f = 0; 7 ff 0000000 nan 不是数:e=1024, f != 0; fff 8000000
舍入误差和多项式 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) 5 x 10 -14 4 3 2 1 0 -1 -2 -3 -4 -5 0. 985 0. 995 1 1. 005 1. 015
ezplot (x-1)7 5 x 10 1 0. 5 0 -0. 5 -1 -1. 5 -2 -2. 5 -3 -3. 5 -4 -6 -4 -2 0 x 2 4 6
- Slides: 59