应用物理学 计算物理实践.pdf

南京邮电大学实验报告课程名称:计算物理实践专 业:应用物理学学 号:姓 名:完成日期:年 月目 录一、简单物理实验的模拟及实验数据处理.11.1问题描述：.11.2单摆运动原理.11.3 模型的建立.21.4 流程图.51.5 MATLAB程序设计仿真.71.6 MATLAB 程序.71.7 单摆演示截图.8二、方程组的数值解法.92.1问题描述：.42.2原理分析.92.2.1二分法理论.92.2.2分析求解.152.3 MATLAB 程序.162.4 MATLAB程序运行结果：.5三、静电场问题的计算.193.1 问题描述：.73.2原理分析.193.2.1简单迭代法.73.2.2有限差分:.83.2.3解题过程.213.3MATLAB程序设计仿真.223.4MATLAB仿真结果.22四、热传导方程和波动方程的差分解法.244.1问题描述.244.2原理分析.244.3具体步骤.254.4 MATLAB程序设计仿真.254.5MATLAB程序运行结果.25结束语.27参考文献.28附录1：.17附录2:.18附录3:.19附录4:21一、简单物理实验的模拟及实验数据处理1.1 问题描述:编写单摆运动演示程序。在不考虑空气阻力和很小的假设下，单位质 量小球做理想简谐运动，此时。取g=9.8,L=l,0=pi/4.1.2 单摆运动原理设在某一时刻，单摆的摆线偏离垂直线的角位移为，将重力mg分解为径向力F 和切向力T,则T的大小为mg sin,切向加速度为a=Ld2d t2.根据牛顿第二定律得方程ma二mLd2d t21二一mgsin从而单摆运动的微分方程为 M。二一分反。dtz L3 5由于sin 6=-+-.当很小时，sin所以单摆的微分方程可表示为夕二一匐 atz L上式表明，当很小时，单摆的角加速度与角位移成正比，但方向相反，且方程的解可表示为8二Acos(a)t+p)1.3模型的建立建立物理模型，假设单摆运动过程中的小球中心点的坐标为X,y。根据几何关系x,y满足如下关系：x2+y2L2 tany2X在不考虑空气阻力和。很小的假设下。二60cos0后 与式比较知 A二%,(p=OM=J|所以，单摆的运动方程为x=L sin=L sin 0cos tgLy=L cosL3cos0costgL4将0=pi/4,g=9.8,L=1 带入得,x=sinn4cost9.85y=_cosn4cost9.861.4流程图图i.i 程序流程图1.5 Mat lab程序设计仿真通过set函数将变现和小球的图像句柄，加入X、Y的动态参量，它们得轨迹变化即为上面分析的轨迹方程。用line函数画出初始的位置，并将句柄分别给 spherev Ip,然后每隔dt时间刷新一次。尝试对dt的设置发现，值在0.0005 的效果比较好。1.6 Matlab程序程序见附录1。71.7单摆演示截图图1.2单摆演示图8二 方程组的数值解法2.1问题描述：二分法求解方程x3+4x2 10=0在区间1,2内的根，精度自设。2.2原理分析2.2.1二分法理论f(x)Ca,b，单调，f(a)f(b)0f(x)=0在(a,b)有唯一根。设f(x)在a,b上连续，6)二0在上存在唯一解，且 f(a)f(b)0,记CLQ CLfb()=bf劭十为%o=第一步,计算f()f().若 f()f(0)0,则x*(900），记1010，否则X*（00），记1010，对两种情形均有X*100,记02第k步，计算f(-1)f(-1)0,则 x*-1-1，记1111，否则x*11，对这两种情形均有x*e，记+2,k=l,2,V k,x*e 12且2,所以数列满足I x*-12121321-1)=12+100)=12+1(b一 a)-014即X*,从而当k充分大，x*=且可由I x*-区12+1(b-a)控制精度。2.2.2分析求解15令6)=-*3+4*210代&)在1,2上连续，且ff包贝lf(x)=O在1,2上有唯一解,记a=l,b=2,x=l.5.然后计算f(a)f，若f(a)f(b)0,则x (1,1.5),此时记a=l,b=l.5,否则x(1.5,2),记a=l.5,b=2,对两种情形 均有 xG a,b,记 x-2,按照同样的方法依次向下计算，直到求出的相邻两个x的值之差绝对误差小于0.00005,即可求出结果。2.3 Mat lab 程序程序见附录2.2.4 Mat lab程序运行结果:x=1.2500 x=1.3750 x=1.3125x=1.3438x=1.3594x=1.367216X=1.3633x=1.3652x=1.3643x=1.3647x=1.3650 x=1.3651x=1.3652x=1.3652x=1.3652x=1.3652x=1.3652x=l.36523,f(x)=0.000011718三、静电场问题的计算3.1问题描述:设两个同轴矩形金属槽如图3-1所示，外金属槽电位为0,内金属槽电位为100V,求内 电位分布，并绘出电位分布图。3.2原理分析1、3.2.1简单迭代法对某一网格点设一初值，这个初值完全可以任意给定，称为初值电位。虽然，问题的最终结果与初值无关，但若初值选择得当，则计算步骤会 得到简化（当利用计算机来实现迭代计算时，为了简化程序，初值点为 一般可取值为零）。初值电位给定后，按一个固定顺序（点的顺序是从左到右，从下到上）依次计算每点的点位，即利用614+1+119+T,+1,），用围绕它的四个点的电压的平均值作为它的新值，当所有的点计算完 后，用他们的新值代替旧值，即完成了一次迭代计算。然后再进行下一 次迭代计算，直到每一点计算的新值和旧值之差小于指定的范围为止。简单迭代法的特点是用之前一次迭代得到的网络点电位作为下一次迭代 的初值。如在（i,j）点在n+1次迭代时计算公式为：1喊产=喝+1+%-+年陌J+矶U）3.2.2有限差分：二维拉普拉斯方程（羽 y）=甲XX+Ky=于 x,y （1）有限差分法的网格划分，通常采用完全有规律的分布方式，这样可使每个离 散点上得到相同形式的差分方程，有效的提高解题速度，经常采用的是正方形网 格划分。设网格节点（i,j）的电位为，其上下左右四个节点的电位分别为在h充分小 的情况下，可以为基点进行泰勒级数展开：9i=Rj+如 h+一双 h+一双 h+L1J dy dy dy黑吟急-M+9g-dy20(pi+i.=/j+-h+-h2 j J dx 2 a%21 d3cp6 dx3川+1 d3(p6 dx3川+把以上四式相加，在相加的过程中，h的所有奇次方项都抵消了。得到的结 果的精度为h的二次项。Q,j+i+%-1+0-1,/+9i+i,j=40,/+*+（2）由于场中任意点都满足泊松方程：一件落皿）式中为场源，则式（2）可变为:1 h2外，=4（%加+九t+Of）一丁方（羽丁）（3）对于无源场，则二维拉普拉斯方程的有限差分形式为:1 z、5i,j=+0f+5i+i,j）（4）上式表示任一点的电位等于围绕它的四个等间距点的电位的平均值，距离h 越小则结果越精确，用式（4）可以近似的求解二维拉普拉斯方程。边界条件：9（犬,y|内槽二100；(x,y)|外槽二03.2.3解题过程解：在直角坐标系中，金属槽中的电位函数小满足拉普拉斯方程:d2（p d2（p-+-=0dx2 dy2其边界条件满足混合型编值问题的边界条件：5(x,y)|久=o=0,dcp前.=仇0(%/)1尸0=0(p(x,y)y=b=100 sin nx 21取步长h=l,X、y方向的网格数为m=16,n=10,共有16*10=160个网孔17*11=187 个节点,其中槽内节点（电位待求点）有15*9=135个，界节点（电位已知点）有187135=52个。设迭代精度为10-63.3Mat I ab程序设计仿真源程序见附录三3.4Matlab仿真结果图3-3运行结果图2223四、热传导方程和波动方程的差分解法4.1 问题描述求热传导方程混合问题：的数值解，其中N、h、k值等参数自取(将计算结果图形化)。4.2原理分析二维热传导方程的初、边值混合问题与一维的相似，在确定差分方程格式并 给出定解条件后，按时间序号分层计算，只是每一层是由二维点阵组成，通常 称为网格。各向同性介质中无热源的二维热传导方程为：，(,，)初始条件是：设时间步长为工，空间步长为h,二维平面xoy分为MXN的网格，并 使，O,则有(),对节点，在时刻(即时刻)有：%,j,kdt Tle-5if f(c)*f(b)l.0e-6)k=k+l;difmax=0.0;for i=2:hy-lfor j=2:hx-l m=meshl(i,j);if(m=2)vold=v(i,j);%计算迭代次数%从2到20行循环%从2到20列循环%取（i,j）点标志值%标志判断%取该点的原值v(i,j)=(l/4)*(v(iT,j)+v(i,j-l)+v(i+l,j)+v(i,j+1);%拉普拉斯 方程差分式dif=v(i,j)-vold;差dif=abs(dif);if(difdifmax)difmax=dif;end end%前后两次迭代值的%取绝对值%所有网格中取最大误值end31end endsubplot(1,2,1),mesh(v)%画三维曲面图axis(-2,hx+3,-2,hy+3,0,100)subplot(1,2,2),contour(v,13)%画等电位线图 hold onx=l:1:hx;y=l:1:hy;xx,yy=meshgrid(x,y);%形成栅格Gx,Gy=gradient(v,0.6,0.6);%计算梯度%quiver(xx,yy,Gx,Gy,-0.5,r)%根据梯度数画箭头axis(1-2,hx+3,-2,hy+3)plot(1,1,hx,hx,1,1,hy,hy,1,1,k)plot(CXI,CXI,CX2,CX2,CXI,CY1,CY2,CY2,CY1,CY1/k)%画外框边线 text(CX1+0.6,CY1+(CY2-CYl)/2,IM00,fontsize,10);%画内框边界线 text(hx/2,hy+1,U=0,fontsize5,10);%外框上边界标注text(hx/2,0,，U=0,J fontsize,10);%外框下边界标注text(-1.7,hy/2,U=O,fontsize5,10);%外框左边界标注text(hx+0.4,hy/2,U=O,fontsize,10);%外框右边界标注hold off32附录4:h=0.1;for k=l:37%从1到37循环37次u(l,k)=0;N=10;a=l/6u(ll,k)=l;for i=2:10%从 2 至(J 10 共 9 步u(i,l)=i*h*i*h;for k=l:36for i-2:10u(i,k+1)=a*u(i+1,k)+(l-2*a)*u(i,k)+a*u(i-l,k);endendendendmesh(u)%画图xlabel(X Axis),ylabel(Y Axis),zlabel(Temperature5),title(ThermalField Distribution5)33