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

南 京 邮 电 大 学实 验 报 告课程名称： 计算物理实践 专 业： 应用物理学 学 号： 姓 名： 完成日期： 年 月 目 录一、简单物理实验的模拟及实验数据处理11.1问题描述：11.2单摆运动原理11.3模型的建立11.4流程图21.5 Matlab程序设计仿真21.6Matlab程序31.7单摆演示截图3二、方程组的数值解法42.1问题描述：42.2原理分析42.2.1二分法理论42.2.2 分析求解42.3 Matlab程序52.4 Matlab程序运行结果：5三、静电场问题的计算73.1问题描述：73.2原理分析73.21简单迭代法73.2.2有限差分：83.2.3 解题过程93.3Matlab程序设计仿真93.4Matlab仿真结果10四 、 热传导方程和波动方程的差分解法114.1问题描述114.2原理分析114.3具体步骤134.4 MATLAB程序设计仿真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的大小为mg sin,切向加速度为a=Ld2dt2.根据牛顿第二定律得方程ma=mLd2dt2=mgsin从而单摆运动的微分方程为 d2dt2=gLsin由于sin=33!+55!.当很小时， sin所以单摆的微分方程可表示为d2dt2=gL 上式表明，当很小时，单摆的角加速度与角位移成正比，但方向相反，且方程的解可表示为=Acos（t+） 1.3模型的建立建立物理模型，假设单摆运动过程中的小球中心点的坐标为(x, y)。根据几何关系x, y满足如下关系：x2+y2=L2tan=yx在不考虑空气阻力和很小的假设下=0costgL 与式比较知 A= 0, =0,=gL所以，单摆的运动方程为x=Lsin=Lsin0costgLy=Lcos=Lcos0costgL将0=pi/4,g=9.8,L=1带入得，x=sin4cost9.8y=-cos4cost9.81.4流程图开 始输入0，g, Lqt=0, dt=0.005= t=1?q YY=-L*cosAngle*cossqrt(g/L*t);X=L*sin(Angle*cossqrt(g/L*t);= N结束t=t+dt图1.1 程序流程图1.5 Matlab程序设计仿真通过set函数将变现和小球的图像句柄，加入X、Y的动态参量，它们得轨迹变化即为上面分析的轨迹方程。用line函数画出初始的位置，并将句柄分别给sphere、lp，然后每隔dt时间刷新一次。尝试对dt的设置发现，值在0.0005的效果比较好。1.6Matlab程序程序见附录1。1.7单摆演示截图图1.2 单摆演示图二、方程组的数值解法2.1问题描述：二分法求解方程x3+4x2100在区间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上存在唯一解，且f(a)f(b)0,记a0=a, b0=b, x0=a0+b02第一步，计算f(a0)f(x0).若f(a0)f(x0)0,则x*(a0, b0),记a1=a0，b1=x0,否则x*(x0, b0),记a1=x0，b1=b0,对两种情形均有x*a0, b0,记x1=a1+b02第k步，计算f(ak-1)f(xk-1)0,则x*ak-1, xk-1,记ak=ak-1，bk=xk-1,否则x*xk-1, bk-1,对这两种情形均有x*xk, bk ,记xk=ak+bk2 ,k=1,2,.k,x*ak, bk且xk=ak+bk2 ，所以数列xk满足x*- xk |12(bk ak)=122(bk-1 ak-1)=.= 12k+1(b0 a0)= 12k+1(b a)0即xkx*,从而当k充分大，x*xk且可由x*- xk|12k+1(b a)控制精度。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) 1e-5if f(c)*f(b)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(difdifmax) 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