上海复旦大学matlab课件

MatlabMathCleveMorler著陈文斌（wbchen@)复旦大学2002年MatlabMathCleveMorler著1/class/cs138/mlmath.html/class/2黄金分割黄金分割3p=[1–1–1];r=roots(p);r=-0.618033988749891.61803398874989C(1)*X^N+...+C(N)*X+C(N+1)=0p=[1–1–1];r=C(1)*X^N+4r=solve(‘1/x=x-1’);r=[1/2*5^(1/2)+1/2][1/2-1/2*5^(1/2)]符号工具箱Mapplephi=r(1)vpa(phi,50)1/2*5^(1/2)+1/21.6180339887498948482045868343656381177203091798058phi=double(phi)1.61803398874989r=solve(‘1/x=x-1’);r=符号5图形的方法f=inline('1/x-(x-1)')ezplot(f,0,4)phi=fzero(f,1)holdonplot(phi,0,’o’)图形的方法f=inline('1/x-(x-1)')600.511.522.533.54-3-2-101234567x1/x-(x-1)结果图00.511.522.533.54-3-2-101234567ff-111黄金分割图ff-111黄金分割图8%GOLDRECTGoldenRectangle%GOLDRECTplotsthegoldenrectanglephi=(1+sqrt(5))/2;x=[0phiphi00];y=[00110];u=[11];v=[01];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')axisequalaxisoffset(gcf,'color','white')%GOLDRECTGoldenRectangle9连分式连分式10连分式Goldfract.mg=1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))g=21/13g=1.61538461538462err=0.0026字符串有效数字2位连分式Goldfract.mg=字符串有效数字2位11Fibonacci

LeonardoPisanoFibonacciwasbornaround1170anddiedaround1250inPisainwhatisnowItaly.HetraveledextensivelyinEuropeandNorthernAfrica.Hewroteseveralmathematicaltextsthat,amongotherthings,introducedEuropetotheHindu-Arabicnotationfornumbers.Eventhoughhisbookshadtobetranscribedbyhand,theywerestillwidelycirculated.Inhisbestknownbook,LiberAbaci,publishedin1202,heposedthefollowingproblem.Amanputapairofrabbitsinaplacesurroundedonallsidesbyawall.Howmanypairsofrabbitscanbeproducedfromthatpairinayearifitissupposedthateverymontheachpairbegetsanewpairwhichfromthesecondmonthonbecomesproductive?FibonacciLeonardoPisanoFibo12兔子的故事兔子从出生要一个月，从出生到成熟要过一个月月123456繁殖111235出生011235成熟001123总数1235813兔子的故事兔子从出生要一个月，从出生到成熟要过一个月13Fibonacci序列Fibonacci.m123581321345589144233递归实现：Fibnum.mFibonacci序列Fibonacci.m1递归实现：Fi14递归实现递归实现是elegantbutexpensivetic,fibnum(24),toctic,fibonacci(24),toc试一下，比较时间！递归实现递归实现是elegantbutexpensive15黄金分割和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)f(2:n)./f(1:n-1)-phi黄金分割和Fibonacci数比较一下goldfract(616Fibonacci的兔子Fibonacci的兔子以黄金分割的速度增长。＋有两个解phi和(1-phi)由初始条件：是有理分式Fibonacci的兔子Fibonacci的兔子以黄金分割的17Fibonacci的兔子注意：没有半只兔子

formatlongen=(1:40)';f=(phi.^(n+1)-(1-phi).^(n+1))/(2*phi-1)

f=round(f)Fibonacci的兔子注意：没有半只兔子formatl18魔方阵A=magic(3)sum(A)sum(A’)sum(diag(A))sum(diag(flipud(A)))sum(1:9)/3816357492魔方阵A=magic(3)sum(1:9)/3819魔方阵的八种组合fork=0:3rot90(A,k)rot90(A’,k)end逆时阵旋转k*90度魔方阵的八种组合fork=0:3逆时阵旋转k*90度20基本代数运算X=inv(A)det(A)-3600.1472-0.14440.0639-0.06110.02220.1056-0.01940.1889-0.1028formatrat53/360-13/9023/360-11/1801/4519/180-7/36017/90-37/360基本代数运算X=inv(A)det(A)-360021基本代数运算r=norm(A)e=eig(A)s=svd(A)e=15.00004.8990-4.8990s=15.00006.92823.4641r=15NORM(X)isthelargestsingularvalueofX,max(svd(X)).基本代数运算r=norm(A)e=s=r=NO22符号计算A=sym(A)sum(A)sum(A’)’det(A)inv(A)eig(A)svd(A)符号计算A=sym(A)23100200300400500100200300400500600图象显示loaddurerwhosimage(X)colormap(map)axisimageMelancolia,aRenaissanceetchingbyAlbrectDurer1002003004005001002003004005002450100150200250300350501001502002503003504阶魔方阵loaddetailimage(X)colormap(map)axisimage501001502002503003505010015020254阶魔方阵16231351110897612414151A=magic(4)A=A(:,[1324])4阶魔方阵：880个5阶魔方阵：275305224个6阶魔方阵：？4阶魔方阵16226秩和奇异阵A=magic(4)det(A)inv(A)矩阵的秩是矩阵中线性无关的行和列的个数，一个n阶矩阵是奇异的当且仅当它的秩小于n

rank(magic(4))秩和奇异阵A=magic(4)矩阵的秩是矩阵中线性27秩和奇异阵forn=1:24,r(n)=rank(magic(n));end[(1:24)’r’]bar(r)tiltle(‘魔方阵的秩’)秩和奇异阵forn=1:24,r(n)=2805101520250510152025n:奇数，满秩n=4k:秩是3

n=2k!=4k:秩是n/2+205101520250510152025n:奇数，满秩29魔方阵的三维显示surf(magic(n))axisoffset(gcf,’doublebuffer’,’on’)cameratoolbarforn=8:11

subplot(2,2,n-7),surf(magic(n))axisoff,view(30,45)%axistightend魔方阵的三维显示surf(magic(n))forn30图形对象图形对象31上海复旦大学matlabppt课件32小鸭子小鸭子33加密－－Matlab字符能力x=reshape(32:127,32,3)’c=char(x)char(32)c=!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~加密－－Matlab字符能力x=reshape(334字符串d=char(48:57)double(d)–‘0’d=0123456789ans=0123456789字符串d=char(48:57)d=35字符串（高八位）char(reshape(160:255,32,3)’)char(160)char(160:161)牎ⅲぅΗī氨渤吹斗腹夯冀究懒旅呐魄壬仕掏蜗醒矣哉肿刭谯茌捱噌忏溴骁栝觌祉铒瘃蝮趱鲼?字符串（高八位）char(reshape(160:255,36模运算x=[37–3737–37]’;y=[1010-10-10]’;r=[xyrem(x,y)mod(x,y)]371077-3710-7337-107-3-37-10-7-7模运算x=[37–3737–37]’;337字符加密取素数：P=97(要表示95个字符）x=double(‘TV’)’-32y=AxmodPA=[712226]‘TV’’1U’

Ify=AxmodPthenx=AymodPmod(A^2,P)=I

Crypto.m‘helloworld’

,?{p]>~Y&=字符加密取素数：P=97(要表示95个字符）x=38分形蕨Fern.m分形蕨Fern.m39数论中未解决的3n+1问题如果n=1，停止如果n是偶数，n=n/2;如果n是奇数，n=3n+1问题：这个过程是否一定会终止？综述文章：AmericanMathematicalMonthly,92(1985),3-23

http://www.cecm.sfu.ca/organics/papers/lagarias数论中未解决的3n+1问题如果n=1，停止问题：这个过程40threenplus1.mthreenplus1.m41functionthreenplus1(n)%"Threenplus1".%Studythe3n+1sequence.%threenplus1(n)plotsthesequencestartingwithn.%threenplus1withnoargumentsstartswithn=1.%uicontrolsdecrementorincrementthestartingn.%Isitpossibleforthistorunforever?functionthreenplus1(n)42if~isequal(get(gcf,'tag'),'3n+1')shgclfresetuicontrol('position',[26052522],...'string','<',...'callback','threenplus1(''<'')');uicontrol('position',[30052522],...'string','>',...'callback','threenplus1(''>'')');set(gcf,'tag','3n+1');endif~isequal(get(gcf,'tag'),'3n43ifnargin==0n=1;elseifisequal(n,'<')n=get(gcf,'userdata')-1;elseifisequal(n,'>')n=get(gcf,'userdata')+1;endifn<1,n=1;endset(gcf,'userdata',n)ifnargin==044y=n;whilen>1ifrem(n,2)==0n=n/2;elsen=3*n+1;endy=[yn];endy=n;45semilogy(y,'.-')axistightymax=max(y);ytick=[2.^(0:ceil(log2(ymax))-1)ymax];iflength(ytick)>8,ytick(end-1)=[];endset(gca,'ytick',ytick)title(['n='num2str(y(1))]);semilogy(y,'.-')46x=32768

y=0L:lcadyshiftright5bitsaddxstoreinxchangesignshiftright5bitsaddystoreinyplotxygotoL画园h=1/32;x=1;y=0;while1x=x+h*y;y=y-h*x;plot(x,y,’.’)drawnowendCirclegen.mx=32768画园h=1/32;47画园画园48画园中的特征值[V,E]=eig(A)E是对角阵h=2*rand,A=[1h;-h1-h^2],lambda=eig(A),abs(lambda)画园中的特征值[V,E]=eig(A)E是对角阵h=49特征值的符号运算symshA=[1h;-h1-h^2]lambda=eig(A)d=det(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)特征值的符号运算symshd=det(A)设th50浮点数1985年，IEEE标准委员会和美国国家标准局对二进制浮点数运算建立ANSI/IEEE745-1985标准。

Matlab采用IEEE双精度浮点运算。f最多用52位表示规范数表示精度与f有关,表示的数的范围与e有关指数e是整数浮点数1985年，IEEE标准委员会和美国国家标准局对二进制51双精度浮点数双精度浮点数以64位存储，f:52位,e:11位,数的符号：1位。指数以e+1023存储(1,211-2)整个浮点数的分数部分不是f,而是1+f，IEEE标准在64位中存储了65位的信息。Floatgui.meps双精度浮点数双精度浮点数以64位存储，f:52位,e:11位52精度和舍入误差eps是1到一个最近浮点数的距离,eps=2^(-t)对IEEE双精度：eps=2^(-52)≈2.2204*10^(-16)最大的相对误差是eps/2，数之间的最大的相对距离是epst=0.1不能被被精确存储第1项的1不用存储，所以重复出现1001，知道第53项，而54项以后的级数大于eps/2,舍入到53项精度和舍入误差eps是1到一个最近浮点数的距离,eps=253精度和舍入误差t=0.1formathext3fb999999999999a3fbh=10191019-1023=-4试一下：-0.1是如何存的？看一下：0.3/0.1和3的存储0:0.1:1的问题formatlongx=4/3–1y=3*xz=1-y精度和舍入误差t=0.13fb9999999999954浮点数

BinaryDecimalHexeps2^(-52)2.2204e-163cb0000000000000realmin2^(-1022)2.2251e-308