數(shù)學(xué)軟件Matlab初步

PAGE1數(shù)學(xué)軟件Matlab初步肖水明南昌大學(xué)數(shù)學(xué)系2009年8月一、MATLAB小傳與外觀(guān)圖1.Matlab界面元素二、初探Matlab－從線(xiàn)性代數(shù)談起從線(xiàn)性代數(shù)談起三、Matlab求解簡(jiǎn)單的數(shù)學(xué)模型1.1“黃金比率”問(wèn)題世界上最有趣的數(shù)字是什么？也許你喜歡，或是e，或是17。在這里我們選擇黃金比率。在數(shù)學(xué)的學(xué)習(xí)過(guò)程中，我們會(huì)遇到很多黃金比率方面的應(yīng)用。為什么叫黃金比率？黃金比率到底是什么呢？黃金比率得名于一個(gè)黃金矩形（如圖1.1）。黃金矩形有這樣的性質(zhì)：將圖中的正方形割去后，剩下來(lái)的小矩形和原矩形形狀相同。圖1.1黃金矩形由兩矩形對(duì)應(yīng)邊成比例，容易得到一個(gè)關(guān)于的方程：這個(gè)方程說(shuō)明了黃金比率比的一個(gè)性質(zhì)：的倒數(shù)等于減1。方程兩邊同乘，整理得關(guān)于的一元二次方程：解這個(gè)方程得其中正根即為黃金比率。在Matlab命令窗口中輸入：phi=(1+sqrt(5))/2得到phi=1.6180多顯示幾位數(shù)字：formatlongphi得到phi=1.61803398874989上面的命令沒(méi)有重新計(jì)算，而是顯示15位有效數(shù)字。這里，我們也可以借助MATLAB來(lái)解方程。在MATLAB中，向量又可以表示一個(gè)降次排列的多項(xiàng)式。如p=[1-1-1]表示多項(xiàng)式p(x)=x2+x+1利用roots函數(shù)可計(jì)算多項(xiàng)式的根。r=roots(p)得到r=-0.618033988749891.61803398874989這兩個(gè)數(shù)就是僅有的滿(mǎn)足前面提到那個(gè)性質(zhì)的數(shù)，即減去1等于自身倒數(shù)的數(shù)。還可以使用MATLAB中的符號(hào)工具箱（SymbolicToolbox）來(lái)解這個(gè)比例方程，而不必把它先化為多項(xiàng)式。方程用一字符串表示，再用solve命令求解方程。命令如下：r=solve('1/x=x-1')得到r=[1/2*5^(1/2)+1/2][1/2-1/2*5^(1/2)]r實(shí)際上就是一個(gè)由多項(xiàng)式p的兩個(gè)解所構(gòu)成的二維向量，其第一個(gè)元素為phi=r(1)得到phi=1/2*5^(1/2)+1/2可以用兩種方法將這個(gè)含有根式的結(jié)果化為十進(jìn)制的小數(shù)形式。一種方法是使用vpa（variable-precisionarithmetic）函數(shù)，它可以計(jì)算出任意多有效位數(shù)的值，如：vpa(phi，50)得有50位有效數(shù)字1.6180339887498948482045868343656381177203091798058做實(shí)驗(yàn)請(qǐng)?jiān)囈辉嚸睿簐pa(pi，50).另一種方法是使用double函數(shù)，它表示MATLAB中最常用的數(shù)字類(lèi)型，即雙精度浮點(diǎn)數(shù)。Phi=double(phi)得到phi=1.61803398874989在MATLAB中可將函數(shù)表達(dá)式用inline函數(shù)來(lái)表示。命令：f=inline('1/x-(x-1)');定義了函數(shù)f(x)=1/x-(x-1)，運(yùn)行后得到：f=Inlinefunction:f(x)=1/x-(x-1)利用下面函數(shù)，我們可以作出函數(shù)f(x)在區(qū)間[0，4]上的圖像。ezplot(f，0，4).函數(shù)名ezplot表示“easyplot”，讀為“e-zedplot”。表達(dá)式phi=fzero(f，1)表示數(shù)值求解f(x)在x=1附近的零點(diǎn)。以下命令表示將這樣得到的的近似解，與ezplot畫(huà)出的圖像在同一張圖中顯示，結(jié)果見(jiàn)圖1－2。holdonplot(phi，0，'o')注意，以上命令中，第一個(gè)0是數(shù)字零，單引號(hào)中的是字母o。圖1.2f()=0下面我們來(lái)用MATLAB作出圖1.1中的黃金矩形。程序?yàn)椋?GOLDRECTPlotthegoldenrectanglephi=(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')以上程序保存在一個(gè)名為goldrect.m的M-file中，在Matlab命令窗口中使用命令goldrect就能運(yùn)行程序并創(chuàng)建圖像。下面，我們來(lái)詳細(xì)解釋該程序。第一行為注釋行。x和y均為含有五個(gè)分量的向量，分別表示圖1.1中長(zhǎng)方形四個(gè)頂點(diǎn)的橫縱坐標(biāo)；向量u和v含有兩個(gè)分量，分別表示圖1.1中虛線(xiàn)兩端點(diǎn)與的橫縱坐標(biāo)。語(yǔ)句plot(x,y,'b',u,v,'b--')表示依次用藍(lán)色直線(xiàn)連接點(diǎn)和藍(lán)色的虛線(xiàn)連接點(diǎn)與，其中的字母b表示藍(lán)色（blue）.接下來(lái)的四條text語(yǔ)句，是表示在各個(gè)指定的位置用符號(hào)作標(biāo)記。字符串'\phi'表示希臘字母。兩條axis語(yǔ)句分別表示將x軸和y軸方向比例調(diào)整為一致，然后再關(guān)閉顯示坐標(biāo)軸。最后一條語(yǔ)句是將圖像的背景顏色設(shè)置為白色，其中的gcf意為getcurrentfigure.下面，我們從連分?jǐn)?shù)的角度來(lái)研究黃金比率。連分?jǐn)?shù)是一種可以無(wú)限延伸擴(kuò)展的表達(dá)式將上式中所有的都設(shè)置為1，那么所得到的連分?jǐn)?shù)就可以視為黃金比率的另一種表示：下面的函數(shù)就可以得到截取連分?jǐn)?shù)上面部分，得到近似于的值。這些語(yǔ)句儲(chǔ)存在名為goldfract.m的M-file中。functiongoldfract(n)%GOLDFRACTGoldenratiocontinuedfraction.%GOLDFRACT(n)displaysnterms.p='1';fork=1:np=['1+1/('p')'];endpp=1;q=1;fork=1:ns=p;p=p+q;q=s;endp=sprintf('%d/%d',p,q)formatlongp=eval(p)formatshorterr=(1+sqrt(5))/2-p得到p=1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))p=21/13p=1.61538461538462err=0.0026這三個(gè)p為近似值的三種不同形式。第一個(gè)p是由連分?jǐn)?shù)截取上面六的層（即包含六個(gè)右括?。┒玫?。這個(gè)p僅用1（1=goldfract(0)）開(kāi)始，并重復(fù)地在前面插入'1+1/（'，而同時(shí)在后面插入'）'而得到的。不管這一串式子多長(zhǎng)，它都是一個(gè)有效的MATLAB表達(dá)式。第二個(gè)p是將第一個(gè)p化簡(jiǎn)為一個(gè)普通分?jǐn)?shù)。這個(gè)變形的基本思路就是：因此這個(gè)重復(fù)變形過(guò)程始于并重復(fù)用新的數(shù)代替原來(lái)的最后再用語(yǔ)句sprint('%d/%d'，p，q)就能打印出最后的結(jié)果，其中p和q分別為十進(jìn)制的整數(shù)，它們之間用'/'分隔。第三個(gè)p實(shí)際上和前面兩個(gè)p是同一個(gè)數(shù)，只不過(guò)是利用第二個(gè)p作了除法而得到的一個(gè)十進(jìn)制的小數(shù)表達(dá)形式。1.2Fibonacci數(shù)列LeonardoPisanoFibonacci（公元1170—1250）生于Pisa即現(xiàn)在的意大利。他幾乎游遍歐洲和北非。除了寫(xiě)了為數(shù)不多的數(shù)學(xué)論文，他還向印度和阿拉伯人介紹歐洲。盡管他的書(shū)是手寫(xiě)的，但是廣為傳閱。其中最有名的一本名為L(zhǎng)iberAbaci出版于1202年，在書(shū)中他提出了一個(gè)有名的問(wèn)題：一對(duì)兔子，假使它們每個(gè)月生一對(duì)兔子，而每對(duì)新出生的兔子要在第二個(gè)月才有生育能力。新出生的兔子從第二月開(kāi)始每個(gè)月生一對(duì)兔子，那么一年之后總共有多少對(duì)兔子？Amanputapairofrabbitsinaplacesurroundedonallsidesbyawall.Howmanypairsofrabbitscanbeproducedfromthatpairinayearifitissupposedthateverymontheachpairbegetsanewpairwhichfromthesecondmonthonbecomesproductive?這個(gè)問(wèn)題被稱(chēng)為Fibonacci序列或者Fibonacci數(shù)列問(wèn)題。在因特網(wǎng)上搜索“Fibonacci”，我們可以找到很多介紹Fibonacci數(shù)列的網(wǎng)站。甚至有Fibonacci協(xié)會(huì)，這個(gè)協(xié)會(huì)還出版了一本專(zhuān)業(yè)性很強(qiáng)的雜志《FibonacciQuarterly》。如果Fibonacci沒(méi)有這樣一個(gè)假設(shè)：一對(duì)剛出生的兔子從第二月開(kāi)始每個(gè)月生一對(duì)兔子，也許這個(gè)數(shù)列就不會(huì)以他的名字命名。因?yàn)?，此時(shí)問(wèn)題變得很簡(jiǎn)單：一個(gè)月之內(nèi)的兔子總數(shù)是上個(gè)月的兩倍，n個(gè)月之后的兔子總數(shù)為對(duì)。仔細(xì)分析可以發(fā)現(xiàn)，每個(gè)月的兔子總數(shù)剛好是這個(gè)月以前的兔子的總數(shù)加上在這個(gè)月由成熟兔子所生的兔子數(shù)，若用表示n個(gè)月之后的兔子總數(shù)，則有：.這個(gè)問(wèn)題的初始條件是：第一個(gè)月有一對(duì)兔子，第二個(gè)月有兩對(duì)兔子即：,.我們可以用Matlab的M文件fibonacci.m來(lái)產(chǎn)生一個(gè)由n個(gè)Fibonacci數(shù)所組成的向量。程序是：functionf=fibonacci(n)%FIBONACCIFibonaccisequence%f=FIBONACCI(n)generatesthefirstnFibonaccinumbers.f=zeros(n,1)；f(1)=1；f(2)=2；fork=3:nf(k)=f(k-1)+f(k-2);end一年后兔子總數(shù)可以用fibonacci(12)來(lái)計(jì)算。結(jié)果為：123581321345589144233結(jié)果是一年后有233對(duì)兔子（如果每月成倍增長(zhǎng)，一年后將有4096對(duì)兔子）。讓我們來(lái)解釋這個(gè)

M文件fibonacci.m中的每一條語(yǔ)句。第一行functionf=fibonacci(n)中function是一個(gè)關(guān)鍵詞，它表明該M文件是一個(gè)函數(shù)，而f=fibonacci(n)是指明這個(gè)函數(shù)的名字及輸入?yún)?shù)n，并把結(jié)果賦予變量f。注意，此處的函數(shù)名應(yīng)與保存該M文件的文件名一致。接下來(lái)的這兩句是起注釋用的，其目的是方便使用者通過(guò)help命令來(lái)了解本函數(shù)的功能和用法。輸入：helpFibonacci結(jié)果：FIBONACCIFibonaccisequencef=FIBONACCI(n)generatesthefirstnFibonaccinumbers.FIBONACCI這個(gè)函數(shù)名字是大寫(xiě)，這是因?yàn)槔习姹镜腗atlab不區(qū)分大小寫(xiě)。第一次用Matlab的用戶(hù)會(huì)對(duì)這種大寫(xiě)字母不習(xí)慣，但這是歷史原因造成的。下一句：f=zeros(n,1);產(chǎn)生一個(gè)n1矩陣，初始值為n1的零矩陣。在Matlab中一個(gè)矩陣只有一行（列），它就是一個(gè)行（列）行量。緊接著的兩句：f(1)=1;f(2)=2;為初始條件。最后三句才是最重要的：fork=3:nf(k)=f(k-1)+f(k-2);end在使用for或者if語(yǔ)句時(shí)，我們通常采用與其他高級(jí)語(yǔ)言類(lèi)似的縮進(jìn)格式書(shū)寫(xiě)，以方便閱讀。我們可以用另一種方式構(gòu)造Fibonacci函數(shù)，它的結(jié)果僅為第n個(gè)Fibonacci數(shù)。functionf=fibnum(n)%FIBNUMFibonaccinumber.%FIBNUM(n)generatesthenthFibonaccinumber.ifn<=1f=1;elsef=fibnum(n-1)+fibnum(n-2);end輸入語(yǔ)句：fibnum(12)結(jié)果為：ans=233這個(gè)fibnum函數(shù)采用了遞歸調(diào)用方式。這個(gè)關(guān)系就是一個(gè)遞歸關(guān)系，是一個(gè)自己調(diào)用自己的遞歸關(guān)系。遞歸非常好用，但是也要付出代價(jià)。你可用Matlab中的tic和toc命令來(lái)估計(jì)執(zhí)行命令所需的時(shí)間。試一下tic,fibnum(24),toc注意，不要嘗試tic,fibnum(50),toc下面，我們來(lái)研究黃金比率和Fibonacci數(shù)列的關(guān)系。首先，輸入goldfract(6)和fibonacci(7).第一個(gè)命令的結(jié)果21/13，第二個(gè)命令的最后兩個(gè)數(shù)是13和21！這不是一個(gè)巧合。連分?jǐn)?shù)部分是由前面相加而得，即不斷執(zhí)行p=p+q;語(yǔ)句。而Fibonacci數(shù)列是由f(k)=f(k-1)+f(k-2);產(chǎn)生的。事實(shí)上，如果我們用來(lái)表示黃金比率點(diǎn)序列的第n個(gè)點(diǎn)，則=。當(dāng)n趨于無(wú)窮，將無(wú)窮小地趨于黃金比率點(diǎn)即：=。為了看清這一事實(shí)，先計(jì)算40個(gè)Fibonacci數(shù)，n=40;f=fibonacci(n);然后計(jì)算比率：f(2:n)./f(1:n-1);這條語(yǔ)句產(chǎn)生了一個(gè)向量（，，…，），前9個(gè)分量如下；2.000000000000001.500000000000001.666666666666671.600000000000001.625000000000001.615384615384621.619047619047621.617647058823531.61818181818182最后5個(gè)分量為：1.618033988749901.618033988749891.618033988749901.618033988749891.61803398874989你有沒(méi)有看出為什么我們選n=40呢？我們把f(2:n)./f(1:n-1);改為：f(2:n)./f(1:n-1)-phi;敲一下回車(chē)鍵，最后一個(gè)分量會(huì)是多少？兔子的總數(shù)并不是成倍增長(zhǎng)，而是每個(gè)月是以接近黃金比率增長(zhǎng)，我們用一個(gè)函數(shù)表達(dá)式來(lái)表示Fibonacci數(shù)，設(shè) =c這里c和是常數(shù)，根據(jù)=+可得：=+1這是一個(gè)二元一次方程，這個(gè)方程有兩個(gè)根和，于是有=+這里的和是由初始條件所決定的常數(shù)，代入初始條件:=1，=1可得：=+=1=+（1-）=1我們可以用Matlab的“\”(backslashoperator)來(lái)解這個(gè)二元一次線(xiàn)性方程組。我們也可以筆算這個(gè)方程組，可得：==。把和代入中可得：=（—）。這是一個(gè)令人困惑的方程，右邊既有乘方又有除法，似乎不是一個(gè)整數(shù)。用Matlab來(lái)驗(yàn)證，