2023/10/111Matlab在物理中的應用－案例研究2010年4月11日2023/10/112一、Matlab概述MATLAB是美國Mathworks公司推出的一套高性能的數值分析和計算軟件。最初的目的是為了解決矩陣的運算問題而開發的，經過不斷完善，已成為將矩陣運算、數值分析、圖形處理、編程技術結合在一起的功能強大的工具。它的優點在于快速開發計算方法，而不在于計算速度。2023/10/113Matlab的學習方法必須做大量的練習，熟悉其中的函數聯系和自己的課題相關，深刻的體會多看幫助文件，又一本好的參考書記住：Matlab可以做很多事情2023/10/1141.1Matlab的出現

70年代中期，CleveMoler和他的同事開發了LINPACK和EISPACK的Fortran子程序庫

70年代末期，CleveMoler在新墨西哥大學給學生開線性代數，為學生編寫了接口程序，這程序取名為MATLAB，即MATrixLABoratory1983年春天，工程師JohnLittle與Moler、SteveBangert一起開發了第二代專業版MATLAB1984年，MathWorks公司成立，MATLAB正是推向市場。2023/10/1151.2

Matlab的版本演化Matlab1.0Pcmatlab->matlab386Matlab3.5+simulinkMatlab4.0：simlink內嵌Matlab5.0:全面的面向對象Matlab5.1~5.3Matlab6.0Matlab6.5：購并了MATRIXxMatlab7.0:20042023/10/1161.3Matlab的特點1高效方便的矩陣和數組運算2編程效率高（庫函數非常豐富）3結構化、面向對象4方便的繪圖工具5用戶使用方便6功能強大的工具箱（核心部分中有數百個核心內部函數）總之，Matlab在信號處理、建模、系統識別、控制、優化等領域的問題處理上優勢明顯于其它高級語言。2023/10/1171.4

解決物理問題常用的軟件Mathmatica數學計算功能。復雜符號計算，主要在理論界MAGIC電磁粒子仿真軟件仿真軟件：HFSS、CST等等2023/10/118二、MATLAB基礎使用2.1MATLAB操作界面簡介操作桌面的缺省外貌與先前版本相比，6.x版引入大量的交互工作界面，如通用操作界面、幫助界面、演示界面等指令窗、歷史指令窗、工作空間瀏覽器和只有窗名的當前目錄窗“開始”按鈕(matlab6.5版本新加)通用操作界面2023/10/119通用交互界面介紹內存數組編輯器(ArrayEditor)在工作空間瀏覽器中對變量操作時才啟動交互界面分類目錄窗（LauchPad)。可展開成樹狀結構，列著所有交互界面。M文件編輯/調試器（Editor/Debugger)當進行“打開文件”操作時才啟動。幫助導航/瀏覽器(Helpavigator/Browser)

詳盡展示由超文本寫成的在線幫助。2023/10/11102.2Matlab命令的兩種執行方式一.交互式的命令執行方式用戶在命令窗口逐條輸入命令，matlab逐條執行，這種方式操作簡單直觀，但速度慢，中間過程無法保留。通用操作界面二.M命令文件的程序設計方式將有關命令編成程序存儲在一個文件(擴展名為.m)中，matlab自動依次執行，可調試復雜的程序(調試時只顯示第一個錯誤)，是實際應用中主要的執行方式2023/10/11112.3幫助的使用

MATLAB中有以下幾種方法獲得幫助：幫助命令、幫助窗口、MATLAB幫助臺、在線幫助頁或用Web菜單直接鏈接到MathWorks公司的主頁。2023/10/1112怎樣獲得詳盡的幫助信息首先是MATLAB窗口中的help命令，它最常用、最快，往往也是最有效的。給出關鍵字就能找到相應的說明。Help命令簡單易用，但是它要求準確給出關鍵字，如果記憶不準，就很難找到。這時就應該用第二種方案。HelpDesk：這是MATLAB的HTML格式的幫助文檔，有著很好的組織。可以在命令行鍵入helpdesk，進入“MatlabHelpDesk”。HelpDesk的搜索可以使用不完全關鍵字，這樣即使遇到記憶模糊的情況也可以很快查到。前面兩種方法所得到的幫助信息常常是不夠細致的，需要用第三種途徑來補充，這就是電子版的使用手冊。MATLAB里包含了大量關于MATLAB安裝、使用、編程以及各種工具箱等的電子版手冊，全部為PDF格式,具有搜索功能，因此也可以在指定的手冊中搜索關鍵字，使用非常方便。2023/10/11132.4Matlab變量變量要求變量名、函數名區分字母的大小寫。如SIN、Sin與sin是不同的對6.5以前的版本，變量名字符不超過31個，對6.5版本，變量名最多可包含63個字符。超出長度的字符被忽略變量必須以字母開頭，之后可以是任意字母、數字或者下滑線，但不能是空格和標點Matlab不支持漢字，漢字不能出現在變量名和文件名中2023/10/1114

特殊變量特殊變量取值ans用于結果的缺省變量名pi圓周率eps計算機的最小數inf無窮大如1/0Nan或nan不定量如0/0iji=j=realmin最小的可用正實數realmax最大的可用正實數2023/10/11152.5

復數不需要特殊的處理。復數可以表示為：>>a=10-9i>>a=10-9j%i=j=sqrt(-1)注意下面的例子中的表達式

>>a4=5+sin(.5)*i如果把sin(.5)*i中的*去掉，就會出錯，sin(.5)i對MATLAB是無意義的，直接以字符i和j結尾，只適用于簡單數值，不適用于表達式。

復數運算不需要特殊處理，可以直接進行。另外real(z),imag(z),abs(z)和angle(z)等2023/10/1116三、MATLAB的數組(矩陣)及其運算3.1數組(矩陣)的創建

1.直接輸入法整個輸入數組必須以方括號“[]”為其首尾數組的行與行之間必須用”;”或者【Enter】隔離數組元素必須由逗號”,”或空格分隔2.一維數組的冒號生成法x=a:inc:ba是數組的第一個元素，inc是間隔，b是最后一個元素a，inc，b之間的冒號須在英文狀態下輸入。中文中的冒號將導致matlab出錯。inc缺省值為1，可以取正數或負數。2023/10/1117語句生成矩陣3.linspace命令a=linspace(1,10,5);%產生線性等間距的列向量4.logspace命令b=logspace(0,2,10);%產生對數等間距的列向量5.矩陣連接c=[ab];6.數組生成函數

zerosoneseyerandrandn2023/10/11183.2矩陣運算矩陣的轉置用符號’來表示和實現非共軛轉置.’矩陣的四則運算+-*/\點乘：.*右除：./左除：.\矩陣的乘方^矩陣乘方.^元素對元素的乘方2023/10/1119四、流程控制為了便于應用，MATLAB提供了一些流程控制的命令。這些命令對腳本編寫帶來了一些方便，但是需要注意的是，盡量不要使用這些流程控制命令，尤其是循環控制命令。1for循環for循環允許一組命令以固定的次數重復，它的一般形式是forx=arraycommandendfor和end之間的命令串按數組array的每一列執行一次，直到n次后終止。2023/10/11201For循環應當注意的一些方面for循環不能使用內部重新賦值循環變量而終止；for循環內部接受任何有效的MATLAB數組；for循環可以嵌套；只要有矩陣形式可以解決的問題，不要使用for循環。使用for循環的算法執行很慢，一個好的MATLAB算法部應當出現循環語句。循環可以使用break跳出，但只跳出所在的循環，不跳出整個嵌套結構。2023/10/11212while循環與for循環以固定的次數求一組指令相反，while循環以不定的次數求一組語句的值。While循環的一般形式為：whileexpressioncommondsend只要表達式expression里的所有元素為真，就執行命令串commands。通常表達式求值給一個標量值，單數組值也同樣有效。2023/10/11223if-else-end結構很多情況下，命令的序列必須根據關系的檢驗有條件的執行，它由if-else-end結構提供。它的結構如下：ifexpression1commands1elseifexpression2commands2elseif…

…elseifcommandsend2023/10/1123五、數據的可視化二維圖形

在二維曲線的繪制命令中，最重要、最基本的命令時plot，其他許多特殊繪圖命令都以它為基礎形成。

1plot(x,’s’)x為實向量時，以該向量元素的下標為橫坐標、元素值為縱坐標畫出一條連續曲線。X為實矩陣時，則按列繪制每列元素值相對其下標的曲線，圖中曲線數等于X陣列數。s用來指定線型。（正方形）2023/10/11242plot(x,y,’s’)x,y為同維向量時，繪制以x、y元素為橫、縱坐標的曲線。3plot(x1,y1,’s1’,x2,y2,’s2’,…)

在此格式中，每個繪線“三元組”（x,y,’s’）的結構和作用，與plot(x,y,’s’)相同，不同的三元組之間沒有約束。2023/10/1125線型和標記

Plot命令中，在每一對數組數據后面，給plot一個附加參量，就可以指定所要的顏色和線型。

>>plot(x,y,’rd’,x,z,’k--’,x,y,’b-.’,x,z,’co’)b藍色,c青色，g綠色，k黑色，m紫紅色，r紅色，w白色，y黃色加格線和標注，gridon命令在當前圖形的單位標記處加格線。gridoff命令取消格線。用text命令可以在圖形的任意指定位置增加標記和其他文本信息，格式是text(x,y,’string’).2023/10/1126>>plot(x,y,x,z)>>grid>>title(‘sin(x)andcos(x)’)>>xlabel(‘independentvariablex’)>>ylabel(‘dependentvariableyandz’)>>text(2,0.8.’sin(x)’)2023/10/1127plot應用舉例——畫一個圓程序如下：R=2;t=0:pi/20:2*pi;x=R*sin(t);y=R*cos(t)plot(x,y,'b:'),grid2023/10/1128例1——無阻力拋射體的運動設目標相對于射點的高度為，給定初速，試計算物體在真空中飛行的時間和距離。◆建模無阻力拋射體的飛行是中學物理就解決了的問題，本題的不同點是目標和射點不在同一高度上，用MATLAB可使整個計算和繪圖過程自動化。其好處是快速地計算物體在不同初速和射角下的飛行時間和距離。關鍵是在求落點時間時，需要解一個二次線性代速方程。2023/10/1129MATLAB程序clear;y0=0;x0=0;

%初始位置vMag=input(‘輸入初始速度(m/s):’);%輸入初始速度的大小和方向vDir=input('輸入初速方向(度)：');

yf=input(‘輸入目標高度(m):’);%輸入目標高度vx0=vMag*cos(vDir*(pi/180));%計算x,y方向的初始速度vy0=vMag*sin(vDir*(pi/180));wy=-9.81;wx=0;%重力加速度（m/s^2)tf=roots([wy/2,vy0,y0-yf]);%解方程，計算落點時間tf=max(tf);%去掉庸解t=0:0.1:tf;y=y0+vy0*t+wy*t.^2/2;%計算軌跡x=x0+vx0*t+wx*t.^2/2;xf=max(x),plot(x,y);%計算射程grid2023/10/1130初始速度50m/s，初速方向為40度，目標高度為8m時的輸入結果2023/10/1131例2——有空氣阻力拋射體的運動步驟一：給定小球質量m,重力加速度g，拖曳系數C和小球面積A（半徑為r），根據D=rou*C*A/2求出D步驟二：給定小球的初始速度v（1）（大小和方向角theta），小球的位置（x（1）和y（1）。設定時間間隔detat步驟三：選定最大的間隔數N（或者最大的時間）2023/10/1132求解步驟步驟四：當n<N或者t<時循環步驟五到步驟九步驟五：計算加速度分量

步驟六：打印或者畫2023/10/1133步驟七：計算新的速度分量步驟八：計算新的位移坐標步驟九：時間上增加一個間隔步驟十：結束

2023/10/1134Matlab程序g=9.8;m=0.145;C=0.5;r=0.0366;A=pi*r^2;rou=1.2;D=(rou*C*A)/2;detat=0.01;theta=(35/180*pi);x(1)=0;y(1)=0;V(1)=50;Vx(1)=V(1)*cos(theta);Vy(1)=V(1)*sin(theta);N=700;2023/10/1135forn=1:NV(n)=sqrt(Vx(n)^2+Vy(n)^2)ax(n)=-(D/m)*V(n)*Vx(n);ay(n)=-g-(D/m)*V(n)*Vy(n);Vx(n+1)=Vx(n)+ax(n)*detat;Vy(n+1)=Vy(n)+ay(n)*detat;x(n+1)=x(n)+Vx(n)*detat+0.5*ax(n)*detat^2;y(n+1)=y(n)+Vy(n)*detat+0.5*ay(n)*detat^2;px(n)=Vx(1)*n*detat;py(n)=Vy(1)*n*detat-0.5*g*(n*detat)^2;endplot(x,y,'r',px,py,'g');gridxlabel('x'),ylabel('y');title('Projectilemotionwithairresistance');2023/10/1136結果2023/10/1137例3——有空氣阻力的落體運動。一個質量為75kg的人在某一高度靜止下落，空氣阻力F與v的關系滿足F=-k，其中常數k=0.22kg/m，根據牛頓第二定律知

若令g=9.80m/,求下落過程中速度v(t)和位移x(t)的表達式，并畫圖。2023/10/1138推導2023/10/1139Matlab程序v=dsolve('m*D2v+2*k*v*Dv=0','v(0)=0,Dv(0)=g','t')x=dsolve('m*g-k*(Dx)^2=m*D2x','x(0)=0,Dx(0)=0,D2x(0)=g','t')m=75;k=0.22;g=9.8t=0:0.5:30;v=tanh(t.*(g*m*k)^(1/2)/m)*(g*m*k)^(1/2)/k%x=(m/k)*log(cosh(t.*(k*g/m)^1/2))plot(t,v),gtext('v(t)')holdonplot(t,x),gtext('x(t)')2023/10/1140結果2023/10/1141例4－小球碰撞問題質量為m的小球以速度正面碰撞質量為M的靜止小球，