南郵應(yīng)用物理計(jì)算物理實(shí)踐報(bào)告

1、南京郵電大學(xué)實(shí)驗(yàn)報(bào)告課程名稱：計(jì)算物理實(shí)踐專 業(yè)：應(yīng)用物理學(xué)學(xué) 號：姓 名：完成日期： 至a目 錄 TOC o 1-5 h z 一、簡單物理實(shí)驗(yàn)的模擬及實(shí)驗(yàn)數(shù)據(jù)處理 1 HYPERLINK l bookmark32 o Current Document 問題描述： 1單擺運(yùn)動原理 1 HYPERLINK l bookmark16 o Current Document 模型的建立 1 HYPERLINK l bookmark18 o Current Document 流程圖 2 HYPERLINK l bookmark20 o Current Document MAtlabS序設(shè)計(jì)仿真 2 HY

2、PERLINK l bookmark22 o Current Document Matlab 序 2 HYPERLINK l bookmark24 o Current Document 單擺演示截圖 3二、方程組的數(shù)值解法 4 HYPERLINK l bookmark46 o Current Document 問題描述： 4 HYPERLINK l bookmark48 o Current Document 原理分析 4二分法理論 4分析求解 4 HYPERLINK l bookmark26 o Current Document MAtlabS序 5 HYPERLINK l bookmark2

3、8 o Current Document MAtlabS序運(yùn)行結(jié)果： 5 HYPERLINK l bookmark30 o Current Document 三、靜電場問題的計(jì)算 7問題描述： 7 HYPERLINK l bookmark34 o Current Document 原理分析 7 HYPERLINK l bookmark36 o Current Document . 1簡單迭代法 7 HYPERLINK l bookmark38 o Current Document 有限差分： 8解題過程 9 HYPERLINK l bookmark40 o Current Document M

4、ATLA訊序設(shè)計(jì)仿真 10 HYPERLINK l bookmark42 o Current Document MATLAB仿真結(jié)果 10 HYPERLINK l bookmark44 o Current Document 四、熱傳導(dǎo)方程和波動方程的差分解法 11問題描述 11原理分析 11 HYPERLINK l bookmark56 o Current Document 具體步驟 13 HYPERLINK l bookmark58 o Current Document MATLA程序設(shè)計(jì)仿真 13 HYPERLINK l bookmark60 o Current Document MATL

5、ABS序運(yùn)彳亍結(jié)果 13 HYPERLINK l bookmark62 o Current Document 結(jié)束語 15 HYPERLINK l bookmark64 o Current Document 參考文獻(xiàn) 16 HYPERLINK l bookmark66 o Current Document 附錄1: 17 HYPERLINK l bookmark68 o Current Document 附錄2: 18 HYPERLINK l bookmark70 o Current Document 附錄3: 19 HYPERLINK l bookmark72 o Current Docum

6、ent 附錄4: 21一、簡單物理實(shí)驗(yàn)的模擬及實(shí)驗(yàn)數(shù)據(jù)處理問題描述：編寫單擺運(yùn)動演示程序。在不考慮空氣阻力和 日很小的假設(shè)下，單位質(zhì)量小球做理想簡諧運(yùn)動，此時(shí)取 g=9.8 ,L=1,U =pi/4.單擺運(yùn)動原理設(shè)在某一時(shí)刻，單擺的擺線偏離垂直線的角位移為，將重力mg分解為徑向力F和切向力T ,則T的大小為mgsin日,切向加速度為二L 0 .根據(jù)牛頓第二定律 得方程 ma=m i= mg i 從而單擺運(yùn)動的微分方程為0=e由于6=6 -+- .當(dāng)日很小時(shí)， 0 0所以單擺的微分方程可表示為 0=0 上式表明，當(dāng)日很小時(shí)，單擺的角加速度與角位移成正比，但方向相反，且方程的解可表本為日=A（ 小

7、）模型的建立建立物理模型，假設(shè)單擺運(yùn)動過程中的小球中心點(diǎn)的坐標(biāo)為 y根據(jù)幾何關(guān)系x, y滿足如下關(guān)系：.口 一在不考慮空氣阻力和e很小的假設(shè)下日二日- 與式比較知A=日 小 CD - TOC o 1-5 h z 66所以，單擺的運(yùn)動方程為018兀將日=pi/4,g=9.8,L=1 帶入得，流程圖圖1.1程序流程圖Matlab程序設(shè)計(jì)仿真通過set函數(shù)將變現(xiàn)和小球的圖像句柄，加入 X、Y的動態(tài)參量，它們得軌跡 變化即為上面分析的軌跡方程。用line函數(shù)畫出初始的位置，并將句柄分別給 sphere、lp ,然后每隔dt時(shí)間刷新一次。嘗試對dt的設(shè)置發(fā)現(xiàn)，值在0.0005的效 果比較好。Matlab

8、 程序程序見附錄11.7單擺演示截圖圖1.2單擺演示圖二、方程組的數(shù)值解法問題描述：二分法求解方程x3+4x210= 0在區(qū)間1, 2內(nèi)的根，精度自設(shè)原理分析二分法理論f(x)單調(diào)，f(a)f(b)0f(x)=0 在(a,b)有唯一根。設(shè) f(x)在a,b上連續(xù)，f(x)=0在a,b上存在唯一解，且f(a)f(b)0, 記第一步，計(jì)算f( )f().若f(,否則x*,記)f()0,則 x*，記,對兩種情形均有 x*，記 TOC o 1-5 h z 第k步，計(jì)算f( )f()0,則x*，記，否則x*，對這兩種情形均有x*，記,k=1,2,所以數(shù)列滿足I x*-|一( 一 )=)=.= 一 )=

9、一即 從而當(dāng)k充分大， 且可由I x*- 控制精度分析求解令 f(x)=-x 3+4x210,f(x)在1, 2上連續(xù)，且 ff(2) 0, 則 f(x)=0 在1,2上有唯一解，記 a=1,b=2,x=1.5.然后計(jì)算 f(a)f(b), 若 f(a)f(b)0,x ，此時(shí)記a=1,b=1.5,否貝U x，記a=1.5 ,b=2.對兩種情形均有x ,記x=,按照同樣的方法依次向下計(jì)算，直到求出的相鄰兩個(gè)x的值之差絕對誤差小于0.00005,即可求出結(jié)果Matlab 程序程序見附錄2.Matlab程序運(yùn)行結(jié)果：x = 1.2500 x =1.3750 x = 1.3125x = 1.3438x

10、 = 1.3594x =1.3672x = 1.3633x = 1.3652x = 1.3643x = 1.3647x = 1.3650 x =1.3651x =1.3652x = 1.3652x = 1.3652x =1.3652x = 1.3652x=1.36523,f(x)=0.00001三、靜電場問題的計(jì)算問題描述：設(shè)兩個(gè)同軸矩形金屬槽如圖 3-1所示，外金屬槽電位為 0,內(nèi)金屬槽電位為100V,求內(nèi) 電位分布，并繪出電位分布圖。圖3-1原理分析. 1簡單迭代法1、對某一網(wǎng)格點(diǎn)設(shè)一初值，這個(gè)初值完全可以任意給定，稱為初值電位。 雖然，問題的最終結(jié)果與初值無關(guān)，但若初值選擇得當(dāng)，則計(jì)算步

11、驟會 得到簡化（當(dāng)利用計(jì)算機(jī)來實(shí)現(xiàn)迭代計(jì)算時(shí)，為了簡化程序，初值點(diǎn)為 一般可取值為零）。2、初值電位給定后，按一個(gè)固定順序（點(diǎn)的順序是從左到右，從下到上）依次計(jì)算每點(diǎn)的點(diǎn)位，即利用小-用圍繞它的四個(gè)點(diǎn)的電壓的平均值作為它的新值，當(dāng)所有的點(diǎn)計(jì)算完后,用他們的新值代替舊值，即完成了一次迭代計(jì)算。然后再進(jìn)行下一次迭 代計(jì)算，直到每一點(diǎn)計(jì)算的新值和舊值之差小于指定的范圍為止。簡單迭代法的特點(diǎn)是用之前一次迭代得到的網(wǎng)絡(luò)點(diǎn)電位作為下一次迭代的初值。如在（i,j ）點(diǎn)在n+1次迭代時(shí)計(jì)算公式為：3.2.2有限差分:二維拉普拉斯方程_ 2_R 中（x, y）=中 XX + 中 yy = f （x, y）（ 1

12、）有限差分法的網(wǎng)格劃分，通常采用完全有規(guī)律的分布方式，這樣可使每個(gè)離 散點(diǎn)上得到相同形式的差分方程，有效的提高解題速度，經(jīng)常采用的是正方形網(wǎng) 格劃分。圖3-2迭代法網(wǎng)格劃分，h Lyi,j=i,j:y6 ；:y3h3j j中h.x36 ：xh3：i 1,j = i,jh-X1 : h22 |2 ex把以上四式相加，在相加的過程中, 果的精度為h的二次項(xiàng)。h的所有奇次方項(xiàng)都抵消了。得到的結(jié)cp設(shè)網(wǎng)格節(jié)點(diǎn)（i,j ）的電位為i，j ,其上下左右四個(gè)節(jié)點(diǎn)的電位分別為Q （D （D Q（pi,j% i,j,j i*j在h充分小的情況下，可以“為基點(diǎn)進(jìn)行泰勒級數(shù)展開:=h - hi j - i j二 y

13、 ：y(2)二 2二：2d4 P +h2( -)i,j + i ,j 1 TjFj4 i, j h (取2cy2由于場中任意點(diǎn)（i, j）都滿足泊松方程：2二2二2二= F（x,y） TOC o 1-5 h z : x2 y式中F（x, y）為場源，則式（2）可變?yōu)椋篽2九=4Wi,jJ*,j J*,j+j）-zF（x,y）（3）對于無源場，F(xiàn)（x, y） = 0 ,則二維拉普拉斯方程的有限差分形式為：1（二）j . （ i，j 書 i i，j i 書,j ）（4）4上式表示任一點(diǎn)的電位等于圍繞它的四個(gè)等間距點(diǎn)的電位的平均值，距離h越小則結(jié)果越精確，用式（4）可以近似的求解二維拉普拉斯方程。邊

14、界條件：中（x, y） | 內(nèi)槽=100;甲（x, y） |外槽二03.2.3解題過程解：在直角坐標(biāo)系中，金屬槽中的電位函數(shù)小滿足拉普拉斯方程：其邊界條件滿足混合型編值問題的邊界條件：取步長h=1,x、y方向的網(wǎng)格數(shù)為 m=16,n=10，共有16*10=160個(gè)網(wǎng)孔17*11=187個(gè)節(jié)點(diǎn)，其中槽內(nèi)節(jié)點(diǎn)（電位待求點(diǎn)）有 15*9=135個(gè)，界節(jié)點(diǎn)（電位已知點(diǎn)）有187135=52個(gè)。設(shè)迭代精度為。3.3Matlab程序設(shè)計(jì)仿真源程序見附錄三3.4Matlab仿真結(jié)果等電位圖圖3-3運(yùn)行結(jié)果圖四、熱傳導(dǎo)方程和波動方程的差分解法4.1問題描述2cu u u彳=7T20 x 1,0 t二 t .

15、X求熱傳導(dǎo)方程混合問題：U（X,0）=X20EXE1的數(shù)值解，其u（0,t）=0,u（1,t）=10tI中N、h、k值等參數(shù)自取（將計(jì)算結(jié)果圖形化）。4.2原理分析二維熱傳導(dǎo)方程的初、邊值混合問題與一維的相似，在確定差分方程格式并 給出定解條件后，按時(shí)間序號分層計(jì)算，只是每一層是由二維點(diǎn)陣組成，通常稱 為網(wǎng)格。各向同性介質(zhì)中無熱源的二維熱傳導(dǎo)方程為：(0 xl , 0ys,0ys,0tT).初始條件是：U（X，y，0） =P（X，y） （2設(shè)時(shí)間步長為p,空間步長為 h ,二維平面Xoy分為M N的網(wǎng)格，并使Nh = l , Mh = s, t = kw （ k = 0,1,2.），則有 x

16、= ih （ i = 0,1,2，N）,對節(jié) 點(diǎn)（i，j）,在k時(shí)刻（即已時(shí)刻）有：-Ui,j,k _ Ui, j,k 1 - Ui,j ,k.：t.2二 Ui,j,k _ Ui 1,j,k - 2ui, j,k Ui,j,k TOC o 1-5 h z _2 -r2xh二 Ui, j,k = Ui,j 1,k - 2Ui, j,k , Ui, j,k2 一2yh將差分格式代入偏微分方程中，可得:Ui,j,k 十 一(14a )Ui,j,k ( (Ui +1, j ,k UTj,k Ui,j”k Ui, j -1,k )h2就所設(shè)定的具體問題來討論邊界條件在x = 0邊界的0y MlM2hyM

17、h區(qū)域以及整個(gè)x=Nh 0 y Mh邊界均為絕熱壁；而在x = 0邊界的M1hy1e-5if f*f(b)1.0e-6)k=k+1;%difmax=0.0; for i=2:hy-1%for j=2:hx-1%m=mesh1(i,j);%if(m=2)%vold=v(i,j);%v(i,j)=(1/4)*(v(i-1,j)+v(i,j-1)+v(i+1,j)+v(i,j+1);%方程差分式dif=v(i,j)-vold;%差dif=abs(dif);%if(difdifmax) difmax=dif; end %endendendendsubplot(1,2,1),mesh(v)% 畫三維曲面

18、圖axis(-2,hx+3,-2,hy+3,0,100)subplot(1,2,2),contour(v,13)%畫等電位線圖hold onx=1:1:hx;y=1:1:hy;xx,yy=meshgrid(x,y);%形成柵格Gx,Gy=gradient(v,0.6,0.6);% 計(jì)算梯度%quiver(xx,yy,Gx,Gy,-0.5,r)%根據(jù)梯度數(shù)畫箭頭axis(-2,hx+3,-2,hy+3)畫外框邊線畫內(nèi)框邊界線plot(1,1,hx,hx,1,1,hy,hy,1,1,k)plot(CX1,CX1,CX2,CX2,CX1,CY1,CY2,CY2,CY1,CY1,k)% text(CX1+0.6,CY1+(CY2-CY1)/2,U=100,fontsize,10);%外框上邊界標(biāo)注外框下邊界標(biāo)注外框左邊界標(biāo)注外框右邊界標(biāo)注text(hx/2,hy+1,U=0,fontsize,10);% text(hx/2,0,U=0,fontsize,10);% text(-1.7,hy/2,U=0,fontsize,10);% text(hx+04hy