南郵應用物理計算物理實踐

....南京郵電大學實驗報告課程名稱：計算物理實踐專業(yè)：應用物理學學號：姓名：完成日期：年月

目錄一、簡單物理實驗的模擬及實驗數(shù)據(jù)處理 11.1問題描述： 11.2單擺運動原理 11.3模型的建立 11.4流程圖 21.5Matlab程序設計仿真 21.6Matlab程序 31.7單擺演示截圖 3二、方程組的數(shù)值解法 42.1問題描述： 42.2原理分析 42.2.1二分法理論 42.2.2分析求解 42.3Matlab程序 52.4Matlab程序運行結(jié)果： 5三、靜電場問題的計算 73.1問題描述： 73.2原理分析 73.2．1簡單迭代法 73.2.2有限差分： 83.2.3解題過程 93.3Matlab程序設計仿真 93.4Matlab仿真結(jié)果 10四、熱傳導方程和波動方程的差分解法 114.1問題描述 114.2原理分析 114.3具體步驟 134.4MATLAB程序設計仿真 134.5MATLAB程序運行結(jié)果 13結(jié)束語 15參考文獻 16附錄1： 17附錄2： 18附錄3： 19附錄4： 21..一、簡單物理實驗的模擬及實驗數(shù)據(jù)處理1.1問題描述：編寫單擺運動演示程序。在不考慮空氣阻力和很小的假設下.單位質(zhì)量小球做理想簡諧運動.此時。取g=9.8.L=1,0=pi/4.1.2單擺運動原理設在某一時刻.單擺的擺線偏離垂直線的角位移為.將重力mg分解為徑向力F和切向力T.則T的大小為mgsin,切向加速度為a=Ld2dt2.根據(jù)牛頓第二定律得方程ma=mLd從而單擺運動的微分方程為d2dt由于sin=—33!+55當很小時.sin≈所以單擺的微分方程可表示為d2dt2=上式表明.當很小時.單擺的角加速度與角位移成正比.但方向相反.且方程的解可表示為=Acos（ωt+φ）②1.3模型的建立建立物理模型.假設單擺運動過程中的小球中心點的坐標為x,y。根據(jù)幾何關系x,y滿足如下關系：x在不考慮空氣阻力和很小的假設下=0costgL與②所以.單擺的運動方程為x=L將0=pi/4,g=9.8,L=1帶入得.x=1.4流程圖開始開始輸入0輸入0t=0,dt=0.005t=0,dt=0.005t=1?Yt=1?Y=-L*cosAngle*cosY=-L*cosAngle*X=L*sin(Angle*結(jié)束t=t+dt結(jié)束t=t+dt圖1.1程序流程圖1.5Matlab程序設計仿真通過set函數(shù)將變現(xiàn)和小球的圖像句柄.加入X、Y的動態(tài)參量.它們得軌跡變化即為上面分析的軌跡方程。用line函數(shù)畫出初始的位置.并將句柄分別給sphere、lp.然后每隔dt時間刷新一次。嘗試對dt的設置發(fā)現(xiàn).值在0.0005的效果比較好。1.6Matlab程序程序見附錄1。1.7單擺演示截圖圖1.2單擺演示圖二、方程組的數(shù)值解法2.1問題描述：二分法求解方程x3+4x2－10＝0在區(qū)間［1.2］內(nèi)的根.精度自設。2.2原理分析2.2.1二分法理論f(x)∈Ca,b,單調(diào).f(a)f(b)<0f(x)=0在（a,b）有唯一根。設f(x)在[a,b]上連續(xù).f(x)=0在[a,b]上存在唯一解.a第一步.計算f(a0)f(x0).若f(a0)f(x0)<0,則 x*∈(a0,b0)第k步.計算f(ak-1)f(xk-1)<0,則x*∈[ak-1,xk-1],記ak?k,x*∈[ak,bk]且xk=ak+bk2.所以數(shù)列｛即xk→x*,從而當k充分大.x*≈xk且可由2.2.2分析求解令f(x)=-x3+4x2－10,f(x)在［1.2］上連續(xù).且f(1)f(2)<0,則f(x)=0在［1.2］上有唯一解.記a=1,b=2,x=1.5.然后計算f(a)f(b),若f(a)f(b)<0.則x∈(1,1.5),此時記a=1,b=1.5,否則x∈(1.5,2),記a=1.5,b=2.對兩種情形均有x∈[a,b],記x=a+b2,按照同樣的方法依次向下計算.直到求出的相鄰兩個x的值之差絕對誤差小于0.00005.2.3Matlab程序程序見附錄2.2.4Matlab程序運行結(jié)果：x=1.2500x=1.3750x=1.3125x=1.3438x=1.3594x=1.3672x=1.3633x=1.3652x=1.3643x=1.3647x=1.3650x=1.3651x=1.3652x=1.3652x=1.3652x=1.3652x=1.3652x=1.36523,f(x)=0.00001三、靜電場問題的計算3.1問題描述：設兩個同軸矩形金屬槽如圖3-1所示.外金屬槽電位為0.內(nèi)金屬槽電位為100V.求內(nèi)電位分布.并繪出電位分布圖。圖3-13.2原理分析3.2．1簡單迭代法對某一網(wǎng)格點設一初值.這個初值完全可以任意給定.稱為初值電位。雖然.問題的最終結(jié)果與初值無關.但若初值選擇得當.則計算步驟會得到簡化（當利用計算機來實現(xiàn)迭代計算時.為了簡化程序.初值點為一般可取值為零）。初值電位給定后.按一個固定順序（點的順序是從左到右.從下到上）依次計算每點的點位.即利用φi,j=14(φi,j+1+φi,j-1簡單迭代法的特點是用之前一次迭代得到的網(wǎng)絡點電位作為下一次迭代的初值。如在（i,j）點在n+1次迭代時計算公式為：φ3.2.2有限差分：二維拉普拉斯方程（1）有限差分法的網(wǎng)格劃分.通常采用完全有規(guī)律的分布方式.這樣可使每個離散點上得到相同形式的差分方程.有效的提高解題速度.經(jīng)常采用的是正方形網(wǎng)格劃分。圖3-2迭代法網(wǎng)格劃分設網(wǎng)格節(jié)點（i,j）的電位為.其上下左右四個節(jié)點的電位分別為在h充分小的情況下.可以為基點進行泰勒級數(shù)展開：把以上四式相加.在相加的過程中.h的所有奇次方項都抵消了。得到的結(jié)果的精度為h的二次項。（2）由于場中任意點都滿足泊松方程：式中為場源.則式（2）可變?yōu)椋海?）對于無源場..則二維拉普拉斯方程的有限差分形式為：（4）上式表示任一點的電位等于圍繞它的四個等間距點的電位的平均值.距離h越小則結(jié)果越精確.用式（4）可以近似的求解二維拉普拉斯方程。邊界條件：|內(nèi)槽=100；|外槽=03.2.3解題過程解：在直角坐標系中.金屬槽中的電位函數(shù)φ滿足拉普拉斯方程：?其邊界條件滿足混合型編值問題的邊界條件：φ?φ取步長h=1,x、y方向的網(wǎng)格數(shù)為m=16,n=10,共有16*10=160個網(wǎng)孔17*11=187個節(jié)點.其中槽內(nèi)節(jié)點（電位待求點）有15*9=135個.界節(jié)點（電位已知點）有187—135=52個。設迭代精度為103.3Matlab程序設計仿真源程序見附錄三3.4Matlab仿真結(jié)果圖3-3運行結(jié)果圖四、熱傳導方程和波動方程的差分解法4.1問題描述求熱傳導方程混合問題：的數(shù)值解.其中N、ｈ、k值等參數(shù)自?。▽⒂嬎憬Y(jié)果圖形化）。4.2原理分析二維熱傳導方程的初、邊值混合問題與一維的相似.在確定差分方程格式并給出定解條件后.按時間序號分層計算.只是每一層是由二維點陣組成.通常稱為網(wǎng)格。各向同性介質(zhì)中無熱源的二維熱傳導方程為：.(,,,).①初始條件是：②設時間步長為τ.空間步長為h.二維平面xoy分為M×N的網(wǎng)格.并使..（）.則有（）.對節(jié)點.在時刻（即時刻）有：③將差分格式③代入偏微分方程①中.可得：④（i=1,2,…,N-1,j=1,2,就所設定的具體問題來討論邊界條件:在邊界的和區(qū)域以及整個.邊界均為絕熱壁；而在邊界的區(qū)域為與恒溫熱源相連的口。和兩邊界溫度始終為0.實際上也是與恒溫源相連的。也就是說.對于絕熱壁應滿足：,（上述邊界條件的差分近似式為：即：（）對于與恒溫源相連的邊界.在熱傳導過程中始終有恒定的熱流.?？扇w一化值.例如高溫熱源可取“1”.而低溫熱源可取“0”。對于這種情況邊界條件還有：綜合上述初值、邊值混合問題.并設初始時刻各點溫度均為零.則上述差分格式可歸納為：u(i,k+1)=a*u(i+1,k)+(1-2*a)*u(i,k)+a*u(i-1,k)u(0,k)=0u(1,k)=1可以證明.對于二維熱傳導方程.若滿足.則差分格式式④或式⑥就是穩(wěn)定的差分格式.一般的講.對于n維拋物線型微分方程差分格式穩(wěn)定的充分條件是：。4.3具體步驟1.給定、、和.取h=0.1.2.計算初值和邊值：；；；；；用差分格式計算；4.4MATLAB程序設計仿真源程序見附錄四MATLAB程序運行結(jié)果圖4-1運行結(jié)果圖結(jié)束語計算物理學是利用電子計算機進行數(shù)據(jù)采集、數(shù)值計算和數(shù)字仿真來發(fā)現(xiàn)和研究物理現(xiàn)象與物理規(guī)律的一門現(xiàn)代交叉學科。經(jīng)過本次計算物理學實驗周的學習.我學會了怎樣利用計算機來進行數(shù)據(jù)分析和用Matlab來進行數(shù)學以及物理問題的仿真而得到問題的結(jié)果。我認識到自己對于以前學習過的一些課程掌握得還不夠.Matlab編程語言的運用也不夠熟練。通過這次實驗也很好的鞏固了以前學習的一些知識點.將理論運用于實踐。這次實驗讓我認識到數(shù)理方程的實用性.掌握了利用差分代替微分來求解波動方程、熱傳導方程、拉普拉斯方程等的基本原理和方法。本次實踐涉及到的二維拉普拉斯方程以及二維熱傳導方程的解題方法.都是先將連續(xù)的方程以及邊界條件離散化.再用計算機進行計算.因為計算機智能對離散的數(shù)值進行計算。對于非線性方程的求解往往是采用迭代的方法求解.本次實踐主要涉及了Newton迭代法的重要思想.也是將連續(xù)的方程離散化后再進行計算。通過本次實踐.我對Matlab的使用更加熟練了.解決問題以及整理數(shù)據(jù)等能力有很大提高.希望我們在學?？梢杂懈嗟臋C會實踐.有更多的收獲。參考文獻[1]陳鍾賢.計算物理學.哈爾濱工業(yè)大學出版社.2001.3[2]楊振華.酈志新.數(shù)學實驗.科學出版社.2010.2[3]林亮.吳群英.數(shù)值分析方法與實驗：基于MATLAB實現(xiàn).高等教育出版社.2012.9[4]李慶楊.王能超.易大義.數(shù)值分析.華中科技大學出版社.2006.7[5]鐘季康.鮑鴻吉.大學物理習題計算機解法—MATLAB編程應用.機械工業(yè)出版社.2008.1[6]何紅雨.電磁場數(shù)值計算法與MATLAB實現(xiàn).華中科技大學出版社附錄1：plot([-0.2;0.2],[0;0],'color','y','linestyle','-','linewidth',10);g=9.8;l=1;theta0=pi/4;x0=l*sin(theta0);y0=(-1)*l*cos(theta0);axis([-0.75,0.75,-1.25,0]);axis('on');head=line(x0,y0,'color','r','linestyle','.','erasemode','xor','markersize',40);body=line([0;x0],[0;y0],'color','b','linestyle','-','erasemode','xor');t=0;dt=0.005;while1t=t+dttheta=theta0*cos(sqrt(g/l)*t);x=l*sin(theta);y=(-a)*l*cos(theta);set(head,'xdata',x,'ydata',y);set(body,'xdata',[0;x],'ydata',[0;y]);drawnow;end附錄2：a=1;b=2;f=@(x)x^3+4*x^2-10;c=(a+b)/2;whileabs(b-a)>1e-5iff(c)*f(b)<0a=c;elseb=c;endc=(a+b)/2;x=cendfprintf('\nx=%.5f,f(x)=%.5f\n',x,f(x));附錄3：FI=100;hx=23;hy=23;%定外框網(wǎng)絡；x軸1~23,1~23CX1=9;CX2=15;CY1=9;CY2=15;%定內(nèi)框位置：x軸9~15.y軸9~15DX=1+CX2-CX1;DY=1+CY2-CY1;v=ones(hy,hx)*50;%每個網(wǎng)格點的初值mesh1=ones(hy,hx)*2;%每個網(wǎng)格點的計算標志值v(1,:)=zeros(1,hx);%第一行的初值（下邊界）mesh1(1,:)=zeros(1,hx);%第一行標志值為零（下邊界）v(hy,:)=zeros(1,hx);%上邊界初值mesh1(hy,:)=zeros(1,hx);%上邊界標志值為零v(:,1)=zeros(hy,1);%左邊界初值mesh1(:,1)=zeros(hy,1);%左邊界標志值為零v(:,hx)=zeros(hy,1);%右邊界初值mesh1(:,hx)=zeros(hy,1);%右邊界標志值為零v(CY1:CY2,CY1:CY2)=ones(DX,DY)*FI;%內(nèi)框區(qū)初值mesh1(CY1:CY2,CX1:CX2)=ones(DX,DY);%內(nèi)框區(qū)標志值k=0;difmax=1.0;%最大誤差while(difmax>1.0e-6)k=k+1;%計算迭代次數(shù)difmax=0.0;fori=2:hy-1%從2到20行循環(huán)forj=2:hx-1%從2到20列循環(huán)m=mesh1(i,j);%?。╥,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(dif>difmax)difmax=dif;e