數(shù)學軟件Matlab初步

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

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