上海复旦大学matlab课件ppt.ppt

《上海复旦大学matlab课件ppt.ppt》由会员分享，可在线阅读，更多相关《上海复旦大学matlab课件ppt.ppt（59页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Matlab MathCleve Morler著陈文斌（)复旦大学2002年http:/www.stanford.edu/class/cs138/mlmath.html黄金分割 p=1 1 1;r=roots(p);r=-0.61803398874989 1.61803398874989 C(1)*XN+.+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)符号工具箱MappleMapple phi=r(1)vpa(phi,50)1/2*5(1/2)+1/21.6180339887498948482045868343

2、656381177203091798058 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)00.511.522.533.54-3-2-101234567x1/x-(x-1)结果图 ff-111黄金分割图%GOLDRECT Golden Rectangle%GOLDRECT plots the golden rectanglephi=(1+sqrt(5)/2;x=0 phi phi 0 0;y=0 0 1 1 0;u=1 1;v=0 1;

3、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 equalaxis offset(gcf,color,white)连分式连分式Goldfract.mg=1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1)g=21/13g=1.61538461538462err=0.0026字符串有效数字2位Fibonacci Leonardo Pisano Fibonacci was born around 1170 and died aroun

4、d 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 hi

5、s 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 awall.How many pairs of rabbits can be produced from that pair in ayear if it is supposed that every month each pair begets a new pairwhich from the second

6、 month on becomes productive?兔子的故事兔子从出生要一个月，从出生到成熟要过一个月 月123456繁殖111235出生011235成熟001123总数1235813Fibonacci序列Fibonacci.m123581321345589144 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

7、-1)+f(k-2);n=40;f=fibonacci(n);f(2:n)./f(1:n-1)f(2:n)./f(1:n-1)-phiFibonacci的兔子Fibonacci的兔子以黄金分割的速度增长。有两个解phi和(1-phi)由初始条件：是有理分式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)/38 1 63 5 74 9 2魔方

8、阵的八种组合 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.8990s=15.0000 6.9282 3.4641r=15

9、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)100200300400500100200300400500600图象显示 load durer whos image(X)colormap(map)axis imageMelancolia,a Renaissance etching by Albrect Durer50100150200250300350501001502002503003504阶魔方阵 load detail image(

10、X)colormap(map)axis image4阶魔方阵 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(魔方阵的秩)05101520250510152025 n:奇

11、数，满秩 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 tightend图形对象小鸭子加 密Matlab字符能力 x=reshape(32:127,32,3)c=char(x)char(32)c=!#$%&()*+,-./0123456789:;?ABCDEFGHIJKLMNOPQRSTUVWXYZ_abcdefghijk

12、lmnopqrstuvwxyz|字符串 d=char(48:57)double(d)0d=0123456789ans=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=d

13、ouble(TV)-32 y=Ax mod P A=71 2 2 26TV1U If y=Ax mod P then x=Ay mod P mod(A2,P)=I Crypto.mhello world,?pY&=分形蕨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/lagariasthreenplus1.mfunction three

14、nplus1(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,2

15、60 5 25 22,.string,.callback,threenplus1(,.callback,threenplus1();set(gcf,tag,3n+1);endif nargin=0 n=1;elseif isequal(n,)n=get(gcf,userdata)+1;endif n 1 if rem(n,2)=0 n=n/2;else n=3*n+1;end y=y n;endsemilogy(y,.-)axis tightymax=max(y);ytick=2.(0:ceil(log2(ymax)-1)ymax;if length(ytick)8,ytick(end-1)=

16、;endset(gca,ytick,ytick)title(n=num2str(y(1);x=32768 y=0L: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 endCirclegen.m画 园画园中的特征值V,E=eig(A)E是对角阵 h=2*rand,A=1 h;-h 1-h2,lambda=eig

17、(A),abs(lambda)特征值的符号运算 syms h A=1 h;-h 1-h2 lambda=eig(A)d=det(A)d=simple(prod(lambda)设 theta=acos(1-h2/2）;Lambda=cos(theta)-i*sin(theta);cos(theta)+i*sin(theta)diff=simple(lambda-Lambda)浮点数1985年，IEEE标准委员会和美国国家标准局对二进制浮点数运算建立ANSI/IEEE 745-1985标准。Matlab采用IEEE双精度浮点运算。f最多用52位表示规范数表示精度与f有关,表示的数的范围与e有关指数

18、e是整数双精度浮点数双精度浮点数以64位存储，f:52位,e:11位,数的符号：1位。指数以e+1023存储(1,211-2)整个浮点数的分数部分不是f,而是1+f，IEEE标准在64位中存储了65位的信息。Floatgui.meps精度和舍入误差 eps是1到一个最近浮点数的距离,eps=2(-t)对IEEE双精度：eps=2(-52)2.2204*10(-16)最大的相对误差是eps/2，数之间的最大的相对距离是eps t=0.1不能被被精确存储第1项的1不用存储，所以重复出现1001，知道第53项，而54项以后的级数大于eps/2,舍入到53项精度和舍入误差 t=0.1 format h

19、ex t 3fb999999999999a3fbh=10191019-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 Hexeps 2(-52)2.2204e-16 3cb0000000000000realmin 2(-1022)2.2251e-308 0010000000000000realmax (2-eps)*21023 1.7977e+308 7fefffffffffffff上溢，下溢 inf 无穷大:e=1024,f=0;7ff0000

20、000000000 nan 不是数：e=1024,f!=0;fff8000000000000subnormal浮点数许多机器支持subnormal浮点数：realmin eps*realmin标志：e=-1023,例如eps*realmin 0000000000000001 Matlab也用浮点数表示整数，我们用flint来指这些整数，在flint浮点运算不引入舍入误差，加减乘的结果还是flint，当除和平方根的结果是整数时，也产生flint舍入误差和方程求解 A=17 5;1.7 0.5 b=22;2.2 x=Abx=-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)/170.9850.990.99511.0051.011.015-5-4-3-2-1012345x 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-20246-4-3.5-3-2.5-2-1.5-1-0.500.51x 105x(x-1)7ezplot