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

标准文档实 验 报 告课程名称： 计算物理实践 专 业： 应用物理学 学 号： 姓 名： 完成日期： 年 月 目 录一、简单物理实验的模拟及实验数据处理1问题描述：1单摆运动原理1模型的建立1流程图21.5 Matlab程序设计仿真2程序3单摆演示截图3二、方程组的数值解法4问题描述：4原理分析4二分法理论42.2.2 分析求解42.3 Matlab程序52.4 Matlab程序运行结果：5三、静电场问题的计算7问题描述：7原理分析71简单迭代法7有限差分：83.2.3 解题过程9程序设计仿真9仿真结果10四 、 热传导方程和波动方程的差分解法11问题描述11原理分析11具体步骤134.4 MATLAB程序设计仿真13程序运行结果13结束语15参考文献16附录1：17附录2：18附录3：19附录4：21实用文案标准文档一、简单物理实验的模拟及实验数据处理：编写单摆运动演示程序。在不考虑空气阻力和很小的假设下，单位质量小球做理想简谐运动，此时。取g=9.8，L=1, =pi/4.设在某一时刻，单摆的摆线偏离垂直线的角位移为，将重力mg 分解为径向力 F 和切向力T ，则T的大小为mg sin,切向加速度为=L.根据牛顿第二定律得方程ma=m=mg从而单摆运动的微分方程为 =由于=+.当很小时， 所以单摆的微分方程可表示为= 上式表明，当很小时，单摆的角加速度与角位移成正比，但方向相反，且方程的解可表示为=A 建立物理模型，假设单摆运动过程中的小球中心点的坐标为(x, y)。根据几何关系x, y满足如下关系：在不考虑空气阻力和很小的假设下= 与式比较知 A= 所以，单摆的运动方程为将=pi/4,g=9.8,L=1带入得，开 始输入q5= t=1?q YY=-L*;X=L*;= N结束t=t+dt图1.1 程序流程图 Matlab程序设计仿真通过set函数将变现和小球的图像句柄，加入X、Y的动态参量，它们得轨迹变化即为上面分析的轨迹方程。用line函数画出初始的位置，并将句柄分别给sphere、lp，然后每隔dt时间刷新一次。尝试对dt的设置发现，值在的效果比较好。Matlab程序程序见附录1。图1.2 单摆演示图二、方程组的数值解法2.1问题描述：二分法求解方程x3+4x2100在区间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()f()<0,则x*,记，,否则x*,记，,对两种情形均有x*,记第k步，计算f()f()<0,则x*,记，,否则x*,对这两种情形均有x* ,记 ,k=1,2,.x*且所以数列满足x*- |()=)=.=)=即从而当k充分大，且可由x*- 控制精度。2.2.2 分析求解令f(x)=-x3+4x210,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,此时记a=1,b=1.5,否则 x,记a=1.5 ,b=2. 对两种情形均有x,记x=,按照同样的方法依次向下计算，直到求出的相邻两个x的值之差绝对误差小于0.00005，即可求出结果。2.3 Matlab程序程序见附录2.2.4 Matlab程序运行结果：三、静电场问题的计算：设两个同轴矩形金属槽如图3-1所示，外金属槽电位为0，内金属槽电位为100V，求内电位分布，并绘出电位分布图。 图3-13.21简单迭代法1、 对某一网格点设一初值，这个初值完全可以任意给定，称为初值电位。虽然，问题的最终结果与初值无关，但若初值选择得当，则计算步骤会得到简化（当利用计算机来实现迭代计算时，为了简化程序，初值点为一般可取值为零）。2、 初值电位给定后，按一个固定顺序（点的顺序是从左到右，从下到上）依次计算每点的点位，即利用,用围绕它的四个点的电压的平均值作为它的新值，当所有的点计算完后，用他们的新值代替旧值，即完成了一次迭代计算。然后再进行下一次迭代计算，直到每一点计算的新值和旧值之差小于指定的范围为止。简单迭代法的特点是用之前一次迭代得到的网络点电位作为下一次迭代的初值。如在（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个，界节点（电位已知点）有187135=52个。设迭代精度为源程序见附录三图3-3 运行结果图四 、 热传导方程和波动方程的差分解法1.问题描述求热传导方程混合问题： 的数值解，其中N、k值等参数自取（将计算结果图形化）。二维热传导方程的初、边值混合问题与一维的相似，在确定差分方程格式并给出定解条件后，按时间序号分层计算，只是每一层是由二维点阵组成，通常称为网格。各向同性介质中无热源的二维热传导方程为：，(,).初始条件是： 设时间步长为，空间步长为 h ，二维平面 xoy 分为MN的网格，并使，（），则有（），对节点，在时刻（即时刻）有： 将差分格式代入偏微分方程中，可得： （）就所设定的具体问题来讨论边界条件:在边界的和区域以及整个，边界均为绝热壁；而在边界的区域为与恒温热源相连的口。和两边界温度始终为0，实际上也是与恒温源相连的。也就是说，对于绝热壁应满足：上述边界条件的差分近似式为： 即： 对于与恒温源相连的边界，在热传导过程中始终有恒定的热流，常可取归一化值，例如高温热源可取“1”，而低温热源可取“0”。对于这种情况边界条件还有：综合上述初值、边值混合问题，并设初始时刻各点温度均为零，则上述差分格式可归纳为： 可以证明，对于二维热传导方程，若满足，则差分格式式或式就是稳定的差分格式，一般的讲，对于n维抛物线型微分方程差分格式稳定的充分条件是：。、和，取h=0.1，2.计算初值和边值：；3 用差分格式计算；4.4 MATLAB程序设计仿真源程序见附录四4.5 MATLAB程序运行结果图4-1 运行结果图结束语计算物理学是利用电子计算机进行数据采集、数值计算和数字仿真来发现和研究物理现象与物理规律的一门现代交叉学科。经过本次计算物理学实验周的学习，我学会了怎样利用计算机来进行数据分析和用Matlab来进行数学以及物理问题的仿真而得到问题的结果。我认识到自己对于以前学习过的一些课程掌握得还不够，Matlab编程语言的运用也不够熟练。通过这次实验也很好的巩固了以前学习的一些知识点，将理论运用于实践。这次实验让我认识到数理方程的实用性，掌握了利用差分代替微分来求解波动方程、热传导方程、拉普拉斯方程等的基本原理和方法。本次实践涉及到的二维拉普拉斯方程以及二维热传导方程的解题方法，都是先将连续的方程以及边界条件离散化，再用计算机进行计算，因为计算机智能对离散的数值进行计算。对于非线性方程的求解往往是采用迭代的方法求解，本次实践主要涉及了Newton迭代法的重要思想，也是将连续的方程离散化后再进行计算。通过本次实践，我对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;5;while 1 t=t+dt theta=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)x3+4*x2-10;c=(a+b)/2;while abs(b-a)>1e-5if f(c)*f(b)<0a=c;elseb=c;endc=(a+b)/2;x=cendfprintf('nx=%.5f,f(x)=%.5fn',x,f(x);附录3：FI=100;hx=23;hy=23;%定外框网络；x轴123,123CX1=9;CX2=15;CY1=9;CY2=15;%定内框位置：x轴915，y轴915DX=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;for i=2:hy-1 %从2到20行循环for j=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; end %所有网格中取最大误值endend endendsubplot(1,2,1),mesh(v)%画三维曲面图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+0.4,hy/2,'U=0','fontsize',10);%外框右边界标注hold off附录4：h=0.1;for k=1:37 %从1到37循环37次u(1,k)=0;N=10;a=1/6 u(11,k)=1;for i=2:10 %从2到10共9步 u(i,1)=i*h*i*h; for k=1:36 for i=2:10 u(i,k+1)=a*u(i+1,k)+(1-2*a)*u(i,k)+a*u(i-1,k); end end end endmesh(u)%画图xlabel('X Axis'),ylabel('Y Axis'),zlabel('Temperature'),title('Thermal Field Distribution')实用文案