计算物理学实习报告

一、作业实例
1、编写单摆运动演示程序。

在不考虑空气阻力和θ很小的假设下，单位质量
小球做理想简谐运动，此时⎪⎪⎭
⎫ ⎝⎛=L g t cos 0θθ。

在此，取1,8.9,60===L g πθ。

分析过程：
利用一个无穷循环不断的改变移动对象的位置（本问题中是摆锤和摆杆对象），只要时间间隔取得足够小，就不会产生闪烁感，取dt=0.001足够。

加上运动对象都设置“erasemode ”属性值为“xor ”，也不会破坏对象。

所以产生的效果应该是不错的。

以下就是实现这种算法的程序：
t=0; dt=0.001; %dt 是刷新屏幕的时间间隔 while 1 %无穷循环 t=t+dt;
theta=theta0*cos(sqrt(g/l)*t); x=1*sin(theta);
y=-1*cos(theta); %计算摆锤的新位置 set(head,'xdata',x,'ydata',y); %改变摆锤的位置 set(body,'xdata',[0;x],'ydata',[0;y]); %改变摆干的位置 drawnow; %重画屏幕 end 结果图表：
-0.6
-0.4
-0.200.20.4
0.6
-1.2
-1
-0.8
-0.6
-0.4-0.2
水平方向移动距离
竖直方向移动距离
图1 单摆运动快照
2、用牛顿法解方程1-x xe =0，要求绝对误差小于0.00001。

分析过程：
()1-=x xe x f ()x x xe e x f +=' ()()
n n n n x f x f x x '
1
-=+，即x
x x n n xe e xe x x +--=+1
1，然后进行迭代，设迭代精度为510- ，主要算法程序如下：
eps=10^(-5);
k=0; %迭代次数计数置0 while 1
f1=x*exp(x)-1; %计算函数值
f2=exp(x)+x*exp(x); %计算函数的导数值 if f2==0
fprintf('导数为零，错误！') break end
x1=x-f1/f2; %牛顿迭代
k=k+1; %迭代次数增1 err=abs(x1-x); %误差估计 if err<eps
fprintf( '近似根为ｘ=%10.9f,迭代次数为ｋ=%2d\n',x1,k) break end if k==N
fprintf( '迭代次数超限') break end
x=x1 %准备下一轮的迭代初值 end
结果如下：
3、设两个同轴矩形金属槽如图3-3所示，外金属槽电位为0，内金属槽电
位为100V ，求槽内电位分布，并绘出电位分布图。

题3图
分析过程：
再直角坐标系中，同轴矩形金属槽中的电位函数ϕ，满足拉普拉斯方程：
02
222=∂∂+∂∂y x ϕϕ
其边界条件满足第一类边界条件（Dirichlet 边界条件）：
()()0
,100,==外槽内槽y x y x ϕϕ
取步长h=1，x 、y 方向的网格数为m=23，n=23，共有23*23=529个网孔，设迭代精度为610- ，主要算法程序如下：
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 %所有网格中取最大误差值 end end end end
结果图表：
图2 同轴的矩形金属槽内的电位分布图
4、求热传导方程混合问题：⎪⎪⎩⎪
⎪
⎨⎧==-=∂∂=∂∂0),1(,0),0()1(4)0,(22t u t u x x x u x u
t u t x t x ≤≤≤<<<0100,10的数值解，
并与分离变量法求解结果进行图示比较，取N =10，h =0.1, k 值自己根据需要选取。

分析过程： u k i , i=1,2,3,4,5 k=0,1,…,36 由初始条件u(x,0)=4x(1-x) u 0,i =(ih)2=(0.1i)2
分别令i=1,2,3,4,5,则初始条件为 u 0,0=0 u 0,1=0.01 u 0,2=0.04 u 0,3=0.09 u 0,4=0.16 u 0,5=0.25 取α=1/6，则600/1=τ，由式可得
u 1,+k i =1/6 u k i ,1++2/3 u k i ,+1/6 u k i ,1-
由此取k=0,i=1,2,3,4,5, 一直到k=37止。

主要算法程序如下： for k=1:37 %从1到37循环37次 u(1,k)=0; u(11,k)=0;
for i=2:10 %从2到10共9步 u(i,1)=4*i*h*(1-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 end
结果图表：
40
t
热传导三维图
x
u
图3 热传导数值解图
分离变量法求解：求解结果为u(x,t)=x n e k t n n k πππsin 16)1(2
23
30)
1(-∞
=+∑-
程序代码：N=10;
t=1e-5:0.00001:0.005;x=0:0.1:1; w=rcdf(N,t(1));
h=plot(x,w,'linewidth',5,'color','r'); xlabel(' x '); ylabel(' u'); axis([0,1,0,500]) for n=2;length(t) w=rcdf(N,t(n)); set(h,'ydata',w); drawnow; end
function u=rcdf(N,t) x=0:0.1:1;u=0; for k=1:1:2*N
cht=(-1).^(k+1)*16*sin(k*pi*x)/k.^3*pi.^3; u=u+cht*exp(-k^2*pi^2*t); end
结果图表：
0.1
0.2
0.3
0.4
0.50.6
0.7
0.8
0.9
1
050100150200250300350400450
500 热传导二维图
x
u
图4 热传导分离变量法求解图
将两图进行比较，得出的结果很相近。

二、心得小结
通过这次实习，我受益匪浅。

经过这两个星期的训练，我学会了MATLAB 在物理学和数学方面的初步应用，尤其是对理论力学、热学、电磁场等方面有了更加深入的认识，让我们的理论知识得到了具体的应用，同时也加强了对MATLAB 的了解。

在老师和同学们的巨大帮助下，我勉强地完成了任务，我认识到了团队合作的重要性，一个人的力量是渺小的，只有在得到大家的帮助，才能升华自己的力量，我相信，这种团结合作的精神，我会终生受用的。

三、参考文献
[1] 计算物理学，哈尔滨工业大学出版社，程锺贤编著 [2] 数值分析，华中科技大学出版社，李庆杨等编著
[3] 工程电磁场数值计算，高等教育出版社，盛霓虹等编著
[4] MATLAB 及其在理工科课程中的应用，西安电子科技大学出版社，陈怀琛编著 [5] 电磁场数值计算法与MATLAB 实现，华中科技大学出版社，何红雨等编著 [6] MA TLAB 实用教程（第二版），电子工业出版社，苏金明 阮沈勇编著。