MATLAB數(shù)學(xué)實(shí)驗(yàn)整套教學(xué)課件

1、MATLAB數(shù)學(xué)實(shí)驗(yàn)2022/7/142第一章 Matlab入門第一章 MATLAB入門 1.1 MATLAB桌面 1.2 數(shù)據(jù)和變量1.3 數(shù)組及其運(yùn)算 1.4 字符串、元胞和結(jié)構(gòu) 2022/7/143第一章 Matlab入門1.1 MATLAB桌面最小安裝：MATLAB 2013a，Curve Fitting Toolbox ， Optimization Toolbox ，Symbolic Math Toolbox工具條：Home, Plots, Apps等窗口: 指令(Command Window) 當(dāng)前文件夾(Current Directory) 工作空間(Workspace) 指令歷

2、史(Command History)例: a=1;b=2;c=a+b (不輸入提示符)2022/7/144第一章 Matlab入門1.2 數(shù)據(jù)和變量P4 例1.1 球的體積計(jì)算表達(dá)式 分號(hào)(;),逗號(hào)(, ),省略號(hào)(.)。歷史指令調(diào)用科學(xué)記數(shù)法數(shù)據(jù)顯示格式Shortlongrational顯示格式與計(jì)數(shù)精度區(qū)別2022/7/145第一章 Matlab入門1.2 數(shù)據(jù)和變量復(fù)數(shù)i, j注意(-8)(1/3)預(yù)定義變量pi 圓周率3.1415eps 浮點(diǎn)數(shù)識(shí)別精度2.2210-16realmin 最小正實(shí)數(shù)2.225110 -308realmax 最大正實(shí)數(shù)1.797710308 Inf 無窮大

3、 NaN 不定值注意e在MATLAB中只是普通變量 2022/7/146第一章 Matlab入門1.2 數(shù)據(jù)和變量用戶變量命名規(guī)則：字母開頭，由字母、數(shù)字或下劃線組成（不得使用減號(hào)、空格等 ），區(qū)分大小寫防止與系統(tǒng)的預(yù)定義變量名(如i, j, pi, eps等)，函數(shù)名(如who, length等)，保留字(for, if , while, end等)沖突。特殊變量ans 是系統(tǒng)本身一個(gè)特殊變量名，若運(yùn)算結(jié)果沒有賦于任何變量，系統(tǒng)將其賦予ans clear 清除(注意Clear Workspace與Clear Command Window的區(qū)別. )2022/7/147第一章 Matlab入門

4、1.2 數(shù)據(jù)和變量數(shù)據(jù)Mat文件實(shí)現(xiàn)與外部數(shù)據(jù)文件交換: mat, txt等菜單方式：Save Workspace和Import Data例 : save-clear-import指令方式：save 和 loadMATLAB還提供了方便的工具來導(dǎo)入外部數(shù)據(jù)文件，包括文本文件、Excel(Spreadsheets)文件、圖像文件等，以便與其它應(yīng)用程序交換數(shù)據(jù)，我們將在第二章介紹。 2022/7/148第一章 Matlab入門1.3 數(shù)組及其運(yùn)算a=1 2 3;4 5 6;7 8 0變量編輯器(Variables)和剪貼板的使用數(shù)組的輸入和分析 中括號(hào) 表示矩陣，同行無素間用空格或逗號(hào)分隔，不同行

5、間用分號(hào)或回車分隔。冒號(hào)運(yùn)算函數(shù)linspace(x1,x2,n) 生成x1與x2間的n維等距行向量,即將x1,x2 n-1等分length, size編址：不能為0，按列編址, 如a(6)2022/7/149第一章 Matlab入門1.3 數(shù)組及其運(yùn)算數(shù)組的抽取和統(tǒng)計(jì) 分塊矩陣a(1,3,1:3)sum,prod,min,max按列計(jì)算例如 x,y=max(-1 3;5 -6;2 4)x = 5 4y = 2 32022/7/1410第一章 Matlab入門1.3 數(shù)組及其運(yùn)算數(shù)組運(yùn)算 A+B與A-B 加與減 k*A或A*k 數(shù)乘矩陣 k+A與k-A k加(減)A的每個(gè)元素A.k , k.A

6、 數(shù)組乘方 A.*B 數(shù)組乘數(shù)組 k./A 數(shù)除以數(shù)組 左除A.B=右除B/.A 數(shù)組除法 點(diǎn)運(yùn)算就是對應(yīng)元素的運(yùn)算P12例（注意點(diǎn)運(yùn)算與矩陣運(yùn)算的區(qū)別）2022/7/1411第一章 Matlab入門1.3 數(shù)組及其運(yùn)算數(shù)學(xué)函數(shù)矩陣的數(shù)學(xué)函數(shù)也是按元素的運(yùn)算，使用通常的函數(shù)號(hào)，如sin(A), cos(A), asin(A) , acos(A), tan(A),cot(A), exp(A), sqrt(A)等。fix 向0取整 floor 向下取整ceil 向+取整 mod 模余rem 除法余數(shù) abs 絕對值(模）real 復(fù)數(shù)實(shí)部 imag 復(fù)數(shù)虛部angle 復(fù)數(shù)幅角 conj 復(fù)數(shù)共軛

7、log 自然對數(shù)ln log10 以10為底對數(shù)2022/7/1412第一章 Matlab入門1.3 數(shù)組及其運(yùn)算4、關(guān)系與邏輯運(yùn)算 、 、 = 大于、大于等于 = = 、 = 等于、不等于&（與）、|（或）、 （非）any、all 、find“真(true)”用1表示，“假(false)”用0.邏輯運(yùn)算中，所有非零元素作為1處理.P15例子小結(jié)：Matlab中的三類運(yùn)算函數(shù)按元素運(yùn)算（數(shù)學(xué)函數(shù)）：如sin(A), log(A), find(A)按列運(yùn)算（統(tǒng)計(jì)函數(shù)）：如sum(A), prod(A), min(A), max(A), any(A), all(A)按矩陣運(yùn)算（矩陣函數(shù)）：如siz

8、e(A), length(A), det(A)2022/7/1414第一章 Matlab入門1.4 字符串、元胞和結(jié)構(gòu)數(shù)據(jù)類型：數(shù)值(Double)邏輯(Logical)字符(Char)元胞(Cell)結(jié)構(gòu)(Structure)2022/7/1415第一章 Matlab入門1.4 字符串、元胞和結(jié)構(gòu)字符串（P17例）單引號(hào)(英文半角輸入狀態(tài)！)中文字符不要在word中輸入后copy引號(hào)內(nèi)字符顯示應(yīng)為淡紫色字符串拼接字符串轉(zhuǎn)化double, char, num2str, str2num比較：a=12,b=double(a),c=str2num(a)Eval執(zhí)行字符串書寫的指令(P18例)2022

9、/7/1416第一章 Matlab入門1.4 字符串、元胞和結(jié)構(gòu)元胞和結(jié)構(gòu)（P18例） 字符矩陣char數(shù)值與字符混合元胞 結(jié)構(gòu): 域的概念struct2cell和cell2struct 2022/7/1417第一章 Matlab入門習(xí)題P20 ex1, ex2, ex3, ex4MATLAB數(shù)學(xué)實(shí)驗(yàn)第二章 MATLAB編程與作圖2022/7/1419第二章 MATLAB編程與作圖第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì)2.2 作圖 2.3 在線幫助和文件管理2022/7/1420第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì)循環(huán)語句for 循環(huán)變量=初值:增量:終值, 語句；end w

10、hile（條件式）, 語句；end分支語句if（條件式）, 語句; endif（條件式1), 語句1; elseif （條件式2), 語句 2; ;else, 語句;endswitch(分支變量)case(值1), 語句1;case(值2), 語句2;otherwise 語句;end其它：pause, break, return, error2022/7/1421第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì) s=0; for n=1:100 , s=s+1/n/n;end；s clear;s=0;n=1; while nm=1000;str= sum(1./(1:m).2;eval(str)

11、2022/7/1426第二章 MATLAB編程與作圖inline函數(shù)與匿名函數(shù)Inline函數(shù)(即將被淘汰)fun=inline(函數(shù)表達(dá)式, 自變量) fname=inline(sum(1./(1:m).2),m) fname(1000) 匿名函數(shù)fun=(自變量)函數(shù)表達(dá)式匿名函數(shù)的參數(shù)傳遞更靈活 k=2; fname=(m)sum(1./(1:m).k) fname(1000)2022/7/1427第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì)注釋：%開頭，對本行后面字符起作用，不參與運(yùn)算。對話：input，disp 全程變量與局部變量 nargin、nargout和varargin子函

12、數(shù)和嵌套函數(shù)提高速度2022/7/1428第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì)普通編程function s=f(m)s=0;for n=1:m s=s+1/n/n;end向量化編程function s=f(m)n=1:m;s=sum(1./n.2);盡量少用for語句2022/7/1429第二章 MATLAB編程與作圖本書配套工具箱的使用主要內(nèi)容教材所有程序；數(shù)學(xué)實(shí)驗(yàn)常用Matlab指令的中文幫助安裝步驟1. 將ME.rar解壓縮至matlabtoolbox; 2. 啟動(dòng)Matlab,利用File菜單中的Set path將 matlabtoolboxme增至path中, 放在最底部，

13、并保存設(shè)置； 3. 回到你的工作目錄?，F(xiàn)在ME已成為一個(gè)普通的工具箱了。2022/7/1430第二章 MATLAB編程與作圖2.1 程序設(shè)計(jì)例2.4 編一M函數(shù),對任意輸入的向量x, 可計(jì)算分段函數(shù)值構(gòu)成的向量。分量方式 eg2_4a，慢向量方式 eg2_4b, eg2_4c, 快 數(shù)組預(yù)分配y=zeros(size(x)2022/7/1431第二章 MATLAB編程與作圖2.2 作圖曲線圖 plot(x,y) 以數(shù)據(jù)(x(i), y(i)為節(jié)點(diǎn)的折線圖, 其中x, y為同長度的向量 plot(x1,y1,x2,y2,.) 多組數(shù)據(jù)折線圖 fplot(fun,a,b) 函數(shù)fun在區(qū)間a,b

14、上的函數(shù)圖 plot3(x,y,z) 空間曲線圖, 其中x, y, z為同長度的向量圖形導(dǎo)出到word線型與標(biāo)記 P31 表 eg2_5 曲線圖y=x3-x-1和y=|x|0.2sin(5x) 2022/7/1432第二章 MATLAB編程與作圖2.2 作圖曲面圖x,y=meshgrid(xa,ya)當(dāng)xa, ya分別為m維和n維行向量，得到x和y均為n行m列矩陣。meshgrid常用于生成X-Y平面上的網(wǎng)格數(shù)據(jù)。x,y=ndgrid(xa,ya)與meshgrid類似，但得到的x和y均為m行n列矩陣 mesh(x,y,z) 繪制網(wǎng)面圖，是最基本的曲面圖形命令, 其中x, y, z是同階矩陣，

15、表示曲面三維數(shù)據(jù)。surf(x,y,z) 繪制曲面圖，與mesh用法類似。2022/7/1433第二章 MATLAB編程與作圖meshgrid解釋xa=6:8;ya=1:4; x,y=meshgrid(xa,ya) %生成X-Y面上網(wǎng)格 z=x.2+y.2 %計(jì)算XY面各格點(diǎn)的z軸高度 x y z6 7 8 1 1 1 37 50 656 7 8 2 2 2 40 53 686 7 8 3 3 3 45 58 736 7 8 4 4 4 52 65 802022/7/1434第二章 MATLAB編程與作圖ngrid解釋xa=6:8;ya=1:4; x,y=ngrid(xa,ya) %生成X-Y

16、面上網(wǎng)格 z=x.2+y.2 %計(jì)算XY面各格點(diǎn)的z軸高度 x y z6 6 6 6 1 2 3 4 37 40 45 527 7 7 7 1 2 3 4 50 53 58 658 8 8 8 1 2 3 4 65 68 73 80eg2_6 二元函數(shù)圖 z = xexp(-x2-y2 )2022/7/1435第二章 MATLAB編程與作圖2022/7/1436第二章 MATLAB編程與作圖2.2 作圖圖形說明和定制title 標(biāo)題說明；xlabel，ylabel，zlabel 說明坐標(biāo)軸x,y,z;hold on/hold off 保留/釋放現(xiàn)有圖形axis(a,b,c,d) 確定坐標(biāo)軸范圍

17、axb,cyGrade,Name,Raw=xlsread(c:grade.xls,Sheet1,a2:d3) Total= Grade (:,1)+ Grade (:,2)+ Grade (:,3); xlswrite(c:grade.xls,Total,Sheet1,e2:e3) xlswrite(c:grade.xls,總分,Sheet1,e1) 2022/7/1442第二章 MATLAB編程與作圖習(xí)題P46 ex2, ex3, ex5, ex6(單號(hào))MATLAB數(shù)學(xué)實(shí)驗(yàn)第三章 矩陣代數(shù) 2022/7/1444第三章 矩陣代數(shù)第三章 矩陣代數(shù)3.1 預(yù)備知識(shí)：線性代數(shù)3.2 矩陣代數(shù)的M

18、ATLAB指令3.3 計(jì)算實(shí)驗(yàn)：線性方程組的通解3.4 建模實(shí)驗(yàn)：投入產(chǎn)出分析和基因遺傳2022/7/1445第三章 矩陣代數(shù)3.1 預(yù)備知識(shí)：線性代數(shù)線性方程組記為 A x = b2022/7/1446第三章 矩陣代數(shù)3.1 預(yù)備知識(shí)：線性代數(shù)線性方程組若秩(A) 秩(A,b)，則無解；若秩(A) = 秩(A,b) = n, 存在唯一解；若秩(A) = 秩(A,b) A=1 1/4 0;0 1/2 0;0 1/4 1;P,T=eig(A) AP=PT P-1AP=T A=PTP-1 An=PTP-1PTP-1PTP-1=PTnP-12022/7/1460第三章 矩陣代數(shù)3.4 建模實(shí)驗(yàn)設(shè)有n

19、個(gè)經(jīng)濟(jì)部門，xi為部門i的總產(chǎn)出，cij為部門j單位產(chǎn)品對部門i產(chǎn)品的消耗，di為外部對部門i的需求，fj為部門j新創(chuàng)造的價(jià)值。分配平衡方程組(部門i產(chǎn)品=內(nèi)部需求+外部需求)消耗平衡方程組(部門j產(chǎn)值=生產(chǎn)消耗+新創(chuàng)造價(jià)值) 2022/7/1461第三章 矩陣代數(shù)列昂杰夫 Leontief, Nobel經(jīng)濟(jì)學(xué)獎(jiǎng) 1973年W. W. Leontief was a Russian-American economist notable for his research on how changes in one economic sector may have an effect on othe

20、r sectors. Leontief won the Nobel Memorial Prize in Economic Sciences in 1973, and three of his doctoral students have also been awarded the prize (Paul Samuelson 1970, Robert Solow 1987, Vernon Smith 2002).2022/7/1462第三章 矩陣代數(shù)投入產(chǎn)出分析令 C =（cij），X = (x1, , xn),D = (d1, , dn)，F(xiàn)= (f1, , fn), 則分配平衡方程組 X=C

21、X+D令 A = EC，E為單位矩陣，則 AX = DC稱為直接消耗矩陣A稱為列昂杰夫矩陣。2022/7/1463第三章 矩陣代數(shù)Y = 1,1,1 B （B各列的和）Y表示各部門的總投入(消耗)。新創(chuàng)造價(jià)值 F=X YB=CB表示各部門間的投入產(chǎn)出關(guān)系，稱為投入產(chǎn)出矩陣。注：bij=cijxj2022/7/1464第三章 矩陣代數(shù)投入產(chǎn)出關(guān)系表(行：分配平衡，列：消耗平衡)消耗部門外界需求總產(chǎn)出123生產(chǎn)部門1b11b12b13d1x12b21b22b23d2x23b31b32b33d3x3新創(chuàng)造價(jià)值f1f2f3總產(chǎn)出x1x2x3注：bij=cijxj2022/7/1465第三章 矩陣代數(shù)投

22、入產(chǎn)出平衡行：分配平衡X=B各行之和+外界需求列：消耗平衡X=B各列之和+新創(chuàng)造價(jià)值2022/7/1466第三章 矩陣代數(shù)投入產(chǎn)出分析 例3.4 某地有三個(gè)產(chǎn)業(yè)，一個(gè)煤礦，一個(gè)發(fā)電廠和一條鐵路，開采一元錢的煤，煤礦要支付0.25元的電費(fèi)及0.25元的運(yùn)輸費(fèi); 生產(chǎn)一元錢的電力，發(fā)電廠要支付0.65元的煤費(fèi)，0.05元的電費(fèi)及0.05元的運(yùn)輸費(fèi); 創(chuàng)收一元錢的運(yùn)輸費(fèi),鐵路要支付0.55元的煤費(fèi)和0.10元的電費(fèi)，在某一周內(nèi)煤礦接到外地金額50000元定貨，發(fā)電廠接到外地金額25000元定貨，外界對地方鐵路沒有需求。2022/7/1467第三章 矩陣代數(shù)解：這是一個(gè)投入產(chǎn)出分析問題。設(shè)x1為本周內(nèi)

23、煤礦總產(chǎn)值，x2為電廠總產(chǎn)值, x3為鐵路總產(chǎn)值, 則問三個(gè)企業(yè)間一周內(nèi)總產(chǎn)值多少才能滿足自身及外界需求？三個(gè)企業(yè)間相互支付多少金額？三個(gè)企業(yè)各創(chuàng)造多少新價(jià)值?2022/7/1468第三章 矩陣代數(shù)直接消耗矩陣C= 外界需求向量 D =產(chǎn)出向量X = 則原方程為 (E-C)X=D 投入產(chǎn)出矩陣為 B=C*diag(X)總投入向量 Y= ones(1,3)*B 新創(chuàng)造價(jià)值向量 F=X-Y2022/7/1469第三章 矩陣代數(shù)表3.3 投入產(chǎn)出分析表(單位：元)消耗部門外界需求總產(chǎn)出煤 礦電 廠鐵 路生產(chǎn)部門煤礦0365061558250000102088電廠2552228082833250005

24、6163鐵路2552228080028330新創(chuàng)造價(jià)值51044140419915總產(chǎn)出10208856163283302022/7/1470第三章 矩陣代數(shù)例5 設(shè)金魚某種遺傳病染色體的正?；?yàn)锳，異?；?yàn)閍, 那么AA，Aa，aa分別表示正常金魚，隱性患者，顯性患者。設(shè)初始分布為90%正常金魚，10%的隱性患者，無顯性患者?？紤]下列兩種配種方案對后代該遺傳病基因型分布的影響方案一：同類基因結(jié)合，均可繁殖；方案二：顯性患者不允許繁殖，隱性患者必須與正常金魚結(jié)合繁殖2022/7/1471第三章 矩陣代數(shù)后代基因型的概率后代是從父母體的基因?qū)χ懈骼^承一個(gè)基因，形成自己的基因型。2022/7/

25、1472第三章 矩陣代數(shù)第一種方案同類基因結(jié)合，均可繁殖；2022/7/1473第三章 矩陣代數(shù)第二種方案顯性患者不允許繁殖，隱性患者必須與正常金魚結(jié)合繁殖2022/7/1474第三章 矩陣代數(shù)基因遺傳矩陣表示2022/7/1475第三章 矩陣代數(shù)解 設(shè)初始分布X(1)=(0.9 0.1 0),第n代分布為X(n)=方案一： X(n) = Mn-1X(1)方案二： X(n) = Nn-1X(1)2022/7/1476第三章 矩陣代數(shù)結(jié)論計(jì)算方法：取n足夠大，如n=20方案一：Mn-1X(1) 0.95, 0, 0.05存在5%的顯性患者方案二： Nn-1X(1) 1, 0, 0顯性患者和隱性患

26、者都消失本例說明了雜交的優(yōu)勢另一計(jì)算方法：用矩陣相似對角化2022/7/1477第三章 矩陣代數(shù)習(xí)題P65 ex2, ex3, ex4, ex5, ex6, ex102022/7/1478第三章 矩陣代數(shù)習(xí)題4 轉(zhuǎn)移矩陣表示2022/7/1479第三章 矩陣代數(shù)習(xí)題5 投入產(chǎn)出關(guān)系消耗部門外界需求總產(chǎn)出123生產(chǎn)部門1c11x1c12x2c13x3d1x12c21x1c22x2c23x3d2x23c31x1c32x2c33x3d3x3注：bij=cijxjMATLAB數(shù)學(xué)實(shí)驗(yàn)第四章 函數(shù)和方程第四章 函數(shù)和方程4.1 預(yù)備知識(shí)：零點(diǎn)、極值和最小二乘法4.2 函數(shù)零點(diǎn)、極值和最小二乘擬合的MAT

27、LAB指令4.3 計(jì)算實(shí)驗(yàn)：迭代法4.4 建模實(shí)驗(yàn)：購房貸款的利率4.1 預(yù)備知識(shí)：零點(diǎn)非線性方程 f (x) = 0若對于數(shù)有f () = 0, 則稱為方程的解或根，也稱為函數(shù)f (x)的零點(diǎn)非線性方程求解通常用數(shù)值方法求近似解.非線性方程（組）f (x) = 0，其中x=(x1, x2, , xn), f=(f1, f2, , fm) 4.1 預(yù)備知識(shí)：極值設(shè)x為標(biāo)量或向量，y=f(x)是xD上的標(biāo)量值函數(shù)。如果對于包含x=a的某個(gè)鄰域 ,有 f(a)f(x) (f(a)f(x)對任意x成立, 則稱a為f(x)的一個(gè)局部極小(大)值點(diǎn)。如果對任意xD，有f(a)f(x)（f(a)f(x)）

28、成立，則稱a為f(x)在區(qū)域D上的一個(gè)全局極小(大)值點(diǎn)。4.1 預(yù)備知識(shí)：極值4.1 預(yù)備知識(shí)：最小二乘擬合4.1 預(yù)備知識(shí)：最小二乘擬合成績選手國籍日期13秒24米爾布恩美國1972年9月7日13秒21卡薩那斯古巴1977年8月21日13秒16內(nèi)赫米亞美國1979年4月14日13秒00內(nèi)赫米亞美國1979年5月6日12秒93內(nèi)赫米亞美國1981年8月19日12秒92金多姆美國1989年8月16日12秒91杰克遜英國1993年8月20日12秒88劉翔中國2006年7月12日12秒87羅伯斯古巴2008年6月12日12秒80梅里特美國2012年9月7日數(shù)據(jù)：男子110米欄記錄問題：估計(jì)什么時(shí)候

29、突破12秒50？4.1 預(yù)備知識(shí)：最小二乘擬合4.1 預(yù)備知識(shí)：最小二乘擬合假設(shè)已知經(jīng)驗(yàn)公式y(tǒng)=f(c,x)(c為參數(shù), x為自變量), 要求根據(jù)一批有誤差的數(shù)據(jù)(xi,yi), i=0,1,n, 確定參數(shù)c. 這樣的問題稱為數(shù)據(jù)擬合。最小二乘法：求c使得平方誤差最小化 Q(c)=4.2 函數(shù)零點(diǎn)MATLAB指令多項(xiàng)式y(tǒng)=polyval(p,x) 求得多項(xiàng)式p在x處的值y，x可以是一個(gè)或多個(gè)點(diǎn)p3=conv(p1,p2) 返回多項(xiàng)式p1和p2的乘積p3,r=deconv(p1,p2) p3返回多項(xiàng)式p1除以p2的商，r返回余項(xiàng)x=roots(p) 求得多項(xiàng)式p的所有復(fù)根.p=polyfit(x

30、,y,k)用k次多項(xiàng)式擬合向量數(shù)據(jù)(x, y),返回多項(xiàng)式的降冪系數(shù)MATLAB中一個(gè)多項(xiàng)式用系數(shù)降冪排列向量來表示。例2.用2次多項(xiàng)式擬合下列數(shù)據(jù). x 0.1 0.2 0.15 0 -0.2 0.3 y 0.95 0.84 0.86 1.06 1.50 0.72 M文件eg4_2.m例1.求多項(xiàng)式x3 + 2 x2 - 5的根 p=1 2 0 -5; x=roots(p) , polyval(p,x) Fun=Mfun 定義一個(gè)函數(shù)句柄，這里Mfun是函數(shù)的M文件表達(dá)方式Fun=(var)funstr 定義匿名函數(shù)，其中var是變量名, funstr是函數(shù)的表達(dá)式 非線性函數(shù)的MATLAB

31、表達(dá) x=fzero(Fun, x0) 返回一元函數(shù)Fun的一個(gè)零點(diǎn). x0為標(biāo)量時(shí), x 返回函數(shù)在x0附近的零點(diǎn)； x0為向量a, b時(shí), x 返回在a,b中的零點(diǎn)(兩端函數(shù)異號(hào)) 4.2 函數(shù)零點(diǎn)MATLAB指令 x,f,h=fsolve(Fun, x0) x: 返回多元函數(shù)Fun在x0附近的一個(gè)零點(diǎn)，其中x, x0均為向量; f: 返回Fun在零點(diǎn)的函數(shù)值, 應(yīng)該接近0; h: 返回值如果大于零，說明計(jì)算結(jié)果可靠，否則計(jì)算結(jié)果不可靠。 例3 求函數(shù)y=xsin(x2-x-1)在(-2, -0.1)內(nèi)的 零點(diǎn) fun=(x)x*sin(x2-x-1) fzero(fun,-2 -0.1)

32、 fplot(fun,-2,-0.1),grid on; fzero(fun,-2,-1.2), fzero(fun,-1.2,-0.1) fzero(fun,-1.6), fzero(fun,-0.6) 例4 求方程組在原點(diǎn)附近的解xx(1)y x(2)function f=eg4_4fun(x)f(1)=4*x(1)-x(2)+exp(x(1)/10-1;f(2)=-x(1)+4*x(2)+x(1)2/8; x,f,h=fsolve(eg4_4fun,0 0) 使用函數(shù)句柄%使用匿名函數(shù) fun=(x)4*x(1)-x(2)+exp(x(1)/10-1, -x(1) +4*x(2)+x(1

33、)2/8; x,f,h=fsolve(fun,0 0) min(y) 返回向量y的最小值max(y) 返回向量y的最大值x,f=fminbnd(fun,a,b) x返回一元函數(shù)fun在a,b內(nèi)的 局部極小值點(diǎn), f返回局部極小值x,f=fminsearch(fun,x0) x返回多元函數(shù)fun在初始值x0 附近的局部極小值點(diǎn), f返回局部極小值. x, x0均為向量。4.2 函數(shù)極值MATLAB指令例 5 .求二元函數(shù)f(x,y)= 5-x4-y4+4xy在原點(diǎn)附近的極大值。解：max fmin(-f) x x(1), y x(2) fun=(x)x(1)4+x(2)4-4*x(1)*x(2)

34、-5; x,g=fminsearch(fun,0,0)注：在使用fsolve, fminsearch等指令時(shí)， 多變量必須合寫成一個(gè)向量變量，如用x(1), x(2),。 4.2 最小二乘擬合MATLAB指令假設(shè)已知經(jīng)驗(yàn)公式y(tǒng)=f(c,x)(c為參數(shù), x為自變量), 要求根據(jù)一批有誤差的數(shù)據(jù)(xi,yi), i=0,1,n, 確定參數(shù)c.這樣的問題稱為數(shù)據(jù)擬合。最小二乘法就是求c使得平方誤差最小化 Q(c)=基本格式: c=lsqcurvefit(Fun,c0, x, y) 這里Fun(c,x)為兩個(gè)輸入變量的函數(shù)句柄或匿名函數(shù)，c0為參數(shù)c的預(yù)估值，作為迭代初值，x, y為數(shù)據(jù)向量 完整格

35、式：c,Q=lsqcurvefit(Fun,c0, x, y, lb,ub) lb和ub分別表示c的下界和上界。c返回參數(shù)值，Q返回誤差平方和。自變量x可以是多變量，這時(shí)第三輸入?yún)?shù)x應(yīng)為矩陣。 最小二乘擬合函數(shù)類型的選擇根據(jù)物理或數(shù)學(xué)機(jī)理作圖看趨勢常用直線、二次函數(shù)，指數(shù)函數(shù)等參數(shù)初始值的選擇根據(jù)參數(shù)的實(shí)際意義由數(shù)據(jù)估計(jì)網(wǎng)格化搜索常用原點(diǎn)男子110欄問題Matlab程序clear;close all;data=13.24 13.21 13.16 13.00 12.93 12.92 12.91 12.88 12.87 12.80;year=1972 1977 1979 1979 1981 19

36、89 1993 2006 2008 2012;plot(year,data,o);y=data;x=year-1972;fun=(c,x)c(1)*exp(-c(2)*x);c=lsqcurvefit(fun,13,(13.24-12.80)/40,x,y)fun2=(x)c(1)*exp(-c(2)*x)-12.5;x2=fsolve(fun2,50)%初估50年后， %計(jì)算得73年以后，即2045年xx=0:75;yy=fun(c,xx);figure;plot(year,data,o,xx+1972,yy);grid on;男子110欄問題的求解模型估計(jì)大約2045年突破12秒50 迭代

37、法是從解的初始近似值x0(簡稱初值)開始，利用某種迭代格式x k+1 = g (x k ),求得一近似值序列x1, x2, , xk, xk+1, 逐步逼近于所求的解(稱為不動(dòng)點(diǎn))。最常用的迭代法是牛頓迭代法，其迭代格式為1 迭代法4.3 計(jì)算實(shí)驗(yàn)：迭代法Newton法幾何意義：切線法切線代替曲線牛頓法程序newton.mfunction x=newton(fname,dfname,x0,e)if nargine x0=x; x=x0-fname(x0)/dfname(x0);end例6 求方程 x 2 - 3 x + e x = 2 的正根 (要求精度 = 10 -6)解 令f (x) =

38、x 2 - 3 x + e x - 2, f(0)=-1,f (2) 0, 所以根在0,2內(nèi)。M文件eg4_6先用圖解法找初值，再用牛頓法程序newton.m求解。2線性化擬合例7 .用函數(shù)y=aebx 擬合例2的數(shù)據(jù)方法一：非線性擬合，記a=c(1), b=c(2) fun = (c,x)c(1)*exp(c(2)*x; x=;y=; c,Q=lsqcurvefit(fun,0,0,x,y)方法二：線性化擬合，兩邊取對數(shù)得z=lny=lna+bx轉(zhuǎn)化為線性擬合。 M文件eg4_7不難算出，你向銀行總共借了25.2萬，30年內(nèi)共要還51.696萬，這個(gè)案例中貸款年利率是多少呢？ 例8 .下面是

39、新民晚報(bào)2000年3月30日上的一則房產(chǎn)廣告：4.4 建模實(shí)驗(yàn)：購房貸款的利率解 設(shè)xk為第k個(gè)月的欠款數(shù)， a為月還款數(shù), r為月利率。xk+1 = (1+r) xk- a那么 xk = (1+r) xk-1- a = (1+r)2 xk-2 (1+r)a a = = (1+r)k x0 a1+(1+r)+(1+r)k-1 = (1+r)k x0 a(1+r)k-1/r根據(jù) a=0.1436, x0=25.2, x360=0得到 25.2(1+r)360 0.1436(1+r)360-1/r=0很難用roots求解！常識(shí)上，r應(yīng)比當(dāng)時(shí)活期存款月利率略高一些。我們用活期存款月利率0.0198/

40、12 作為迭代初值，用fzero求解clear; fun=(r)25.2*(1+r)360-(1+r)360- 1)/r*0.1436;r=fzero(fun,0.0198/12);R=12*r得年利率為5.53%. (你知道最新利率嗎？)分期付款月還款公式 方程 (1+r)N x0 a(1+r)N-1/r=0月還款計(jì)算x0剩余借款額；N剩余月數(shù)；r月利率=年利率/12線性迭代 xk+1=axk+b 收斂：當(dāng)|a|1, 趨于無窮大。4.4 建模實(shí)驗(yàn)：混沌不收斂的非線性迭代可能會(huì)趨于無窮大，趨于一個(gè)周期解，在一個(gè)有限區(qū)域內(nèi)雜亂無章地游蕩，這類由確定性運(yùn)動(dòng)導(dǎo)致的貌似隨機(jī)的現(xiàn)象稱為混沌現(xiàn)象4.4 建

41、模實(shí)驗(yàn)：混沌例（混沌）：昆蟲數(shù)量的Logistic模型xk+1 = a x k (1 - x k), 0a4 xk表示第k代昆蟲數(shù)量(1表示理想資源環(huán)境最大可能昆蟲數(shù)量)。 a為資源系數(shù).0a4保證了xk在區(qū)間(0,1)上封閉。當(dāng)0a1, 在0，1內(nèi)有一個(gè)不動(dòng)點(diǎn)0， 且由|g(0)| =a1, 不動(dòng)點(diǎn)0不再穩(wěn)定;當(dāng)1a3，由|g (1-1/a)| =|2-a|3，出現(xiàn)兩個(gè)周期2解，且3aa 迭代序列幾乎雜亂無章， 即所謂混沌。 分岔圖程序clear;close;a=0:0.01:4;M=length(a);K=1000;x=zeros(K,M);x(1,:)=rand(1,M);for m=1

42、:M, for k=1:K-1 x(k+1,m)=a(m)*x(k,m)*(1-x(k,m);end, endfor k=1:20, plot(a,x(k,:),.);title(k=,int2str(k);pause(1);end;plot(a,x(900:K,:),.);hold off;分岔圖*混沌的特征 (i)初值敏感性: 兩個(gè)任意近的點(diǎn)出發(fā)的兩條軌跡 遲早會(huì)分得很開;蝴蝶效應(yīng)(Lorenz 1963)：一只南美熱帶雨林中的蝴蝶，偶爾扇動(dòng)幾下翅膀，可能在兩周后在美國德克薩斯引起一場龍卷風(fēng)。 (ii)遍歷性: 任意點(diǎn)出發(fā)的軌跡總會(huì)進(jìn)入 0,1內(nèi)任意小的開區(qū)間。例9（蛛網(wǎng)圖）我們用蛛網(wǎng)圖來

43、顯示混沌的遍歷性。 yk = a x k (1 - x k), xk+1 = yk蛛網(wǎng)圖正好顯示迭代計(jì)算x0, y0, x1, y1,的一系列變化過程。 eg4_9.m蛛網(wǎng)圖M函數(shù) eg4_9.mfunction f=eg4_9(a,x0,m,n)x=0:0.01:1;y=a*x.*(1-x);plot(x,x,r,x,y,r);hold on;clear x,y;x(1)=x0;y(1)=a*x(1)*(1-x(1);x(2)=y(1);if mm, plot(x(i),x(i),x(i+1),y(i-1),y(i),y(i);endendhold off;習(xí)題ex1, ex2, ex6,

44、ex8(1)ex9, ex10(1)(3)ex12，ex13ex13圖thR10MATLAB數(shù)學(xué)實(shí)驗(yàn)第五章 應(yīng)用微積分第五章 應(yīng)用微積分5.1 預(yù)備知識(shí)：微積分的基本概念5.2 數(shù)值微積分MATLAB指令5.3 計(jì)算實(shí)驗(yàn)：數(shù)值微積分5.4 建模實(shí)驗(yàn)：奶油蛋糕、作案時(shí)間1.極限和連續(xù)數(shù)列極限: 0, N0 ，使當(dāng)nN時(shí) 有xn -a0,函數(shù)在x0點(diǎn)附近是上升的；當(dāng)f(x0)0,函數(shù)在x0點(diǎn)附近是下降的；當(dāng)f(x0)=0, x0為駐點(diǎn)（很可能是極值點(diǎn)）函數(shù)可導(dǎo)連續(xù)幾何意義：切線的斜率當(dāng)n=0 得微分中值定理 f(x) - f(x0) = f() (x- x0) 其中是x0與x之間的某個(gè)點(diǎn)。Tayl

45、or公式:當(dāng)f(x)在含有x0某個(gè)開區(qū)間內(nèi)具有直到n+1階的導(dǎo)數(shù)，Taylor展開的幾何解釋任意函數(shù)的局部多項(xiàng)式逼近3.多元函數(shù)微分學(xué) 設(shè)f(x,y)在點(diǎn)(x0,y0)附近有定義，當(dāng)(x,y)以任何方式趨向于(x0,y0)時(shí)，f(x,y)趨向于一個(gè)確定的常數(shù)A，則二重極限若 A=f(x0,y0), 稱f(x,y)在(x0,y0) 點(diǎn)連續(xù)f(x,y)在點(diǎn)(x0,y0)的偏導(dǎo)數(shù)多元函數(shù)Taylor公式 函數(shù)極值一元可微f(x)在內(nèi)點(diǎn)x0取得局部極值必要條件是f (x0) =0. 局部極大充分條件是f (x0) =0, f (x0) 0.函數(shù)極值二元可微f(x,y)在內(nèi)點(diǎn)(x0,y0) 取得局部極大

46、或極小的必要條件是梯度f (x0,y0)=0, 局部極大(或局部極小)充分條件是 f (x0,y0) =0且下列Hesse矩陣負(fù)定（或正定）。梯度f =( ), Hesse矩陣多元函數(shù)極值也有類似的結(jié)果。4 . 積分函數(shù)f(x)在區(qū)間a,b上的定積分其中 a=x0 x1xn=b, xi=xi-xi-1, i(xi-1,xi), i=1,2,n幾何意義：曲邊梯形的面積牛頓-萊布尼茲公式：若在a,b上, F(x)=f(x), 則二重積分定義曲線&曲面平面曲線(x(t), y(t), atb的長度為空間曲線(x(t), y(t), z(t), at z=quadl(x)exp(-x.2),-1,1)

47、 z= integral (x)exp(-x.2),-1,1) 注意：積分中fun里用數(shù)組運(yùn)算! fun = (x) exp(-x.2).*log(x).2; q = integral(fun,0,Inf) % integral能求反常積分 5. 重積分z=dblquad(Fun,a,b,c,d) 求得二元函數(shù) Fun(x,y) 在矩形區(qū)域的重積分.z=triplequad(Fun,a,b,c,d,e,f)求得三元函數(shù)Fun(x,y,z) 在長方體區(qū)域上的三重積分 。z=quad2d(Fun,a,b,cx,dx) 求得二元函數(shù)Fun(x,y)的一般區(qū)域重積分。a, b為變量x的下、上限；cx,

48、 dx為變量y的下、上限函數(shù)(自變量為x)z= integral2 (Fun,a,b,cx,dx) 類似quad2dz= integral3 (Fun,a,b,cx,dx,exy,fxy) 求得三元函數(shù)Fun(x,y,z)一般區(qū)域的三重積分 例5.2 計(jì)算重積分 fun = (z,x,y) log(3+x+y+z) + z.2.*cos(x);xmin = (z)-sqrt(1 - z.2);xmax = (z) sqrt(1 - z.2);ymin = (z,x)-sqrt(1 - x.2 - z.2);ymax = (z,x)sqrt(1 - x.2 - z.2);q = integral

49、3(fun,0,1,xmin,xmax,ymin,ymax) fun = (x,y) 1./( sqrt(x + y) .* (1 + x + y).2 );q = integral2(fun,0,1,0,(x)1-x) 1.數(shù)值微分若f(x)在x=a可導(dǎo), 設(shè)h0且足夠小稱為向前差商向后差商中心差商5.3 計(jì)算實(shí)驗(yàn)：數(shù)值微積分差商公式比較中心差商精度比較高自動(dòng)步長的中心差商程序：deriv.mfunction d=deriv(fname,a,h0,e)h=h0;d=(fname(a+h)-fname (a-h)/2/h;d0=d+2*e; %為了保證abs(d-d0)e成立while abs

50、(d-d0)e d0=d;h0=h;h=h0/2 d=(fname(a+h)-fname (a-h)/2/h;end deriv(x)sin(x),pi/6,0.1,1e-4)例 5.4 （奶油蛋糕）某數(shù)學(xué)家的學(xué)生要送一個(gè)特大的蛋糕來慶賀他90歲生日。為了紀(jì)念他提出的口腔醫(yī)學(xué)的懸鏈線模型，學(xué)生們要求蛋糕店老板將蛋糕邊緣半徑作成下列懸鏈線函數(shù) r = 2-(exp(2h)+exp(-2h)/5, 0h1(單位：米) 。問如何計(jì)算重量？5.4 建模實(shí)驗(yàn)：解 設(shè)高為H，半徑 r, 比重為k若蛋糕是單層圓盤的，則蛋糕的重量為： W=kHr2rHr1r2若蛋糕是雙層的，每層高H/2 ，下層半徑r1, 上

51、層半徑r2，則 W=kH(r12+ r22)/2如果蛋糕是n層的，每層高H/n，半徑分別r1, , rn, 則fun=(h)(2-(exp(2*h)+exp(-2*h)/5) ;n=5;h=0.1:0.2:0.9;r=fun(h);W=pi*1/n*sum(r.2)計(jì)算結(jié)果5.4327當(dāng)n，fun2=(h)(2-(exp(2*h)+exp(-2*h)/5).2;pi*integral(fun2,0,1)計(jì)算結(jié)果5.4171即旋轉(zhuǎn)體體積公式解：t=0時(shí)z=R. 考慮在時(shí)間dt內(nèi)水面變化dz，漏水的體積為 uAdt= - x2dz其中u為出水速度，A為小孔面積，x為高度z水面的半徑5m0.5m0例

52、5.5 （作案時(shí)間）鍋爐房中發(fā)生一起命案，球形鍋爐底部有個(gè)彈孔在漏水。鍋爐半徑R=5m，平時(shí)充滿了水，測得底部小孔半徑b=0.1m，此時(shí)水面已將下降至離底部0.5m。推測槍擊是多長時(shí)間以前發(fā)生的?分析：求時(shí)間t與水面高度z的對應(yīng)關(guān)系。 z表示從球心 測量的水面高度A=b2, x2=R2-z2注意：水從孔漏出的速度u不是常數(shù)！zxt=0時(shí)z=R. 在頂部水降到0.5m時(shí)，z=0.5-R，從而 t =0+5m0.5m0u由下列能量方程決定 mg(z+R)=mu2/2 m是質(zhì)量，g為重力加速度。%M腳本eg5_6.mclear;R=5;b=0.1;g=9.81;z1=0.5-R;z2=R;n=100

53、;h=(z2-z1)/n;z=z1:h:z2;f=(R2-z.2)./(b2*sqrt(2*g*(z+R);I=trapz(z,f)/60/60結(jié)果為0.5144小時(shí)。Page101 習(xí)題ex4ex5(2)(4)(6)(7)(8)ex6ex10ex12 ex13 MATLAB數(shù)學(xué)實(shí)驗(yàn) 第六章 常微分方程第六章 常微分方程6.1 預(yù)備知識(shí)：常微分方程6.2 解常微分方程的MATLAB指令6.3 計(jì)算實(shí)驗(yàn)：Euler法和剛性方程組6.4 建模實(shí)驗(yàn): 導(dǎo)彈系統(tǒng)的改進(jìn) 1.微分方程的概念常微分方程: f (t,y,y,y,y(n)=0常微分方程的解: 函數(shù)y(t).6.1 預(yù)備知識(shí)：常微分方程微分方程

54、組: 聯(lián)系一些未知函數(shù)x(t), y(t), z(t), 的一組微分方程.偏微分方程: 含有偏導(dǎo)數(shù)的微分方程. 其解為多元函數(shù)u(t,x,y,z)。6.1 預(yù)備知識(shí)：常微分方程2.常微分方程解析求解：某些特殊的方程(1)初等積分法例如 y=ty, y(0)=1.分離變量dy/y=tdt, lny=t2/2+c,由初始值得c=0, 所以y(t)=exp(t2/2)（2）線性常系數(shù)微分方程的特征根法線性方程: y(n) + a1 (t) y(n-1) + + an-1 (t) y + an (t) y = b(t)常系數(shù)方程: 若ai (t) (i =1, ,n) 與t無關(guān)。齊次方程: 若b(t)

55、=0。 y(n) + a1 y(n-1) + + an-1 y + an y = 0線性常系數(shù)齊次微分方程的解可用特征根法求得. n+a1 n-1+ + an-1 +an=0非齊次方程的解為一個(gè)特解和相應(yīng)的齊次方程通解的疊加。變系數(shù)方程可嘗試常數(shù)變易法。例6.1 求x+ 0.2 x+3.92x = 0的通解解 特征方程為2 + 0.2 +3.92=0 roots(1 0.2 3.92 求得共軛復(fù)根 i=-0.11.9774i, 通解為 x(t) = Aet cos(t) +Bet sin(t)= Ae-0.1t cos(1.9774t) +Be-0.1t sin(1.9774t)4.常微分方程

56、(組)初值問題數(shù)值解數(shù)值解法:尋求解y(t)在一系列離散節(jié)點(diǎn)t0 t1 tn odefun= (t,y)y-2*t/y ; t,y=ode45(odefun,0,4,1); t,y plot(t,y,o) ode45(odefun,0:1:4,1); t,y=ode45(odefun,0:1:4,1);t,y例2 解微分方程y -y+2t/y=0, y(0)=1, 0t4 例3 解微分方程組 0t t,x=ode45(eg6_3fun,0 30,1;0.5)編程器窗口指令窗口作圖 subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:); subplot(1,2,2);

57、plot(x(:,1),x(:,2);例4 求解微分方程組 已知當(dāng)=0時(shí), f=0, T=1, 解 首先引入輔助變量， t=, y1=f, , , y4=T, 化為一階方程組function f=eg6_4fun(t,y)f(1)=y(2);f(2)=y(3);f(3)=-3*y(1)*y(3)+2*y(2)2-y(4);f(4)=y(5);f(5)=-2.1*y(1)*y(5);f=f(:); y0=0, 0, 0.68, 1,-0.5; t,y=ode45(eg6_4fun,0 5,y0); plot(t,y(:,1),t,y(:,4),:);2. 邊值問題解法 常微分方程邊值問題 （1）

58、sinit=bvpinit(tinit,yinit) 由在粗略節(jié)點(diǎn)tinit的預(yù)估解yinit生成粗略解網(wǎng)絡(luò)sinit這里y是2維向量，表示解函數(shù)x(t)及其導(dǎo)函數(shù)x(t)。求解三部曲：粗網(wǎng)絡(luò)、解結(jié)構(gòu)、數(shù)值化其Matlab標(biāo)準(zhǔn)形式為 （2）sol=bvp5c(odefun,bcfun,sinit) 從粗網(wǎng)絡(luò)sinit產(chǎn)生解結(jié)構(gòu) odefun是微分方程組函數(shù) bcfun為邊值條件函數(shù) sol是一個(gè)結(jié)構(gòu), sol.x為求解節(jié)點(diǎn). sol.y是y(t)的數(shù)值解（3）sx=deval(sol,ti)由解結(jié)構(gòu)sol 計(jì)算在ti的數(shù)值解。例5 求解邊值問題解 首先改寫為標(biāo)準(zhǔn)形式。 令y(1)= z, y(

59、2)= z, 則方程為y(1)=y(2), y(2) = y(2)sin(y(1) 邊界條件為ya(1)=0, yb(1)+2=0eg6_5.m%根據(jù)z初始值預(yù)估：z=-1，z=0clear;close;sinit=bvpinit(0:4,-1;0) %-1;0是常數(shù)猜測值z=-1, z=0odefun=(t,y)y(2);-y(2)*sin(y(1);bcfun=(ya,yb)ya(1);yb(1)+2;sol=bvp5c(odefun,bcfun,sinit) %注意sol的域名t=linspace(0,4,101);y=deval(sol,t);plot(t,y(1,:),sol.x,s

60、ol.y(1,:),o,sinit.x,sinit.y(1,:),s)legend(解曲線,解點(diǎn),粗略解)1. Euler法Euler法：在節(jié)點(diǎn)處用差商近似代替導(dǎo)數(shù)Euler格式k=0,1,26.3 計(jì)算：Euler法和剛性方程組M函數(shù)euler.mfunction t, y = euler(odefun, tspan, y0, h)t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1, y(i+1)=y(i)+h*odefun(t(i),y(i);endt=t;y=y;odefun: f(t, y)的函數(shù)句柄或Inline函數(shù)，tspan=t0