MATLAB應用實例分析例分析

MATLAB應用實例分析例分析Matlab應用例題選講僅舉一些運用MATLAB的例子，這些問題在數(shù)學建模中時常遇到，希望能幫助同學們在短時間內(nèi)方便、快捷的使用MATLAB解決數(shù)學建模中的問題，并善用這一工具。常用控制命令:clc:%清屏；clear:%清變量；save:%保存變量；load:%導入變量一、利用公式直接進行賦值計算本金P以每年n次，每次i%的增值率(n與i的乘積為每年增值額的百分比)增加，當增加到rXP時所花費的時間T為：(利用復利計息公式可得到下式)lnrnT()r，P,P(1，0.01i),T,r,2,i,0.5,n,12nln(1，0.01i)MATLAB的表達形式及結果如下:>>r=2；i=0.5；n=12；%變量賦值>>T=log(r)/(n*log(1+0.01*i))計算結果顯示為:T=11.5813即所花費的時間為T=11.5813年。分析:上面的問題是一個利用公式直接進行賦值計算問題，實際中若變量在某個范圍變化取很多值時，使用MATLAB，將倍感方便，輕松得到結果，其繪圖功能還能將結果輕松的顯示出來，變量之間的變化規(guī)律將一目了然。若r在［1,9］變化，i在［0.5,3.5］變化;我們將MATLAB的表達式作如下改動，結果如圖1。 r=1:0.5:9；i=0.5:0.5:3.5;n=12;p=1./(n*log(1+0.01*i));T=log(r')*p;plot(r,T)xlabel('r')%給x軸加標題ylabel('T')%給y軸加標題q=ones(1,length(i));text(7*q-0.2,[T(14,1:5)+0.5,T(14,6)-0.1,T(14,7)-0.9],num2str(i'))40350.5302520T1151.51022.5533.50123456789r圖11從圖1中既可以看到T隨r的變化規(guī)律，而且還能看到i的不同取值對T—r曲線的影響(圖中的六條曲線分別代表i的不同取值)。二、已知多項式求根65432已知多項式為，求其根。h,x,10x，31x,10x,116x，200x,96分析:對多項式求根問題，我們常用roots()函數(shù)。MATLAB的表達形式及結果如下:>>h=roots([1-1031-10-116200-96])%中括號內(nèi)為多項式系數(shù)由高階到常數(shù)。計算結果顯示為(其中i為虛數(shù)單位):-2.00004.00003.00002.0000+0.0000i2.0000-0.0000i1.0000如果已知多項式的根，求多項式，用poly()函數(shù)。對上面得到的h的值求多項式，其MATLAB的表達形式及結果如下：>>h=[-2.00004.00003.00002.0000+0.0000i2.0000-0.0000i1.0000];>>c=poly(h)計算結果顯示為：c=1-1031-10-116200-96三、方程組的求解8x，x，6x,7.5,123,求解下面的方程組：3x，5x，7x,4,123,4x，9x，2x,12123,分析：對于線性方程組求解，常用線性代數(shù)的方法，把方程組轉化為矩陣進行計算。,1,x,abax,b,x,a\bMATLAB的表達形式及結果如下：>>a=[816;357;492];%建立系數(shù)矩陣>>b=[7.5;4;12];%建立常數(shù)項矩陣>>x=a\b%求方程組的解計算結果顯示為：1.29310.8972-0.6236四、數(shù)據(jù)擬合與二維繪圖在數(shù)學建模競賽中，我們常會遇到這種數(shù)據(jù)表格問題，如果我們僅憑眼睛觀察，很難看到其中的規(guī)律，也就更難寫出有效的數(shù)學表達式從而建立數(shù)學模型。因此可以利用MATLAB的擬合函數(shù)，即polyfit()函數(shù)，并結合MATLAB的繪圖功能(利用plot()函數(shù)),得到直觀的表示。2例:在化學反應中，為研究某化合物的濃度隨時間的變化規(guī)律，測得一組數(shù)據(jù)如下表:12345678T(分)y46.48.08.49.289.59.79.86910111213141516T(分)y1010.210.3210.4210.510.5510.5810.6分析:MATLAB的表達形式如下:t=[1:16];%數(shù)據(jù)輸入y=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];plot(t,y,'o')%畫散點圖p=polyfit(t,y,2)%二次多項式擬合holdonxi=linspace(0,16,160);%在[0,16]等間距取160個點yi二polyval(p,xi);%由擬合得到的多項式及xi,確定yiplot(xi,yi)%畫擬合曲線圖執(zhí)行程序得到圖2;11109876540246810121416圖2顯示的結果為p=0.04451.07114.3252-2p的值表示二階擬合得到的多項式為:y=-0.04451+1.07111+4.3252下面是用lsqcurvefit()函數(shù)，即最小二乘擬合方法的Matlab表達：t=[1:16];y=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];x0=[0.1,0.1,0.1];zuixiao=inline('x(l)*t「2+x(2)*t+x(3)','x'，‘t');x=lsqcurvefit(zuixiao,x0,t,y)%利用最小二乘擬合其顯示的結果為:x=-0.04451.07114.3252可以看出其得到的結果與polyfit函數(shù)的結果相同。這說明在多項式擬合問題上這兩個函數(shù)的效果是相同的。3下面的一個例子將體現(xiàn)lsqcurvefit()函數(shù)的優(yōu)勢。例2:在物理學中，為研究某種材料應力與應變的關系，測得一組數(shù)據(jù)如下表:925112516252125262531253625應力a0.110.160.350.480.610.710.85應變￡如果假定應力與應變有如下關系(a為應力值，￡為應變值)：￡二a+blna試計算a、b的值。MATLAB的表達形式如下：x=[925,1125,1625,2125,2625,3125,3625];y=[0.11,0.16,0.35,0.48,0.61,0.71,0.85];plot(x,y,'o')[p,resid1]=polyfit(x,y,2)holdonxi=linspace(700,3700,3000);yi=polyval(p,xi);plot(xi,yi)x0=[0.1,0.1];fff=inline('a(1)+a(2)*log(x)','a','x');[a,resid2]=lsqcurvefit(fff,x0,x,y)plot(xi,fff(a,xi),'r')執(zhí)行程序得到圖3,圖中藍色曲線為利用polyfit()函數(shù)得到的曲線，紅色曲線為利用lsqcurvefit()函數(shù)得到的曲線;0.90.80.70.60.50.40.30.20.10-0.15001000150020002500300035004000其顯示的結果為:p=-0.00000.0004-0.2266resid1=R:[3x3double]df:4normr:0.0331a=-3.58100.5344resid2=0.0064其中a的值代表利用lsqcurvefit()函數(shù)得到的關系為:￡=-3.5810+0.5344。4residl、resid2分別代表運用polyfit()函數(shù)、lsqcurvefit()函數(shù)得到的殘差?？梢钥闯隼胠sqcurvefit()函數(shù)殘差更小，即得到了更好的擬合效果。在數(shù)學建模的實際問題中，如果問題的機理不明，我們只能采用polyfit()函數(shù)，即多項式擬合的方法，以獲得近似的數(shù)據(jù)描述函數(shù);但如果通過分析，可以得到一些機理，那么采用最小二乘的方法將得到更好的效果，而且得到的擬合函數(shù)也更有意義。五、隱函數(shù)的圖形繪制1plot()只能繪制顯函數(shù)圖形，對于形如的復雜隱函數(shù)，很難轉化為,lny.ln(,1，y)，x,sin(x),0y顯函數(shù)并利用plot()函數(shù)繪制圖形，這時就可以用ezplot()函數(shù)直接繪制其曲線。MATLAB的表達形式如下:>>ezplot('1/y-log(y)+log(-1+y)+x-sin(x)')執(zhí)行程序得到圖5x=sin(3t)cos(t),y=sin(3t)sin(t)1/y-log(y)+log(-1+y)+x-sin(x)=060.440.220-0.2yy0-0.4-2-0.6-4-0.8-6-0.8-0.6-0.4-0.200.20.40.60.8-6-4-20246xx圖5圖6如果是形如下面的參數(shù)方程，同樣可以利用ezplot()函數(shù)繪制其曲x,sin3tcost,y,sin3tsint,t,(0,,)線。MATLAB的表達形式如下：>>ezplot('sin(3*t)*cos(t)'，'sin(3*t)*sin(t)'，[0，pi])執(zhí)行程序得到圖6。六、三維圖形繪制假設有一個時間向量t,對該向量進行下列運算則可以構成三個坐標值向量x,sint,y,cost,z,t對于上面的方程可以利用ezplot3()函數(shù)或plot3()函數(shù)繪制三維曲線。這里僅列舉ezplot3()函數(shù)的使用。MATLAB的表達形式如下：>>ezplot3('sin(t)','cos(t)','t',[0,6*pi])執(zhí)行程序得到圖7：3繪制下述曲面：z(r,,),rcos(3,),其中0,r,1,0,,,2,MATLAB的表達形式如下:nr=12;nth=50;r=linspace(0,1,nr);theta=linspace(0,2*pi,nth);[R,T]=meshgrid(r,theta)x=cos(theta')*r;5y=sin(theta')*r;surf(x,y,R「3.*cos(3*T))執(zhí)行程序得到圖8。x=sin(t),y=cos(t),z=t201510z5010.510.500-0.5-0.5-1-1yx圖7圖8除了surf()函數(shù)還有surfc()、surfl()、mesh()、waterfall()函數(shù)也用于曲面的繪制，具體效果如圖9所示，可以針對自己的需要選取適合的曲面繪制函數(shù)。圖9MATLAB作圖之三維繪圖示例山體繪制%三維繪圖示例山體繪制%mesh函數(shù)演示x=1.0:0.1:2.0;y=2.0:0.1:3.0;[X,Y]=meshgrid(x,y);z=[5.115.135.145.135.095.044.984.934.894.854.85TOC\o"1-5"\h\z5.395.49 5.51 5.46 5.32 5.14 4.944.744.594.49 4.485.615.77 5.81 5.71 5.51 5.23 4.904.594.364.21 4.195.735.92 5.97 5.86 5.62 5.27 4.884.514.234.05 4.035.745.92 5.97 5.86 7.62 5.27 4.884.514.214.04 4.025.635.79 5.84 6.74 10.539.238.914.594.334.184.165.425.53 5.56 5.49 7.35 5.16 4.934.734.554.45 4.445.145.18 5.19 5.17 11.125.054.974.904.844.814.804.484.80 4.79 4.82 4.87 4.94 5.025.105.165.19 5.204.564.45 4.43 4.49 4.64 4.84 5.065.285.455.55 5.564.364.19 4.16 4.25 4.47 4.76 5.095.415.665.81 5.83];figure(1)%畫網(wǎng)格圖mesh（X,Y,z）;6colormap（［010］）;xlabel（'x軸'）;ylabel（'y軸'）;zlabel（'z軸'）;%畫表面圖figure（2）surf（X,Y,z）;colormap（［100］）;xlabel（'x軸'）;ylabel（'y軸'）;zlabel（'z軸'）;1212101088z軸z軸664433221.81.82.52.51.61.61.41.41.21.22211y軸y軸x軸x軸七、二項分布的使用飛機成功起飛的概率問題:由16架飛機組成的空軍飛行中隊要求做好立即起飛的準備，其中一架飛機不能立即起飛的概率為20%，重新起飛需幾分鐘的時間，因此一架飛機立刻起飛的概率為0.80。12架飛機能夠成功起飛的概率為多少,分析:這是一個概率中的二項分布問題，常用binopdf()函數(shù)。h=binopdf(12,16,0.80)%二項分布函數(shù)的概率值計算結果為h=0.2001另一方面，至少有14架飛機立刻成功起飛的概率為:h=1-binocdf(13,16,0.80)%或h=sum(binopdf(14:16,16,0.80))，其中binocdf()為二項分布的累積概率值，計算結果為h=0.3518。在實際的數(shù)學建模競賽中，僅羅列一個一個的數(shù)據(jù)是枯燥而又不直觀的，很難吸引人們的注意，也不容易打動評委們的心;因此，結合數(shù)值計算結果，并合理的利用MATLAB的繪圖功能會起到事半功倍的效果。下面的程序為運行結果的繪圖(圖10)表示:n=1:16;h=binopdf(n,16,0.80);plot([n;n],[zeros(1,16);h],'k')%二維繪圖函數(shù)text(8-.7:16-.7,h(8:16)+.005,num2str(h(8:16)',3))%在圖中進行注釋函數(shù)axis([01700.27])%坐標軸取值范圍函數(shù)xlabel('Numberofaircraftlaunchedontime')%給x軸加標題ylabel('probability')%給y軸加標題set(gca,'XTick',0:2:16)7set(gca,'XTickLabel',{'0架','2架','4架','6架','8架','10架','12架','14架','16架'})130.2460.2512110.2110.20.21090.1580.120.113probability70.160.05550.050.028140.01970.005533002468101214160架2架4架6架8架10架12架14架