南邮应用物理计算物理实践报告

1、南京邮电大学实验报告课程名称：计算物理实践专 业：应用物理学学 号：姓 名：完成日期： 至a目 录 TOC o 1-5 h z 一、简单物理实验的模拟及实验数据处理 1 HYPERLINK l bookmark32 o Current Document 问题描述： 1单摆运动原理 1 HYPERLINK l bookmark16 o Current Document 模型的建立 1 HYPERLINK l bookmark18 o Current Document 流程图 2 HYPERLINK l bookmark20 o Current Document MAtlabS序设计仿真 2 HY

2、PERLINK l bookmark22 o Current Document Matlab 序 2 HYPERLINK l bookmark24 o Current Document 单摆演示截图 3二、方程组的数值解法 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序运行结果： 5 HYPERLINK l bookmark30 o Current Document 三、静电场问题的计算 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讯序设计仿真 10 HYPERLINK l bookmark42 o Current Document MATLAB仿真结果 10 HYPERLINK l bookmark44 o Current Document 四、热传导方程和波动方程的差分解法 11问题描述 11原理分析 11 HYPERLINK l bookmark56 o Current Document 具体步骤 13 HYPERLINK l bookmark58 o Current Document MATLA程序设计仿真 13 HYPERLINK l bookmark60 o Current Document MATL

5、ABS序运彳亍结果 13 HYPERLINK l bookmark62 o Current Document 结束语 15 HYPERLINK l bookmark64 o Current Document 参考文献 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一、简单物理实验的模拟及实验数据处理问题描述：编写单摆运动演示程序。在不考虑空气阻力和 日很小的假设下，单位质量小球做理想简谐运动，此时取 g=9.8 ,L=1,U =pi/4.单摆运动原理设在某一时刻，单摆的摆线偏离垂直线的角位移为，将重力mg分解为径向力F和切向力T ,则T的大小为mgsin日,切向加速度为二L 0 .根据牛顿第二定律 得方程 ma=m i= mg i 从而单摆运动的微分方程为0=e由于6=6 -+- .当日很小时， 0 0所以单摆的微分方程可表示为 0=0 上式表明，当日很小时，单摆的角加速度与角位移成正比，但方向相反，且方程的解可表本为日=A（ 小

7、）模型的建立建立物理模型，假设单摆运动过程中的小球中心点的坐标为 y根据几何关系x, y满足如下关系：.口 一在不考虑空气阻力和e很小的假设下日二日- 与式比较知A=日 小 CD - TOC o 1-5 h z 66所以，单摆的运动方程为018兀将日=pi/4,g=9.8,L=1 带入得，流程图图1.1程序流程图Matlab程序设计仿真通过set函数将变现和小球的图像句柄，加入 X、Y的动态参量，它们得轨迹 变化即为上面分析的轨迹方程。用line函数画出初始的位置，并将句柄分别给 sphere、lp ,然后每隔dt时间刷新一次。尝试对dt的设置发现，值在0.0005的效 果比较好。Matlab

8、 程序程序见附录11.7单摆演示截图图1.2单摆演示图二、方程组的数值解法问题描述：二分法求解方程x3+4x210= 0在区间1, 2内的根，精度自设原理分析二分法理论f(x)单调，f(a)f(b)0f(x)=0 在(a,b)有唯一根。设 f(x)在a,b上连续，f(x)=0在a,b上存在唯一解，且f(a)f(b)0, 记第一步，计算f( )f().若f(,否则x*,记)f()0,则 x*，记,对两种情形均有 x*，记 TOC o 1-5 h z 第k步，计算f( )f()0,则x*，记，否则x*，对这两种情形均有x*，记,k=1,2,所以数列满足I x*-|一( 一 )=)=.= 一 )=

9、一即 从而当k充分大， 且可由I x*- 控制精度分析求解令 f(x)=-x 3+4x210,f(x)在1, 2上连续，且 ff(2) 0, 则 f(x)=0 在1,2上有唯一解，记 a=1,b=2,x=1.5.然后计算 f(a)f(b), 若 f(a)f(b)0,x ，此时记a=1,b=1.5,否贝U x，记a=1.5 ,b=2.对两种情形均有x ,记x=,按照同样的方法依次向下计算，直到求出的相邻两个x的值之差绝对误差小于0.00005,即可求出结果Matlab 程序程序见附录2.Matlab程序运行结果：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三、静电场问题的计算问题描述：设两个同轴矩形金属槽如图 3-1所示，外金属槽电位为 0,内金属槽电位为100V,求内 电位分布，并绘出电位分布图。图3-1原理分析. 1简单迭代法1、对某一网格点设一初值，这个初值完全可以任意给定，称为初值电位。 虽然，问题的最终结果与初值无关，但若初值选择得当，则计算步

11、骤会 得到简化（当利用计算机来实现迭代计算时，为了简化程序，初值点为 一般可取值为零）。2、初值电位给定后，按一个固定顺序（点的顺序是从左到右，从下到上）依次计算每点的点位，即利用小-用围绕它的四个点的电压的平均值作为它的新值，当所有的点计算完后,用他们的新值代替旧值，即完成了一次迭代计算。然后再进行下一次迭 代计算，直到每一点计算的新值和旧值之差小于指定的范围为止。简单迭代法的特点是用之前一次迭代得到的网络点电位作为下一次迭代的初值。如在（i,j ）点在n+1次迭代时计算公式为：3.2.2有限差分:二维拉普拉斯方程_ 2_R 中（x, y）=中 XX + 中 yy = f （x, y）（ 1

12、）有限差分法的网格划分，通常采用完全有规律的分布方式，这样可使每个离 散点上得到相同形式的差分方程，有效的提高解题速度，经常采用的是正方形网 格划分。图3-2迭代法网格划分，h Lyi,j=i,j:y6 ；:y3h3j j中h.x36 ：xh3：i 1,j = i,jh-X1 : h22 |2 ex把以上四式相加，在相加的过程中, 果的精度为h的二次项。h的所有奇次方项都抵消了。得到的结cp设网格节点（i,j ）的电位为i，j ,其上下左右四个节点的电位分别为Q （D （D Q（pi,j% i,j,j i*j在h充分小的情况下，可以“为基点进行泰勒级数展开:=h - hi j - i j二 y

13、 ：y(2)二 2二：2d4 P +h2( -)i,j + i ,j 1 TjFj4 i, j h (取2cy2由于场中任意点（i, j）都满足泊松方程：2二2二2二= F（x,y） TOC o 1-5 h z : x2 y式中F（x, y）为场源，则式（2）可变为：h2九=4Wi,jJ*,j J*,j+j）-zF（x,y）（3）对于无源场，F（x, y） = 0 ,则二维拉普拉斯方程的有限差分形式为：1（二）j . （ i，j 书 i i，j i 书,j ）（4）4上式表示任一点的电位等于围绕它的四个等间距点的电位的平均值，距离h越小则结果越精确，用式（4）可以近似的求解二维拉普拉斯方程。边

14、界条件：中（x, y） | 内槽=100;甲（x, y） |外槽二03.2.3解题过程解：在直角坐标系中，金属槽中的电位函数小满足拉普拉斯方程：其边界条件满足混合型编值问题的边界条件：取步长h=1,x、y方向的网格数为 m=16,n=10，共有16*10=160个网孔17*11=187个节点，其中槽内节点（电位待求点）有 15*9=135个，界节点（电位已知点）有187135=52个。设迭代精度为。3.3Matlab程序设计仿真源程序见附录三3.4Matlab仿真结果等电位图图3-3运行结果图四、热传导方程和波动方程的差分解法4.1问题描述2cu u u彳=7T20 x 1,0 t二 t .

15、X求热传导方程混合问题：U（X,0）=X20EXE1的数值解，其u（0,t）=0,u（1,t）=10tI中N、h、k值等参数自取（将计算结果图形化）。4.2原理分析二维热传导方程的初、边值混合问题与一维的相似，在确定差分方程格式并 给出定解条件后，按时间序号分层计算，只是每一层是由二维点阵组成，通常称 为网格。各向同性介质中无热源的二维热传导方程为：(0 xl , 0ys,0ys,0tT).初始条件是：U（X，y，0） =P（X，y） （2设时间步长为p,空间步长为 h ,二维平面Xoy分为M N的网格，并使Nh = l , Mh = s, t = kw （ k = 0,1,2.），则有 x

16、= ih （ i = 0,1,2，N）,对节 点（i，j）,在k时刻（即已时刻）有：-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就所设定的具体问题来讨论边界条件在x = 0边界的0y MlM2hyM

17、h区域以及整个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);% 计算梯度%quiver(xx,yy,Gx,Gy,-0.5,r)%根据梯度数画箭头axis(-2,hx+3,-2,hy+3)画外框边线画内框边界线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);%外框上边界标注外框下边界标注外框左边界标注外框右边界标注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