畢業(yè)論文應(yīng)用matlab求解經(jīng)典物理若干典型問題

1、<p>　　天津****大學(xué)本科生畢業(yè)論文</p><p>　　應(yīng)用MATLAB求解經(jīng)典物理若干典型問題</p><p>　　The application of MATLAB in solving some classical physics questions</p><p><b>　　專業(yè)班級：物理 </b></p>

2、;<p>　　學(xué)生姓名： </p><p><b>　　指導(dǎo)教師： </b></p><p>　　學(xué) 院：理學(xué)院</p><p><b>　　2011年5月</b></p><p><b>　　摘 要</b></p><p&

3、gt;　　MATLAB是 MathWorks公司推出的一套科學(xué)計算軟件，MATLAB的意思是矩陣實驗室。MATLAB具有起點低、功能強(qiáng)大、易學(xué)易用以及兼有數(shù)值運(yùn)算和符號運(yùn)算功能的優(yōu)點。利用MATLAB，繪圖十分方便，它既可以繪制各種圖形，包括二維圖形和三維圖形，還可以對圖形進(jìn)行修飾和控制。本文通過在MATLAB環(huán)境下編寫通過科學(xué)計算解決經(jīng)典物理問題，如力學(xué)、熱學(xué)、電磁學(xué)中的一些常見問題。本文的思路主要是，先介紹經(jīng)典物理習(xí)題，然后對習(xí)題進(jìn)

4、行分析，解答，再通過MATLAB軟件進(jìn)行編程，模擬實驗結(jié)果。通過多次驗證。得到所需答案。再通過圖形繪制，形象的描繪出圖形，與預(yù)期結(jié)果進(jìn)行比較、驗證。作出總結(jié)。本文展示的MATLAB軟件在解決物理問題中的應(yīng)用。</p><p>　　關(guān)鍵詞：力學(xué)；熱學(xué)；電磁學(xué)；MATLAB程序</p><p><b>　　ABSTRACT</b></p><p>

5、　　.MathWorks MATLAB is introduced in a scientific computing software, MATLAB means Matrix Laboratory . MATLAB has a low starting point, powerful, easy to use, and both numerical calculation and symbolic operation advanta

6、ges. 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 controlled. This article written

7、by the MATLAB environment to solve by classical physics scient</p><p>　　Key Words：Mechanics;heat;electromagnetism,;MATLAB </p><p><b>　　目 錄</b></p><p><b>　　引言1&

8、lt;/b></p><p><b>　　1 力學(xué)問題3</b></p><p>　　1.1質(zhì)點運(yùn)動學(xué)3</p><p>　　1.1.1已知質(zhì)點的運(yùn)動方程求其速度和加速度3</p><p>　　1.1.2已知質(zhì)點的運(yùn)動方程求質(zhì)點的軌跡4</p><p>　　1.1.3考慮空氣阻力的拋

9、射體運(yùn)動5</p><p>　　1.1.4已知加速度求速度、運(yùn)動方程和軌跡7</p><p>　　1.2盧瑟福散射（Rutherford scattering)研究8</p><p><b>　　2 熱學(xué)問題11</b></p><p>　　2.1理想氣體物態(tài)方程11</p><p>　

10、　2.2范德瓦耳斯方程12</p><p>　　2.2.1范德瓦耳斯氣體等溫線12</p><p>　　2.2.2臨界參數(shù)14</p><p><b>　　3電磁學(xué)問題15</b></p><p>　　3.1求電偶極子在其所在平面產(chǎn)生的電場中任一點P的電位15</p><p>　　3.2

11、由電位的表示式計算電場并畫出等電位線和電場方向16</p><p>　　3.3帶電粒子在電磁場中的運(yùn)動18</p><p><b>　　結(jié)論20</b></p><p><b>　　參考文獻(xiàn)21</b></p><p><b>　　致謝22</b></p>

12、<p><b>　　引言</b></p><p>　　近幾十年來，計算機(jī)技術(shù)的廣泛應(yīng)用已經(jīng)深入地影響到社會的各個方面，大大加快了社會的變革進(jìn)程，計算機(jī)的應(yīng)用離不開計算語言，而計算語言本身也處于不斷的發(fā)展之中。</p><p>　　MATLAB是MATrix LABoratory （矩陣實驗室）的縮寫，它自從1984年由美國MathWorks 公司推出以來

13、，經(jīng)過不斷改進(jìn)和發(fā)展，現(xiàn)已經(jīng)成為國際公認(rèn)的優(yōu)秀的工程應(yīng)用開發(fā)環(huán)境。</p><p>　　MATLAB是一種廣泛應(yīng)用于工程計算及數(shù)值分析領(lǐng)域的新型高級語言。它以矩陣作為數(shù)據(jù)操作的基本單位，使得矩陣運(yùn)算變得非常簡捷、方便、高效。MATLAB提供了十分豐富的數(shù)值計算函數(shù)，而且MATLAB和著名的符號計算語言Maple相結(jié)合，使得MATLAB具有符號計算功能。MATLAB的繪圖功能也很強(qiáng)，它既可以繪制各種二維、三維圖形，

14、還可以對圖形進(jìn)行修飾和控制，以增強(qiáng)圖形的表現(xiàn)效果。MATLAB具有編程語言的基本特征，使用MATLAB也可以像使用BASIC、FORTRAN、C等傳統(tǒng)編程語言一樣，進(jìn)行程序設(shè)計，而且簡單易學(xué)、編程效率高。MATLAB包含基本部分和各種可選的工具箱，其基本部分構(gòu)成了MATLAB的核心內(nèi)容，而MATLAB工具箱擴(kuò)充了其功能。應(yīng)用范圍也越來越廣。</p><p>　　物理模型的建立及其數(shù)學(xué)處理在物理學(xué)的教學(xué)中占有重要地

15、位，而MATLAB在這方面具有獨特的優(yōu)勢。因此，利用MATLAB這一先進(jìn)的科學(xué)計算語言來輔助物理學(xué)的教學(xué)工作必將大大提高教學(xué)效率。另外，MATLAB起點低、功能強(qiáng)、易學(xué)易用以及兼有數(shù)值運(yùn)算和符號運(yùn)算功能的優(yōu)點，可以初步掌握這門科學(xué)計算語言，并在整個物理學(xué)習(xí)過程中不斷反復(fù)使用是完全必要和可行的。</p><p>　　運(yùn)動學(xué)的任務(wù)是描述隨時間的推移物體空間位置的變動，不涉及物體間相互作用與運(yùn)動的關(guān)系。本文在力學(xué)中主要

16、討論如何使用MATLAB描述質(zhì)點理想模型的運(yùn)動，最后引入伽利略變換，它和物理學(xué)一條基本原理即相對性原理密切相關(guān)。質(zhì)點平面運(yùn)動指質(zhì)點在平面上的曲線運(yùn)動。這時，質(zhì)點經(jīng)常改變運(yùn)動方向，速度、加速度等物理量的矢量性更突出。如何選擇坐標(biāo)系的問題更加重要。本文在質(zhì)點運(yùn)動方面，主要討論拋體運(yùn)動，在理想情況下，受空氣阻力、斜拋等得運(yùn)動軌跡如何在MATLAB中體現(xiàn)出來。以及，已知速度、如何求加速度等。</p><p>　　本文在熱

17、學(xué)方面主要處理了理想氣體物態(tài)方程、范德瓦耳斯方程如何用MATLAB描述出來。理想氣體，只要在足夠?qū)拸V的溫度、壓強(qiáng)變化范圍內(nèi)進(jìn)行比較精細(xì)的研究，就可發(fā)現(xiàn)，氣體的物態(tài)方程相當(dāng)復(fù)雜，而且不同氣體所遵循的規(guī)律也有所不同。但在壓強(qiáng)趨于零，其溫度不太高也不太低的情況下，不同種類氣體在物態(tài)方程上的差異可趨于消失，氣體所遵從的規(guī)律也趨于簡單。這種壓強(qiáng)趨于零的極限狀態(tài)下的氣體稱為理想氣體。荷蘭物理學(xué)家范德瓦耳斯在克勞修斯的論文的啟發(fā)下，對理想氣體的兩條基

18、本假定即忽略分子固有體積、忽略除碰撞外分子間相互作用力作出了兩條重要修正，得出了能描述真實氣體行為的范德瓦耳斯方程。</p><p>　　在發(fā)現(xiàn)電現(xiàn)象2000多年之后，人們才開始對電現(xiàn)象進(jìn)行定量的研究。1785年，庫倫通過扭秤實驗總結(jié)出兩個靜止電荷之間電相互作用的定量規(guī)律，通常稱之為庫侖定律。實驗表明，靜電力具有疊加性。原則上，庫侖定律加上靜電力的疊加原理可以求解任意帶電體之間的靜電力。實驗也指出，試探電荷在場中

19、所受的靜電力與試探電荷電量之比反映了電場本身的性質(zhì)，該比值被稱為電場強(qiáng)度。電場強(qiáng)度也具有疊加性，由場強(qiáng)的定義加上場的疊加原理可以求解任意帶電體的場強(qiáng)分布。本文在電磁學(xué)中，主要研究如何用MATLAB求解電偶極子，帶電粒子在電場中運(yùn)動的問題。</p><p>　　本文在物理題目的選取上，主要是普遍、常見的問題，意在將計算語言和物理課程的學(xué)習(xí)結(jié)合起來，起到相輔相成的作用。在程序的編寫中，也力求簡潔。</p>

20、<p><b>　　1 力學(xué)問題</b></p><p><b>　　1.1質(zhì)點運(yùn)動學(xué)</b></p><p>　　在一些問題中，若物體的形狀和大小可以忽略，則可以把該物體視為具有一定質(zhì)量的幾何點，這就是所謂的質(zhì)點。質(zhì)點運(yùn)動學(xué)的基本問題是；已知質(zhì)點的運(yùn)動學(xué)方程求質(zhì)點的軌跡、速度和加速度；已知質(zhì)點的速度或加速度求質(zhì)點的運(yùn)動方程和軌跡。

21、下面，結(jié)合大家熟悉的幾個具體例子來說明如何用MATLAB處理上述問題</p><p>　　1.1.1已知質(zhì)點的運(yùn)動方程求其速度和加速度</p><p>　　例：某質(zhì)點的運(yùn)動學(xué)方程為（單位：m，s），求t=1s時質(zhì)點的速度矢量。</p><p><b>　　解題分析</b></p><p>　　質(zhì)點的位置矢量為，<

22、/p><p>　　質(zhì)點的速度矢量為，</p><p>　　質(zhì)點速度大小和方向余弦分別為</p><p><b>　　程序</b></p><p><b>　　syms t</b></p><p>　　r=[-10,15*t,5*t^3]; %用數(shù)組表示位置矢量</

23、p><p>　　V=diff(r,t); %求速度</p><p>　　v=sqrt(sum(V.^2)) %求速度矢量長度，即速度矢量的大小</p><p>　　alpha=acos(V(1)/v); beta=acos(V(2)/v); gamma=acos(V(3)/v);</p><p>　　%求速度矢

24、量的方向角</p><p>　　v1=subs(v,t,1), alpha=subs(alpha,t,1),beta=subs(beta,t,1),</p><p>　　gamma=subs(gamma,t,1)</p><p>　　%求t=1s時質(zhì)點的速率和方向角，使用了置換命令的函數(shù)subs</p><p><b>　　運(yùn)行結(jié)果

25、：</b></p><p>　　v1=21.2132</p><p>　　alpha=1.5708</p><p>　　beta=0.7854</p><p>　　gamma=0.7854</p><p>　　1.1.2已知質(zhì)點的運(yùn)動方程求質(zhì)點的軌跡</p><p>　　例：設(shè)一物體

26、以拋射角，速度拋出，落點與射點在同一水平面，且不計空氣阻力。求物體在空氣中飛行的時間、落點距離和飛行的最大高度。</p><p>　　解題分析：質(zhì)點運(yùn)動學(xué)，有</p><p>　　解出t，它就是落點時間.有兩個解，只取其中的一個有效解，然后求最大飛行距離。</p><p><b>　　MATLAB程序：</b></p><p

27、>　　clear all </p><p>　　y0=0;x0=5; %取初始位置，為了畫出豎拋運(yùn)動，未將x0取在原點。</p><p>　　v0=input('v0=');theta=input('theta='); %輸入拋射速率和岀射角度</p><p&g

28、t;　　v0x=v0*cosd(theta);</p><p>　　v0y=v0*sind(theta); %輸入初速度的x分量和y分量</p><p>　　ay=-9.81;ax=0; %加速度的y分量和x分量</p><p>　　tf=roots([ay/2,v0y,y0]); %解出方程的根，求飛行

29、時間。有兩個解，只取有效解</p><p>　　tf=max(tf); %落點時間</p><p>　　t=0:0.1:tf; %為了畫圖，取時間數(shù)組</p><p>　　y=y0+v0y*t+ay*t.^2/2;x=x0+v0x*t+ax*t.^2/2; %t時刻，質(zhì)點的位置</p&g

30、t;<p>　　xf=max(x), %飛行達(dá)到的最遠(yuǎn)距離，即射程</p><p>　　yf=max(y), %飛行中達(dá)到的最大高度</p><p>　　grid on, hold on %畫網(wǎng)格，保持圖形</p><p>　　plot(x,y),

31、 %畫圖，</p><p>　　xlabel('x'),ylabel('y') %坐標(biāo)標(biāo)注</p><p><b>　　hold off</b></p><p><b>　　仿真結(jié)果與討論：</b></p><p>　　運(yùn)行該程

32、序，例如，取初速度=30，岀射角分別取35,45,55,65，75,85，90，則可畫出圖1.1所示曲線，并在命令窗口中給出相應(yīng)的射程、飛行時間和最大高度。</p><p>　　圖 1.1拋體的運(yùn)動軌跡</p><p>　　在上述程序中，我們設(shè)定了目標(biāo)高度為零。我們還可對上述程序進(jìn)行修改，使其能預(yù)先設(shè)定目標(biāo)高度。</p><p>　　1.1.3考慮空氣阻力的拋射體運(yùn)

33、動</p><p>　　例：一物體在有阻力的空氣中作拋體運(yùn)動，設(shè)拋體質(zhì)量為m，初速度為（可設(shè)），受到的空氣阻力大小與速率v的一次方成正比，b是比例系數(shù)。求拋體的運(yùn)動軌跡和速度的x、y分量以及速率v隨時間的變化。</p><p><b>　　解題分析：</b></p><p>　　以地面為參考系，以拋出點為原點建立直角坐標(biāo)系。質(zhì)點受重力和空氣阻力

34、的作用，其運(yùn)動微分方程為</p><p>　　令y(1)=x, y(2)=dx/dt, y(3)=y, y(4)=dy/dt, 將方程寫成一階微分方程組的形式</p><p>　　再用命令函數(shù)ode45解此常微分方程組。</p><p><b>　　MATLAB程序</b></p><p>　　m=1;

35、 %為簡單起見，取m=1.</p><p>　　b=input('b='); %輸入b值，例如，b=0.3.</p><p>　　[t,y]=ode45('ex1',[0:0.01:5],[0,5,0,19],[],b,m);</p><p>　　%使用了數(shù)值法解微分方程的命令函數(shù)ode45</p><

36、p>　　v=sqrt(y(:,2).^2+y(:,4).^2);</p><p>　　subplot(2,1,1) %繪制子圖</p><p>　　plot(y(:,1),y(:,3)) %繪制運(yùn)動軌跡，即x-y曲線，注意：y(1)=x,y(3)=y.</p><p>　　subplot(2,1,2)</p><p

37、>　　plot(t,y(:,2),t,y(:,4),t,v) %繪制速度的x，y分量以及速率v時間t的變化曲線。</p><p><b>　　函數(shù)文件</b></p><p>　　function ydot=ex1(t,y,flag,b,m)</p><p>　　ydot=[y(2);</p><p>　　-b./

38、m.*y(2).*sqrt(y(2).^2+y(4).^2);</p><p><b>　　y(4);</b></p><p>　　-9.8-b./m.*y(4).*sqrt(y(2).^2+y(4).^2)];</p><p>　　運(yùn)行結(jié)果如圖1.2所示?？梢愿淖僢值（例如b分別取0.3和0）來觀察運(yùn)動軌跡和速度的x分量、y分量 及速率v隨時

39、間的變化。</p><p><b>　?。╝)b=0.3</b></p><p><b>　　(b)b=0</b></p><p>　　圖1.2有空氣阻力時拋射體的運(yùn)動軌跡及速度隨時間的變化</p><p>　　1.1.4已知加速度求速度、運(yùn)動方程和軌跡</p><p>　　

40、例：一質(zhì)點平面運(yùn)動的加速度為。初始條件為。求質(zhì)點軌跡。</p><p><b>　　解題分析：</b></p><p>　　將加速度對時間求積分可得速度，將速度對時間求積分可得位置坐標(biāo)。在得到參數(shù)方程后，給定時間數(shù)組，就可以畫出運(yùn)動軌跡。</p><p><b>　　MATLAB程序</b></p><

41、p><b>　　clear all</b></p><p>　　syms t A B v0x v0y x0 y0 vx vy ax ay t t1 t2;</p><p>　　v0x=0;v0y=B; x0=A; y0=0; %初始條件</p><p>　　ax=-A*cos(t); ay=-B*sin(t); %加

42、速度的x分量和y分量。</p><p>　　vx=v0x+int(ax,t,0,t1),vy=v0y+int(ay,t,0,t1), %速度的x分量和y分量</p><p>　　x=A+int(vx,t1,0,t2); y=int(vy,t1,0,t2); %求參數(shù)方程</p><p>　　x=vpa(subs(x,{A,B},{1,0.5}),

43、3)</p><p>　　%使用了vpa計算數(shù)值；使用subs用數(shù)據(jù)替換A和B。</p><p>　　y=vpa(subs(y,{A,B},{1,0.5}),3)</p><p><b>　　運(yùn)行結(jié)果：</b></p><p>　　vx=-A*sin(t1)</p><p>　　vy=-B*cos

44、(t1)</p><p><b>　　x=cos(t2)</b></p><p>　　y=500*sin(t2)%</p><p>　　下面 繪制質(zhì)點的軌跡 </p><p><b>　　clear</b></p><p>

45、;　　t2=0:0.1:2*pi;</p><p>　　x=cos(t2);</p><p>　　y=.500*sin(t2);</p><p>　　plot(x,y),xlabel('x'),ylabel('y')</p><p>　　從圖1.3給出了程序的運(yùn)行結(jié)果，可知質(zhì)點的運(yùn)動軌跡是橢圓。</p&g

46、t;<p>　　圖1.3：質(zhì)點的軌跡</p><p>　　1.2盧瑟福散射（Rutherford scattering)研究</p><p>　　例：盧瑟福等人發(fā)現(xiàn)用粒子轟擊金箔時有些入射粒子散射偏轉(zhuǎn)角很大，甚至超過。盧瑟福于1911年提出原子必有以帶正電的核心，即原子核；此即原子結(jié)構(gòu)的行星模型。</p><p>　　已知粒子的質(zhì)量為，以速度接近電荷為

47、Ze的重原子核，瞄準(zhǔn)距離為b，如圖所示。設(shè)原子核質(zhì)量比粒子大很多，可以近似看作靜止。</p><p>　　求粒子接近重原子核最近距離d。</p><p>　　畫出粒子在不同初始條件下的軌道，并通過改變初始條件來研究影響散射角的因素。</p><p><b>　　解題分析</b></p><p>　　粒子受靜電力始終指向重

48、核中心，粒子在一平面內(nèi)運(yùn)動。設(shè)z軸垂直于此平面且通過重核中心，則粒子所受靜電力對z軸的力矩為零，即對z軸的角動量守恒。粒子以速度運(yùn)動，對z軸的角動量是，粒子最接近重核（距離為d）時，速度應(yīng)與粒子至核的連線垂直，角動量為。于是</p><p>　　或（1）</p><p>　　在散射過程中，只有庫侖斥力作用，故能量守恒。</p><p><b

49、>　?。?）</b></p><p>　　其中，左邊第二項是庫侖斥力勢能。聯(lián)解（1）、（2）式，可得d的表達(dá)式。</p><p>　　選擇在直角坐標(biāo)系，原點位于力心重核處。根據(jù)牛頓運(yùn)動定律，粒子的運(yùn)動方程在直角坐標(biāo)中的投影方程為</p><p>　　令，，，則上述方程組可寫為</p><p>　　令粒子沿Ox方向入射，入

50、射速率為，初始條件為。為了能得到多粒子的運(yùn)動軌跡，程序中采用input函數(shù)給出不同初始條件。</p><p><b>　　MATLAB程序</b></p><p>　　求粒子接近重原子核最近距離d。</p><p>　　syms v v0 b k Z e m d;</p><p>　　[d,v]=solve('v

51、=v0*b/d','m*v^2/2+k*Z*e^2/d=m*v0^2/2',d,v)</p><p><b>　　運(yùn)行結(jié)果：</b></p><p>　　d=1/2/m/v0^2*(2*k*Z*e^2+2*(k^2*Z^2*e^4+m^2*v0^4*b^2)^(1/2))</p><p>　　v=2*v0^3*b*m/(

52、2*k*Z*e^2+2*(k^2*Z^2*e^4+m^2*v0^4*b^2)^(1/2))</p><p><b>　　即</b></p><p>　　畫出粒子在不同初始條件下的軌道</p><p>　　y0=input('請輸入初始條件:'); %例如，可輸入：[-20 1 10 0];</p>&l

53、t;p>　　line(0,0,'marker','.','markersize',50,'color','r');</p><p>　　text(2,0,'靶粒子');hold on</p><p>　　[t,y]=ode23('ex2f',[0:.1:42],y0,[],

54、3);</p><p>　　axis([-20 20,-20 20])</p><p>　　plot(y(:,1),y(:,3)),hold on</p><p>　　以下是獨立的函數(shù)文件，文件名為ex2f.m,其中。</p><p>　　function ydot=ex2f(t,y,flag,p)</p><p>　

55、　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];</p><p>　　運(yùn)行結(jié)果如圖1.4所示。</p><p>　　圖1.4：粒子的散射軌道</p><p><b>　　2 熱學(xué)問題</b></p

56、><p>　　2.1理想氣體物態(tài)方程</p><p>　　理想氣體是將實際氣體外推到壓強(qiáng)趨于零的極限情況下得到的一個理想模型。1857年，克勞修斯進(jìn)一步提出了理想氣體的微觀模型，并通過計算氣體的壓強(qiáng)得到了理想氣體的物態(tài)方程。而在此之前，理想氣體物態(tài)方程是由氣體三大實驗定律外推得到的。</p><p>　　例：編寫一個繪制帶有等高線的理想氣體狀態(tài)方程pV=RT的曲面。&l

57、t;/p><p><b>　　解題分析</b></p><p>　　理想氣體的物態(tài)方程為。其中，p，V，分別為氣體的壓強(qiáng)、體積和摩爾數(shù)，R為氣體普適常數(shù)，其值為R=8.31J/mol*K.</p><p>　　圖2.1：理想氣體狀態(tài)方程曲面圖</p><p><b>　　MATLAB程序</b><

58、/p><p><b>　　clear</b></p><p><b>　　R=8.31;</b></p><p>　　p=(1:20).*1e5;</p><p>　　v=(1:20)*1e-3;</p><p>　　[v,p]=meshgrid(v,p);</p>

59、<p>　　T=p.*v./R;</p><p>　　meshc(v,p,T),</p><p>　　xlabel('v'),ylabel('p'),</p><p>　　zlabel('T'),</p><p>　　運(yùn)行結(jié)果如圖2.1所示。</p><p>

60、　　實驗指出，理想氣體狀態(tài)方程在一定程度上反映了真實氣體的性質(zhì)，但對低溫和高密度狀態(tài)下的氣體以及氣體和液體之間的相變卻無能為力，因而是一個理想的“永久氣體”狀態(tài)方程。</p><p>　　2.2范德瓦耳斯方程</p><p><b>　　范德瓦耳斯方程</b></p><p>　　1873年，荷蘭物理學(xué)家范德瓦耳斯（van der Waals)

61、在克勞修斯的理想氣體模型和安德魯斯發(fā)現(xiàn)的臨界點現(xiàn)象的啟發(fā)下，考慮了分子體積和分子間吸引力這兩個因素，對理想氣體進(jìn)行了修正，得到了能描述真實氣體行為的范德瓦耳斯方程：</p><p>　　其中，常數(shù)a和b分別是1mol范氏氣體的壓強(qiáng)修正系數(shù)和體積修正系數(shù)，其數(shù)值隨氣體種類的不同而異。下表1列出了幾種氣體的a，b值及臨界參量。</p><p>　　表1：幾種氣體的a、b值及臨界參量</p

62、><p>　　2.2.1范德瓦耳斯氣體等溫線</p><p>　　例:編寫一個繪制范德瓦爾斯氣體等溫線的程序，要求輸入溫度值后便可畫出相應(yīng)的等溫線。</p><p><b>　　解題分析</b></p><p>　　以二氧化碳為例，從表1查得，，由范德瓦爾斯方程</p><p>　　可繪制等溫線簇，溫

63、度選取如圖所示。</p><p><b>　　MATLAB程序</b></p><p>　　v=(0.06:0.001:1).*1e-3;</p><p>　　T=input('T=');</p><p>　　b=0.0428e-3;</p><p><b>　　a=0.

64、3606;</b></p><p><b>　　R=8.31;</b></p><p>　　p=R.*T./(v-b)-a./v.^2;</p><p>　　grid on,plot(v,p),</p><p>　　axis([0,0.4e-3,-2e7,6e7])</p><p>&

65、lt;b>　　hold on;</b></p><p>　　運(yùn)行結(jié)果如圖2.2所示</p><p>　　圖2.2：范德瓦耳斯氣體等溫線</p><p>　　范德瓦爾斯方程不僅對氣體性質(zhì)的描述優(yōu)于理想氣體物態(tài)方程，而且還能描述液相及氣、液兩相轉(zhuǎn)變的性質(zhì)以及臨界點的特征。</p><p><b>　　2.2.2臨界參數(shù)

66、</b></p><p>　　范德瓦爾斯等溫線中有一個特殊的狀態(tài)——臨界點。臨界點所對應(yīng)的壓強(qiáng)、體積和溫度分別稱為臨界壓強(qiáng)、臨界體積和臨界溫度。在臨界點所發(fā)生的氣液相變與在低于臨界溫度時的相變完全相同，屬于二級相變；而低于臨界點是的氣液相變屬于一級相變。在臨界點以上，氣體是不能夠通過等溫壓縮被轉(zhuǎn)變?yōu)橐合嗟?。系統(tǒng)在臨界點具有許多特殊性質(zhì)，稱為臨界現(xiàn)象。下面來介紹臨界點的確定。</p>&l

67、t;p>　　例:由范德瓦爾斯物態(tài)方程求臨界參量、、。</p><p><b>　　解題分析</b></p><p>　　從圖可以看出，臨界點是一拐點，它同時滿足下列條件：</p><p><b>　　，</b></p><p>　　利用上述拐點條件，將范德瓦爾斯方程對求導(dǎo)并聯(lián)解方程，便可求得

68、范德瓦爾斯氣體的三個臨界參量。</p><p><b>　　MATLAB程序</b></p><p><b>　　clear</b></p><p>　　syms a b R T</p><p>　　D1=diff('(p+a/v^2)*(v-b)-R*T','v')

69、;</p><p>　　D2=diff(D1,'v');</p><p>　　[pc,vc]=solve(D1,D2,'v','p')</p><p>　　Tc=solve((pc+a/vc^2)*(vc-b)-R*T,'T')</p><p><b>　　運(yùn)行結(jié)果&l

70、t;/b></p><p>　　pc=1/27*a/b^2</p><p><b>　　vc=3*b</b></p><p>　　Tc=8/27*a/b/R</p><p><b>　　即；，</b></p><p><b>　　3電磁學(xué)問題</b

71、></p><p>　　3.1求電偶極子在其所在平面產(chǎn)生的電場中任一點P的電位</p><p>　　例：已知電偶極子中兩電荷-q和+q的距離為。計算中可取。</p><p><b>　　解題分析</b></p><p>　　設(shè)場點P到的距離為和，則單獨存在時P點的電位分別為</p><p>

72、　　由電位疊加原理，電偶極子產(chǎn)生的電場在P點的電位為</p><p><b>　　MATLAB程序</b></p><p><b>　　clear;</b></p><p>　　q=1.6e-19; %單位電荷電量</p><p><b>　　C0=9e9;</b>&

73、lt;/p><p>　　l=3.0; %偶極子正負(fù)電荷之間的距離l</p><p>　　x=-5:0.5:5;y=x;</p><p>　　[X,Y]=meshgrid(x,y);</p><p>　　r1=sqrt((X-1./2).^2+(Y-0).^2); %電荷距離空間P（x，y）點的距離</p><

74、;p>　　r2=sqrt((X+1./2).^2+(Y-0).^2);</p><p>　　U=q.*C0.*(1./r1-1./r2); %求出空間任意一點P（x，y）的電位</p><p>　　plot(-1/2,0,'ro',-1/2,0,'r-') %標(biāo)出負(fù)電荷</p><p><b>　　

75、hold on,</b></p><p>　　plot(1/2,0,'ro',1/2,0,'r+') %標(biāo)出正電荷</p><p>　　C=contour(X,Y,U,'k-');</p><p>　　clabel(C); %畫等位線并標(biāo)出電位值</p>&l

76、t;p>　　axis('square')</p><p>　　運(yùn)行結(jié)果如圖3.1所示</p><p>　　3.1電位梯度與電場強(qiáng)度</p><p>　　電位是標(biāo)量，它在空間中每點都有一定的數(shù)值，所以電位是標(biāo)量場。標(biāo)量場在空間中沿不同方向的變化率稱為梯度，對電位場而言稱為電位梯度，用grad U或來表示?？梢宰C明，電位梯度和電場強(qiáng)度E的關(guān)系為&l

77、t;/p><p>　　利用上式，可從已知的電位分布求電場強(qiáng)度。</p><p>　　3.2由電位的表示式計算電場并畫出等電位線和電場方向</p><p><b>　　解題分析</b></p><p>　　如果已知空間的電位分布</p><p><b>　　則空間的電場強(qiáng)度為</b>

78、;</p><p>　　按照本題的要求，可利用讀入字符串的指令input('U'(x,y)=','s')來輸入電位方程。在MATLAB中，梯度函數(shù)的調(diào)用格式為gradient（），它是靠數(shù)值微分得到的。因此，空間觀測點應(yīng)取得密一些，以獲得較高的精度。</p><p><b>　　MATLAB程序</b></p>&

79、lt;p><b>　　clear all</b></p><p>　　U=input('請輸入電位方程，U=(x,y)=','s'); %例如，取U(x,y)=log(x.^2+y.^2)。</p><p>　　xmax=5;ymax=5;Ngrid=20; %繪圖區(qū)從x=-xmax到xmax，網(wǎng)格線數(shù)為20

80、</p><p>　　xplot=linspace(-xmax,xmax,Ngrid); %繪圖用x的數(shù)組</p><p>　　[x,y]=meshgrid(xplot); %x，y取同樣范圍，生成二維網(wǎng)格</p><p>　　Uplot=eval(U); %執(zhí)行輸入的字符串U,計算各點U的值&l

81、t;/p><p>　　[Explot,Eyplot]=gradient(-Uplot); %電場等于電位的負(fù)梯度</p><p>　　clf;subplot(1,2,1),meshc(Uplot); %劃分子圖；繪制含等位線的三維曲面</p><p>　　xlabel('x');ylabel('y');zlabel('U&

82、#39;);</p><p>　　subplot(1,2,2),</p><p>　　axis([-xmax,xmax,-ymax,ymax]); %規(guī)定等位線的范圍</p><p>　　cs=contour(x,y,Uplot); %畫等位線，cs是等位線值</p><p>　　clabel(cs);

83、 %標(biāo)出等位線的值</p><p><b>　　hold on;</b></p><p>　　quiver(x,y,Explot,Eyplot); %保持圖形，在原圖形上疊加矢量場圖</p><p>　　xlabel('x');ylabel('y');&

84、lt;/p><p><b>　　hold off;</b></p><p>　　圖3.2：的電位分布與電場分布</p><p>　　運(yùn)行上述程序，所得結(jié)果如圖3.2所示。</p><p>　　3.3帶電粒子在電磁場中的運(yùn)動</p><p>　　例：設(shè)質(zhì)量為m，帶電量為q的粒子在電磁場中的運(yùn)動微分方程為

85、</p><p>　　選場中某點為原點，以為方向，沿方向，建立坐標(biāo)系。令，上式的投影方程為</p><p>　　令，上述方程可改寫為下列一階微分方程組：</p><p><b>　　MATLAB程序</b></p><p>　　%%符號法求離子運(yùn)動微分方程的特解并繪圖</p><p><b&

86、gt;　　clear</b></p><p>　　syms w x y z t B E m q;</p><p>　　E=input('E=');B=input('B='); %輸入E和B值</p><p>　　[x,y,z]=dsolve('D2x=q*B/m*Dy','D2y=q*E/m

87、-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');</p><p>　　%初始條件取x(0)=y(0)=z(0)=0,Dx(0)=0.01,Dy(0)=0.01</p>

88、<p>　　q=1.6e-2;m=0.02</p><p>　　X=subs([x,y,z]);x=X(1),y=X(2),z=X(3),</p><p>　　ezplot3(X(1),X(2),X(3))</p><p>　　運(yùn)行上述程序，例如，取E=4，B=8可得下列特解并給出圖3.3</p><p>　　x =-15/16

89、*cos(32/5*t)-49/640*sin(32/5*t)+1/2*t+15/16</p><p>　　y =15/16*sin(32/5*t)-49/640*cos(32/5*t)+49/640</p><p>　　z =1/100*t</p><p>　　(a)E=4,B=8(b)E=0.01,B=8(c)E=8,B=1</p>

90、<p>　　圖3.3：帶電粒子在電磁場中的運(yùn)動</p><p>　　下面我們給出一段用數(shù)值方法求解該問題的程序，用以比較。</p><p>　　q=1.6e-2;m=0.02;</p><p>　　B=[2;2;0];E=[1;0;1];</p><p><b>　　figure</b></p>

91、<p>　　strd{1}='E\neq 0,B\neq 0';</p><p>　　strd{2}='E=0,B\neq 0';</p><p>　　strd{3}='E\neq 0,B=0';</p><p><b>　　for i=1:3</b></p><p&

92、gt;　　[t,y]=ode23('ex3f',[0:0.1:20],[0,0.01,0,6,0,0.01],[],q,m,B(i),E(i));</p><p>　　axes('unit','normalized','position',[0.0293+(i-1)*0.325 0.062 0.28 0.658]);</p><p

93、>　　plot3(y(:,1),y(:,3),y(:,5),'linewidth',2);</p><p><b>　　grid on</b></p><p>　　title(strd{i},'fontsize',12,'fontweight','demi');</p><p&g

94、t;　　view([-51,18]);</p><p><b>　　End</b></p><p>　　函數(shù)文件是一個獨立的文件，文件名為ex3f.m</p><p>　　function ydot=ex3f(t,y,flag,q,m,b,e)</p><p>　　ydot=[y(2);q*b*y(4)/m;y(4);q

95、*e/m-q*b*y(2)/m;y(6);0];</p><p>　　運(yùn)行該程序，可得到與上圖3.3相同的結(jié)果。</p><p><b>　　結(jié)論</b></p><p>　　從本文利用MATLAB語言對經(jīng)典物理一些具體問題的分析，并得出最終結(jié)論，首先，應(yīng)用MATLAB求解這些問題，使原來繁瑣的手工計算變得簡便，而且可將物理題中的解及一些特殊函

96、數(shù)以圖形的形式顯示出來，形象、直觀，便于理解。而且MATLAB強(qiáng)大的科學(xué)運(yùn)算、靈活的程序設(shè)計、便捷的與其他程序和語言接口的功能，顯示出很強(qiáng)的優(yōu)越性。其次，應(yīng)用MATLAB解決以上物理題方程的時候解決掉了手工計算式子多，計算繁雜，求解過程復(fù)雜的問題。得出的圖形直觀，對掌握物理問題有一定幫助。</p><p><b>　　參考文獻(xiàn)</b></p><p>　　[1] 孫祥

97、.MATLAB 7.0基礎(chǔ)教程.北京:清華大學(xué)出版社 2006</p><p>　　[2] 劉衛(wèi)國.MATLAB程序設(shè)計教程.北京:中國水利水電出版社.北京:2005</p><p>　　[3] 陳懷琛.MATLAB及其在理工課程中的應(yīng)用指南.西安:西安電子科技大學(xué)出版社，2000</p><p>　　[4] 張志涌.精通MATLAB6.5版.北京:北京航空航天大學(xué)

98、出版社，2004</p><p>　　[5] 彭芳麟.數(shù)學(xué)物理方程的MATLAB解法與可視化.北京:清華大學(xué)出版社，2005</p><p>　　[6] 薛定宇,陳陽泉.高等應(yīng)用數(shù)學(xué)問題的MATLAB求解.北京:清華大學(xué)出版社,2004</p><p>　　[7] 黃忠霖,黃京.MATLAB符號運(yùn)算及其應(yīng)用.北京:國防工業(yè)出版社,2004</p>&l

99、t;p>　　[8] 蘇金明,張蓮花等.MATLAB工具箱應(yīng)用.北京:電子工業(yè)出版社,2004</p><p>　　[9] 陸果.基礎(chǔ)物理學(xué).北京:高等教育出版社,1997</p><p>　　[10] 汪志城.熱力學(xué)?統(tǒng)計物理（第三版）.北京:高等教育出版社,2003</p><p>　　[11] 秦允豪.熱學(xué)(第二版).北京:高等教育出版社,2004.6&l

100、t;/p><p>　　[12] 漆安慎,杜嬋英.力學(xué)(第二版).北京:高等教育出版社,2005.6</p><p>　　[13] 梁燦彬,秦光戎等.電磁學(xué).北京:高等教育出版社,2004.5</p><p>　　[14] 盧德馨.大學(xué)物理學(xué).北京:高等教育出版社,1998</p><p>　　[15] 彭芳麟.理論力學(xué)的計算機(jī)模擬.北京:清華大學(xué)

101、出版社,2002</p><p><b>　　致謝</b></p><p>　　在此論文撰寫過程中，要特別感謝我的導(dǎo)師的指導(dǎo)與督促，感謝他的諒解與包容。還得感謝我的母校、感謝所有授予我知識、幫助我的老師。還有同學(xué)們的幫助。在你們的幫助下完成了這篇論文。謝謝我的父母，沒有他們辛勤的付出也就沒有我的今天，在這一刻，將最崇高的敬意獻(xiàn)給你們！</p><p