1、天津*大學本科生畢業論文應用MATLAB求解經典物理若干典型問題The application of MATLAB in solving some classical physics questions 專業班級：物理0701學生姓名：*指導教師：*學 院：理學院2011年5月摘 要MATLAB是 MathWorks公司推出的一套科學計算軟件，MATLAB的意思是矩陣實驗室。MATLAB具有起點低、功能強大、易學易用以及兼有數值運算和符號運算功能的優點。利用MATLAB，繪圖十分方便，它既可以繪制各種圖形，包括二維圖形和三維圖形，還可以對圖形進行修飾和控制。本文通過在MATLAB環境下編寫通過

2、科學計算解決經典物理問題，如力學、熱學、電磁學中的一些常見問題。本文的思路主要是，先介紹經典物理習題，然后對習題進行分析，解答，再通過MATLAB軟件進行編程，模擬實驗結果。通過多次驗證。得到所需答案。再通過圖形繪制，形象的描繪出圖形，與預期結果進行比較、驗證。作出總結。本文展示的MATLAB軟件在解決物理問題中的應用。關鍵詞：力學；熱學；電磁學；MATLAB程序ABSTRACT.MathWorks MATLAB is introduced in a scientific computing software, MATLAB means Matrix Laboratory . MATLAB h

3、as a low starting point, powerful, easy to use, and both numerical calculation and symbolic operation advantages. Using MATLAB, the drawing is very convenient, both to draw various graphics, including the two-dimensional graphics and three-dimensional graphics, graphics can also be modified and cont

4、rolled. This article written by the MATLAB environment to solve by classical physics scientific computing problems, such as mechanical, thermal, electromagnetics some common problems. The main idea of this paper is to introduce classical physics problems, and then exercises to analyze, answer, and t

5、hen programmed by MATLAB software to simulate the experimental results. Through multiple authentication. Get the answers you need. And through graphics rendering, the image depicts the graphics, compared with the expected results to verify. Conclusion. This article presents the MATLAB software to so

6、lve the problem of physics.Key Words：Mechanics;heat;electromagnetism,;MATLAB 目 錄引言11 力學問題31.1質點運動學31.1.1已知質點的運動方程求其速度和加速度31.1.2已知質點的運動方程求質點的軌跡41.1.3考慮空氣阻力的拋射體運動51.1.4已知加速度求速度、運動方程和軌跡71.2盧瑟福散射（Rutherford scattering)研究92 熱學問題122.1理想氣體物態方程122.2范德瓦耳斯方程132.2.1范德瓦耳斯氣體等溫線132.2.2臨界參數153電磁學問題163.1求電偶極子在其所在平面

7、產生的電場中任一點P的電位163.2由電位的表示式計算電場并畫出等電位線和電場方向173.3帶電粒子在電磁場中的運動19結論21參考文獻22致謝2321天津職業技術師范大學2011屆本科生畢業論文引言近幾十年來，計算機技術的廣泛應用已經深入地影響到社會的各個方面，大大加快了社會的變革進程，計算機的應用離不開計算語言，而計算語言本身也處于不斷的發展之中。MATLAB是MATrix LABoratory （矩陣實驗室）的縮寫，它自從1984年由美國MathWorks 公司推出以來，經過不斷改進和發展，現已經成為國際公認的優秀的工程應用開發環境。MATLAB是一種廣泛應用于工程計算及數值分析領域的新

8、型高級語言。它以矩陣作為數據操作的基本單位，使得矩陣運算變得非常簡捷、方便、高效。MATLAB提供了十分豐富的數值計算函數，而且MATLAB和著名的符號計算語言Maple相結合，使得MATLAB具有符號計算功能。MATLAB的繪圖功能也很強，它既可以繪制各種二維、三維圖形，還可以對圖形進行修飾和控制，以增強圖形的表現效果。MATLAB具有編程語言的基本特征，使用MATLAB也可以像使用BASIC、FORTRAN、C等傳統編程語言一樣，進行程序設計，而且簡單易學、編程效率高。MATLAB包含基本部分和各種可選的工具箱，其基本部分構成了MATLAB的核心內容，而MATLAB工具箱擴充了其功能。應用

9、范圍也越來越廣。物理模型的建立及其數學處理在物理學的教學中占有重要地位，而MATLAB在這方面具有獨特的優勢。因此，利用MATLAB這一先進的科學計算語言來輔助物理學的教學工作必將大大提高教學效率。另外，MATLAB起點低、功能強、易學易用以及兼有數值運算和符號運算功能的優點，可以初步掌握這門科學計算語言，并在整個物理學習過程中不斷反復使用是完全必要和可行的。運動學的任務是描述隨時間的推移物體空間位置的變動，不涉及物體間相互作用與運動的關系。本文在力學中主要討論如何使用MATLAB描述質點理想模型的運動，最后引入伽利略變換，它和物理學一條基本原理即相對性原理密切相關。質點平面運動指質點在平面上

10、的曲線運動。這時，質點經常改變運動方向，速度、加速度等物理量的矢量性更突出。如何選擇坐標系的問題更加重要。本文在質點運動方面，主要討論拋體運動，在理想情況下，受空氣阻力、斜拋等得運動軌跡如何在MATLAB中體現出來。以及，已知速度、如何求加速度等。本文在熱學方面主要處理了理想氣體物態方程、范德瓦耳斯方程如何用MATLAB描述出來。理想氣體，只要在足夠寬廣的溫度、壓強變化范圍內進行比較精細的研究，就可發現，氣體的物態方程相當復雜，而且不同氣體所遵循的規律也有所不同。但在壓強趨于零，其溫度不太高也不太低的情況下，不同種類氣體在物態方程上的差異可趨于消失，氣體所遵從的規律也趨于簡單。這種壓強趨于零的

11、極限狀態下的氣體稱為理想氣體。荷蘭物理學家范德瓦耳斯在克勞修斯的論文的啟發下，對理想氣體的兩條基本假定即忽略分子固有體積、忽略除碰撞外分子間相互作用力作出了兩條重要修正，得出了能描述真實氣體行為的范德瓦耳斯方程。在發現電現象2000多年之后，人們才開始對電現象進行定量的研究。1785年，庫倫通過扭秤實驗總結出兩個靜止電荷之間電相互作用的定量規律，通常稱之為庫侖定律。實驗表明，靜電力具有疊加性。原則上，庫侖定律加上靜電力的疊加原理可以求解任意帶電體之間的靜電力。實驗也指出，試探電荷在場中所受的靜電力與試探電荷電量之比反映了電場本身的性質，該比值被稱為電場強度。電場強度也具有疊加性，由場強的定義加

12、上場的疊加原理可以求解任意帶電體的場強分布。本文在電磁學中，主要研究如何用MATLAB求解電偶極子，帶電粒子在電場中運動的問題。本文在物理題目的選取上，主要是普遍、常見的問題，意在將計算語言和物理課程的學習結合起來，起到相輔相成的作用。在程序的編寫中，也力求簡潔。1 力學問題1.1質點運動學在一些問題中，若物體的形狀和大小可以忽略，則可以把該物體視為具有一定質量的幾何點，這就是所謂的質點。質點運動學的基本問題是；已知質點的運動學方程求質點的軌跡、速度和加速度；已知質點的速度或加速度求質點的運動方程和軌跡。下面，結合大家熟悉的幾個具體例子來說明如何用MATLAB處理上述問題1.1.1已知質點的運

13、動方程求其速度和加速度例：某質點的運動學方程為（單位：m，s），求t=1s時質點的速度矢量。解題分析質點的位置矢量為，質點的速度矢量為，質點速度大小和方向余弦分別為程序syms tr=-10,15*t,5*t3; %用數組表示位置矢量V=diff(r,t); %求速度v=sqrt(sum(V.2) %求速度矢量長度，即速度矢量的大小alpha=acos(V(1)/v); beta=acos(V(2)/v); gamma=acos(V(3)/v);%求速度矢量的方向角v1=subs(v,t,1), alpha=subs(alpha,t,1),beta=subs(beta,t,1),gamma=s

14、ubs(gamma,t,1)%求t=1s時質點的速率和方向角，使用了置換命令的函數subs運行結果：v1=21.2132alpha=1.5708beta=0.7854gamma=0.78541.1.2已知質點的運動方程求質點的軌跡例：設一物體以拋射角，速度拋出，落點與射點在同一水平面，且不計空氣阻力。求物體在空氣中飛行的時間、落點距離和飛行的最大高度。解題分析：質點運動學，有 解出t，它就是落點時間.有兩個解，只取其中的一個有效解，然后求最大飛行距離。MATLAB程序：clear all y0=0;x0=5; %取初始位置，為了畫出豎拋運動，未將x0取在原點。v0=input('v0=

15、');theta=input('theta='); %輸入拋射速率和岀射角度v0x=v0*cosd(theta);v0y=v0*sind(theta); %輸入初速度的x分量和y分量ay=-9.81;ax=0; %加速度的y分量和x分量tf=roots(ay/2,v0y,y0); %解出方程的根，求飛行時間。有兩個解，只取有效解tf=max(tf); %落點時間t=0:0.1:tf; %為了畫圖，取時間數組y=y0+v0y*t+ay*t.2/2;x=x0+v0x*t+ax*t.2/2; %t時刻，質點的位置xf=max(x), %飛行達到的最遠距離，即射程yf=max(

16、y), %飛行中達到的最大高度grid on, hold on %畫網格，保持圖形plot(x,y), %畫圖，xlabel('x'),ylabel('y') %坐標標注hold off仿真結果與討論：運行該程序，例如，取初速度=30，岀射角分別取35,45,55,65，75,85，90，則可畫出圖1.1所示曲線，并在命令窗口中給出相應的射程、飛行時間和最大高度。圖 1.1拋體的運動軌跡在上述程序中，我們設定了目標高度為零。我們還可對上述程序進行修改，使其能預先設定目標高度。1.1.3考慮空氣阻力的拋射體運動例：一物體在有阻力的空氣中作拋體運動，設拋體質量為m，

17、初速度為（可設），受到的空氣阻力大小與速率v的一次方成正比，b是比例系數。求拋體的運動軌跡和速度的x、y分量以及速率v隨時間的變化。解題分析：以地面為參考系，以拋出點為原點建立直角坐標系。質點受重力和空氣阻力的作用，其運動微分方程為令y(1)=x, y(2)=dx/dt, y(3)=y, y(4)=dy/dt, 將方程寫成一階微分方程組的形式再用命令函數ode45解此常微分方程組。MATLAB程序m=1; %為簡單起見，取m=1.b=input('b='); %輸入b值，例如，b=0.3.t,y=ode45('ex1',0:0.01:5,0,5,0,19,b,m

18、);%使用了數值法解微分方程的命令函數ode45v=sqrt(y(:,2).2+y(:,4).2);subplot(2,1,1) %繪制子圖plot(y(:,1),y(:,3) %繪制運動軌跡，即x-y曲線，注意：y(1)=x,y(3)=y.subplot(2,1,2)plot(t,y(:,2),t,y(:,4),t,v) %繪制速度的x，y分量以及速率v時間t的變化曲線。函數文件function ydot=ex1(t,y,flag,b,m)ydot=y(2); -b./m.*y(2).*sqrt(y(2).2+y(4).2); y(4);-9.8-b./m.*y(4).*sqrt(y(2).

19、2+y(4).2);運行結果如圖1.2所示?？梢愿淖僢值（例如b分別取0.3和0）來觀察運動軌跡和速度的x分量、y分量 及速率v隨時間的變化。（a)b=0.3(b)b=0圖1.2有空氣阻力時拋射體的運動軌跡及速度隨時間的變化1.1.4已知加速度求速度、運動方程和軌跡例：一質點平面運動的加速度為。初始條件為。求質點軌跡。解題分析：將加速度對時間求積分可得速度，將速度對時間求積分可得位置坐標。在得到參數方程后，給定時間數組，就可以畫出運動軌跡。MATLAB程序clear allsyms t A B v0x v0y x0 y0 vx vy ax ay t t1 t2;v0x=0;v0y=B; x0=

20、A; y0=0; %初始條件ax=-A*cos(t); ay=-B*sin(t); %加速度的x分量和y分量。vx=v0x+int(ax,t,0,t1),vy=v0y+int(ay,t,0,t1), %速度的x分量和y分量x=A+int(vx,t1,0,t2); y=int(vy,t1,0,t2); %求參數方程x=vpa(subs(x,A,B,1,0.5),3)%使用了vpa計算數值；使用subs用數據替換A和B。y=vpa(subs(y,A,B,1,0.5),3)運行結果：vx=-A*sin(t1)vy=-B*cos(t1)x=cos(t2)y=500*sin(t2)%下面 繪制質點的軌跡

21、 cleart2=0:0.1:2*pi;x=cos(t2);y=.500*sin(t2);plot(x,y),xlabel('x'),ylabel('y')從圖1.3給出了程序的運行結果，可知質點的運動軌跡是橢圓。圖1.3：質點的軌跡1.2盧瑟福散射（Rutherford scattering)研究例：盧瑟福等人發現用粒子轟擊金箔時有些入射粒子散射偏轉角很大，甚至超過。盧瑟福于1911年提出原子必有以帶正電的核心，即原子核；此即原子結構的行星模型。已知粒子的質量為，以速度接近電荷為Ze的重原子核，瞄準距離為b，如圖所示。設原子核質量比粒子大很多，可以近似看作靜止

22、。1 求粒子接近重原子核最近距離d。2 畫出粒子在不同初始條件下的軌道，并通過改變初始條件來研究影響散射角的因素。解題分析1 粒子受靜電力始終指向重核中心，粒子在一平面內運動。設z軸垂直于此平面且通過重核中心，則粒子所受靜電力對z軸的力矩為零，即對z軸的角動量守恒。粒子以速度運動，對z軸的角動量是，粒子最接近重核（距離為d）時，速度應與粒子至核的連線垂直，角動量為。于是或（1）在散射過程中，只有庫侖斥力作用，故能量守恒。（2）其中，左邊第二項是庫侖斥力勢能。聯解（1）、（2）式，可得d的表達式。2 選擇在直角坐標系，原點位于力心重核處。根據牛頓運動定律，粒子的運動方程在直角坐標中的投影方程為令

23、，則上述方程組可寫為令粒子沿Ox方向入射，入射速率為，初始條件為。為了能得到多粒子的運動軌跡，程序中采用input函數給出不同初始條件。MATLAB程序1 求粒子接近重原子核最近距離d。syms v v0 b k Z e m d;d,v=solve('v=v0*b/d','m*v2/2+k*Z*e2/d=m*v02/2',d,v)運行結果：d=1/2/m/v02*(2*k*Z*e2+2*(k2*Z2*e4+m2*v04*b2)(1/2)v=2*v03*b*m/(2*k*Z*e2+2*(k2*Z2*e4+m2*v04*b2)(1/2)即2 畫出粒子在不同初始條件下

24、的軌道y0=input('請輸入初始條件:'); %例如，可輸入：-20 1 10 0;line(0,0,'marker','.','markersize',50,'color','r');text(2,0,'靶粒子');hold ont,y=ode23('ex2f',0:.1:42,y0,3);axis(-20 20,-20 20)plot(y(:,1),y(:,3),hold on以下是獨立的函數文件，文件名為ex2f.m,其中。function ydot=ex2

25、f(t,y,flag,p)ydot=y(2);p*y(1)./sqrt(y(1).*y(1)+y(3).*y(3).3;y(4);p*y(3)./sqrt(y(1).*y(1)+y(3).*y(3).3;運行結果如圖1.4所示。圖1.4：粒子的散射軌道2 熱學問題2.1理想氣體物態方程理想氣體是將實際氣體外推到壓強趨于零的極限情況下得到的一個理想模型。1857年，克勞修斯進一步提出了理想氣體的微觀模型，并通過計算氣體的壓強得到了理想氣體的物態方程。而在此之前，理想氣體物態方程是由氣體三大實驗定律外推得到的。例：編寫一個繪制帶有等高線的理想氣體狀態方程pV=RT的曲面。解題分析理想氣體的物態方程

26、為。其中，p，V，分別為氣體的壓強、體積和摩爾數，R為氣體普適常數，其值為R=8.31J/mol*K.圖2.1：理想氣體狀態方程曲面圖MATLAB程序clearR=8.31;p=(1:20).*1e5;v=(1:20)*1e-3;v,p=meshgrid(v,p);T=p.*v./R;meshc(v,p,T),xlabel('v'),ylabel('p'),zlabel('T'),運行結果如圖2.1所示。實驗指出，理想氣體狀態方程在一定程度上反映了真實氣體的性質，但對低溫和高密度狀態下的氣體以及氣體和液體之間的相變卻無能為力，因而是一個理想的“永

27、久氣體”狀態方程。2.2范德瓦耳斯方程范德瓦耳斯方程1873年，荷蘭物理學家范德瓦耳斯（van der Waals)在克勞修斯的理想氣體模型和安德魯斯發現的臨界點現象的啟發下，考慮了分子體積和分子間吸引力這兩個因素，對理想氣體進行了修正，得到了能描述真實氣體行為的范德瓦耳斯方程：其中，常數a和b分別是1mol范氏氣體的壓強修正系數和體積修正系數，其數值隨氣體種類的不同而異。下表1列出了幾種氣體的a，b值及臨界參量。 表1：幾種氣體的a、b值及臨界參量氣體 aTc(K)Ar0.13490.0332747.9789150.86760.09680.36060.0428072.9079304.0640

28、0.1284CO0.14560.0395434.4925132.89470.1186He0.0034140.023712.24925.19650.07110.024610.0266812.804933.28960.0800Ne0.020910.0169726.892244.4680.05090.13500.0386433.4885126.08970.11590.13630.0318449.7950154.49200.09552.2.1范德瓦耳斯氣體等溫線例:編寫一個繪制范德瓦爾斯氣體等溫線的程序，要求輸入溫度值后便可畫出相應的等溫線。解題分析以二氧化碳為例，從表1查得，由范德瓦爾斯方程可繪制等

29、溫線簇，溫度選取如圖所示。MATLAB程序v=(0.06:0.001:1).*1e-3;T=input('T=');b=0.0428e-3;a=0.3606;R=8.31;p=R.*T./(v-b)-a./v.2;grid on,plot(v,p),axis(0,0.4e-3,-2e7,6e7)hold on;運行結果如圖2.2所示圖2.2：范德瓦耳斯氣體等溫線范德瓦爾斯方程不僅對氣體性質的描述優于理想氣體物態方程，而且還能描述液相及氣、液兩相轉變的性質以及臨界點的特征。2.2.2臨界參數范德瓦爾斯等溫線中有一個特殊的狀態臨界點。臨界點所對應的壓強、體積和溫度分別稱為臨界壓強、

30、臨界體積和臨界溫度。在臨界點所發生的氣液相變與在低于臨界溫度時的相變完全相同，屬于二級相變；而低于臨界點是的氣液相變屬于一級相變。在臨界點以上，氣體是不能夠通過等溫壓縮被轉變為液相的。系統在臨界點具有許多特殊性質，稱為臨界現象。下面來介紹臨界點的確定。例:由范德瓦爾斯物態方程求臨界參量、。解題分析從圖可以看出，臨界點是一拐點，它同時滿足下列條件：，利用上述拐點條件，將范德瓦爾斯方程對求導并聯解方程，便可求得范德瓦爾斯氣體的三個臨界參量。MATLAB程序clearsyms a b R TD1=diff('(p+a/v2)*(v-b)-R*T','v');D2=di

31、ff(D1,'v');pc,vc=solve(D1,D2,'v','p')Tc=solve(pc+a/vc2)*(vc-b)-R*T,'T')運行結果pc=1/27*a/b2vc=3*bTc=8/27*a/b/R即；，3電磁學問題3.1求電偶極子在其所在平面產生的電場中任一點P的電位例：已知電偶極子中兩電荷-q和+q的距離為。計算中可取。解題分析設場點P到的距離為和，則單獨存在時P點的電位分別為由電位疊加原理，電偶極子產生的電場在P點的電位為MATLAB程序clear;q=1.6e-19; %單位電荷電量C0=9e9;l=3.0;

32、 %偶極子正負電荷之間的距離lx=-5:0.5:5;y=x;X,Y=meshgrid(x,y);r1=sqrt(X-1./2).2+(Y-0).2); %電荷距離空間P（x，y）點的距離r2=sqrt(X+1./2).2+(Y-0).2);U=q.*C0.*(1./r1-1./r2); %求出空間任意一點P（x，y）的電位plot(-1/2,0,'ro',-1/2,0,'r-') %標出負電荷hold on,plot(1/2,0,'ro',1/2,0,'r+') %標出正電荷C=contour(X,Y,U,'k-'

33、);clabel(C); %畫等位線并標出電位值axis('square')運行結果如圖3.1所示 3.1電位梯度與電場強度電位是標量，它在空間中每點都有一定的數值，所以電位是標量場。標量場在空間中沿不同方向的變化率稱為梯度，對電位場而言稱為電位梯度，用grad U或來表示。可以證明，電位梯度和電場強度E的關系為利用上式，可從已知的電位分布求電場強度。3.2由電位的表示式計算電場并畫出等電位線和電場方向解題分析如果已知空間的電位分布則空間的電場強度為按照本題的要求，可利用讀入字符串的指令input('U'(x,y)=','s')來輸入電位

34、方程。在MATLAB中，梯度函數的調用格式為gradient（），它是靠數值微分得到的。因此，空間觀測點應取得密一些，以獲得較高的精度。MATLAB程序clear allU=input('請輸入電位方程，U=(x,y)=','s'); %例如，取U(x,y)=log(x.2+y.2)。xmax=5;ymax=5;Ngrid=20; %繪圖區從x=-xmax到xmax，網格線數為20xplot=linspace(-xmax,xmax,Ngrid); %繪圖用x的數組x,y=meshgrid(xplot); %x，y取同樣范圍，生成二維網格Uplot=eval(U)

35、; %執行輸入的字符串U,計算各點U的值Explot,Eyplot=gradient(-Uplot); %電場等于電位的負梯度clf;subplot(1,2,1),meshc(Uplot); %劃分子圖；繪制含等位線的三維曲面xlabel('x');ylabel('y');zlabel('U');subplot(1,2,2),axis(-xmax,xmax,-ymax,ymax); %規定等位線的范圍cs=contour(x,y,Uplot); %畫等位線，cs是等位線值clabel(cs); %標出等位線的值hold on;quiver(x,y

36、,Explot,Eyplot); %保持圖形，在原圖形上疊加矢量場圖xlabel('x');ylabel('y');hold off;圖3.2：的電位分布與電場分布運行上述程序，所得結果如圖3.2所示。3.3帶電粒子在電磁場中的運動例：設質量為m，帶電量為q的粒子在電磁場中的運動微分方程為選場中某點為原點，以為方向，沿方向，建立坐標系。令，上式的投影方程為令，上述方程可改寫為下列一階微分方程組：MATLAB程序%符號法求離子運動微分方程的特解并繪圖clearsyms w x y z t B E m q;E=input('E=');B=input(

37、'B='); %輸入E和B值x,y,z=dsolve('D2x=q*B/m*Dy','D2y=q*E/m-q*B/m*Dx','D2z=0','x(0)=0','y(0)=0','z(0)=0','Dx(0)=0.01','Dy(0)=6','Dz(0)=0.01');%初始條件取x(0)=y(0)=z(0)=0,Dx(0)=0.01,Dy(0)=0.01q=1.6e-2;m=0.02X=subs(x,y,z);x=X(1),y=X(2

38、),z=X(3),ezplot3(X(1),X(2),X(3)運行上述程序，例如，取E=4，B=8可得下列特解并給出圖3.3x =-15/16*cos(32/5*t)-49/640*sin(32/5*t)+1/2*t+15/16y =15/16*sin(32/5*t)-49/640*cos(32/5*t)+49/640z =1/100*t(a)E=4,B=8(b)E=0.01,B=8(c)E=8,B=1圖3.3：帶電粒子在電磁場中的運動下面我們給出一段用數值方法求解該問題的程序，用以比較。q=1.6e-2;m=0.02;B=2;2;0;E=1;0;1;figurestrd1='Eneq 0,Bneq 0'strd2='E=0,Bneq 0'strd3='Eneq 0,B=0'for i=1:3t,y=ode23('ex3f',0:0.1:20,0,0.01,0,6,0,0.01,q,m,B(i),E(i);axes('unit',&#