计算物理学实习报告

一、作业实例

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

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

结果图表：