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

....南京邮电大学实验报告课程名称：计算物理实践专业：应用物理学学号：姓名：完成日期：年月

目录一、简单物理实验的模拟及实验数据处理 11.1问题描述： 11.2单摆运动原理 11.3模型的建立 11.4流程图 21.5Matlab程序设计仿真 21.6Matlab程序 31.7单摆演示截图 3二、方程组的数值解法 42.1问题描述： 42.2原理分析 42.2.1二分法理论 42.2.2分析求解 42.3Matlab程序 52.4Matlab程序运行结果： 5三、静电场问题的计算 73.1问题描述： 73.2原理分析 73.2．1简单迭代法 73.2.2有限差分： 83.2.3解题过程 93.3Matlab程序设计仿真 93.4Matlab仿真结果 10四、热传导方程和波动方程的差分解法 114.1问题描述 114.2原理分析 114.3具体步骤 134.4MATLAB程序设计仿真 134.5MATLAB程序运行结果 13结束语 15参考文献 16附录1： 17附录2： 18附录3： 19附录4： 21..一、简单物理实验的模拟及实验数据处理1.1问题描述：编写单摆运动演示程序。在不考虑空气阻力和很小的假设下.单位质量小球做理想简谐运动.此时。取g=9.8.L=1,0=pi/4.1.2单摆运动原理设在某一时刻.单摆的摆线偏离垂直线的角位移为.将重力mg分解为径向力F和切向力T.则T的大小为mgsin,切向加速度为a=Ld2dt2.根据牛顿第二定律得方程ma=mLd从而单摆运动的微分方程为d2dt由于sin=—33!+55当很小时.sin≈所以单摆的微分方程可表示为d2dt2=上式表明.当很小时.单摆的角加速度与角位移成正比.但方向相反.且方程的解可表示为=Acos（ωt+φ）②1.3模型的建立建立物理模型.假设单摆运动过程中的小球中心点的坐标为x,y。根据几何关系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*结束t=t+dt结束t=t+dt图1.1程序流程图1.5Matlab程序设计仿真通过set函数将变现和小球的图像句柄.加入X、Y的动态参量.它们得轨迹变化即为上面分析的轨迹方程。用line函数画出初始的位置.并将句柄分别给sphere、lp.然后每隔dt时间刷新一次。尝试对dt的设置发现.值在0.0005的效果比较好。1.6Matlab程序程序见附录1。1.7单摆演示截图图1.2单摆演示图二、方程组的数值解法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]上连续.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.所以数列｛即xk→x*,从而当k充分大.x*≈xk且可由2.2.2分析求解令f(x)=-x3+4x2－10,f(x)在［1.2］上连续.且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程序运行结果：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.内金属槽电位为100V.求内电位分布.并绘出电位分布图。图3-13.2原理分析3.2．1简单迭代法对某一网格点设一初值.这个初值完全可以任意给定.称为初值电位。虽然.问题的最终结果与初值无关.但若初值选择得当.则计算步骤会得到简化（当利用计算机来实现迭代计算时.为了简化程序.初值点为一般可取值为零）。初值电位给定后.按一个固定顺序（点的顺序是从左到右.从下到上）依次计算每点的点位.即利用φi,j=14(φi,j+1+φi,j-1简单迭代法的特点是用之前一次迭代得到的网络点电位作为下一次迭代的初值。如在（i,j）点在n+1次迭代时计算公式为：φ3.2.2有限差分：二维拉普拉斯方程（1）有限差分法的网格划分.通常采用完全有规律的分布方式.这样可使每个离散点上得到相同形式的差分方程.有效的提高解题速度.经常采用的是正方形网格划分。图3-2迭代法网格划分设网格节点（i,j）的电位为.其上下左右四个节点的电位分别为在h充分小的情况下.可以为基点进行泰勒级数展开：把以上四式相加.在相加的过程中.h的所有奇次方项都抵消了。得到的结果的精度为h的二次项。（2）由于场中任意点都满足泊松方程：式中为场源.则式（2）可变为：（3）对于无源场..则二维拉普拉斯方程的有限差分形式为：（4）上式表示任一点的电位等于围绕它的四个等间距点的电位的平均值.距离h越小则结果越精确.用式（4）可以近似的求解二维拉普拉斯方程。边界条件：|内槽=100；|外槽=03.2.3解题过程解：在直角坐标系中.金属槽中的电位函数φ满足拉普拉斯方程：∂其边界条件满足混合型编值问题的边界条件：φ∂φ取步长h=1,x、y方向的网格数为m=16,n=10,共有16*10=160个网孔17*11=187个节点.其中槽内节点（电位待求点）有15*9=135个.界节点（电位已知点）有187—135=52个。设迭代精度为103.3Matlab程序设计仿真源程序见附录三3.4Matlab仿真结果图3-3运行结果图四、热传导方程和波动方程的差分解法4.1问题描述求热传导方程混合问题：的数值解.其中N、ｈ、k值等参数自取（将计算结果图形化）。4.2原理分析二维热传导方程的初、边值混合问题与一维的相似.在确定差分方程格式并给出定解条件后.按时间序号分层计算.只是每一层是由二维点阵组成.通常称为网格。各向同性介质中无热源的二维热传导方程为：.(,,,).①初始条件是：②设时间步长为τ.空间步长为h.二维平面xoy分为M×N的网格.并使..（）.则有（）.对节点.在时刻（即时刻）有：③将差分格式③代入偏微分方程①中.可得：④（i=1,2,…,N-1,j=1,2,就所设定的具体问题来讨论边界条件:在边界的和区域以及整个.边界均为绝热壁；而在边界的区域为与恒温热源相连的口。和两边界温度始终为0.实际上也是与恒温源相连的。也就是说.对于绝热壁应满足：,（上述边界条件的差分近似式为：即：（）对于与恒温源相连的边界.在热传导过程中始终有恒定的热流.常可取归一化值.例如高温热源可取“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可以证明.对于二维热传导方程.若满足.则差分格式式④或式⑥就是稳定的差分格式.一般的讲.对于n维抛物线型微分方程差分格式稳定的充分条件是：。4.3具体步骤1.给定、、和.取h=0.1.2.计算初值和边值：；；；；；用差分格式计算；4.4MATLAB程序设计仿真源程序见附录四MATLAB程序运行结果图4-1运行结果图结束语计算物理学是利用电子计算机进行数据采集、数值计算和数字仿真来发现和研究物理现象与物理规律的一门现代交叉学科。经过本次计算物理学实验周的学习.我学会了怎样利用计算机来进行数据分析和用Matlab来进行数学以及物理问题的仿真而得到问题的结果。我认识到自己对于以前学习过的一些课程掌握得还不够.Matlab编程语言的运用也不够熟练。通过这次实验也很好的巩固了以前学习的一些知识点.将理论运用于实践。这次实验让我认识到数理方程的实用性.掌握了利用差分代替微分来求解波动方程、热传导方程、拉普拉斯方程等的基本原理和方法。本次实践涉及到的二维拉普拉斯方程以及二维热传导方程的解题方法.都是先将连续的方程以及边界条件离散化.再用计算机进行计算.因为计算机智能对离散的数值进行计算。对于非线性方程的求解往往是采用迭代的方法求解.本次实践主要涉及了Newton迭代法的重要思想.也是将连续的方程离散化后再进行计算。通过本次实践.我对Matlab的使用更加熟练了.解决问题以及整理数据等能力有很大提高.希望我们在学校可以有更多的机会实践.有更多的收获。参考文献[1]陈锺贤.计算物理学.哈尔滨工业大学出版社.2001.3[2]杨振华.郦志新.数学实验.科学出版社.2010.2[3]林亮.吴群英.数值分析方法与实验：基于MATLAB实现.高等教育出版社.2012.9[4]李庆杨.王能超.易大义.数值分析.华中科技大学出版社.2006.7[5]钟季康.鲍鸿吉.大学物理习题计算机解法—MATLAB编程应用.机械工业出版社.2008.1[6]何红雨.电磁场数值计算法与MATLAB实现.华中科技大学出版社附录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;%定外框网络；x轴1~23,1~23CX1=9;CX2=15;CY1=9;CY2=15;%定内框位置：x轴9~15.y轴9~15DX=1+CX2-CX1;DY=1+CY2-CY1;v=ones(hy,hx)*50;%每个网格点的初值mesh1=ones(hy,hx)*2;%每个网格点的计算标志值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;%内框区初值mesh1(CY1:CY2,CX1:CX2)=ones(DX,DY);%内框区标志值k=0;difmax=1.0;%最大误差while(difmax>1.0e-6)k=k+1;%计算迭代次数difmax=0.0;fori=2:hy-1%从2到20行循环forj=2:hx-1%从2到20列循环m=mesh1(i,j);%取（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(dif>difmax)difmax=dif;e