计算物理基础-曲线作图、动画
在实际问题中，很多变量随时间、空间变化。和公式相比，曲线作图、动画在展示上更直观清晰，也有助于我们对物理过程的理解。这里举一维、二维的案例来说明。
（1） 正弦函数和泰勒展开近似
从展开式看，正弦函数的不同阶近似需要计算幂函数和阶乘。在Octave中，计算阶乘可用函数prod(1:N)，对应N! x是1维数组时，计算乘法和幂次需要用., 比如x.^5是对应x的每个分量都是5次方。对应程序中前7行是直接用plot函数画曲线并标记x、y轴，第8行hold on是保持不同曲线在同一张图，第9行的y2是1阶近似，第10行到第13行是根据求和公式不断增加项并画图，最后一行是设置图的x、y方向最小值、最大值。
（2） 曲线平移和旋转
对于圆、椭圆等闭合曲线，我们不能按照显式函数的方法来作图，比较方便的方法是借助参数方程。下面程序的第二、三行就是单位圆的画法。第四行生成椭圆的横、纵坐标，并保存在矩阵XY中，通过矩阵加法将中心移到指定的[x_0;y_0]，或者根据 R 的矩阵乘法实现旋转。
在循环中加入data(:,i)=getframe; hold off 语句可以保存每次的图，并逐帧连续播放得到动画，对应的语句是 movie(data,1)。具体效果参考视频讲解。
（3） 二维显示
在波的干涉模拟中，我们要注意image函数的用法，可以将0到60的数据用不同颜色显示。这里y_1+y_2的最大值、最小值分别为2、-2，我们用了image((z+2)*15)。phase1和phase2对应平面各点到两个波源的相位差，数值为负时，说明振动没有传播到，就用phase1>0实现类似阶跃函数的效果。动图效果参考视频。norm([a b])返回是a、b的平方和开根号，对应的是计算距离。
程序附录：
A. 正弦函数和多项式近似
clear
x=linspace(0,2*pi);
y=sin(x);
plot(x,y,'.')
xlabel('x')
ylabel('y')
legend('y=sin(x)')
hold on
y2=x;
for i=1:7
plot(x,y2)
y2=y2+(-1)^i*x.^(2*i+1)/prod(1:2*i+1);
end
axis([0 2*pi -1 1])
B. 振动合成和莉萨如图形
clear
t=linspace(0,10,300);
w_1=3;w_2=3;
x=cos(w_1*t);
subplot(2,3,1)
dphi=pi/4;
y=cos(w_2*t+dphi);
plot(x,y)
axis equal
subplot(2,3,2)
dphi=pi/2;
y=cos(w_2*t+dphi);
plot(x,y)
axis equal
subplot(2,3,3)
dphi=3*pi/4;
y=cos(w_2*t+dphi);
plot(x,y)
axis equal
C. 圆和椭圆
clear
t=linspace(0,2*pi,101);
a=cos(t)/5;b=sin(t)/3;
for i=1:length(t)
r=[cos(t(i)) -sin(t(i)); sin(t(i)) cos(t(i))];
XY=[a;b];
XY=r*XY+[cos(t(i)); sin(t(i))];
plot(XY(1,:),XY(2,:),'b')
hold on
plot(cos(t),sin(t),'r')
axis([-1.5 1.5 -1.5 1.5])
axis equal
data(:,i)=getframe;
hold off
end
movie(data,1)
D. 波的传播
clear
N=101;c=(N+1)/2;z=zeros(N,N);t=linspace(0,2);
for k=1:length(t)
for i=1:N
for j=1:N
phase=15*t(k)-norm([i-c j-c])/2;
z(i,j)=cos(phase)*(phase>0);
end
end
image((z+1)*30); colormap(jet)
axis equal; data(:,k)=getframe;
end
movie(data,1)
E. 波的干涉
clear
N=101;c=(N+1)/2;z=zeros(N,2*N);t=linspace(0,3);
for k=1:length(t)
for i=1:N
for j=1:2*N
phase1=15*t(k)-norm([i-c j-1.75*c])/2;
phase2=15*t(k)-norm([i-c j-2.25*c])/2;
z(i,j)=cos(phase1)*(phase1>0)+cos(phase2)*(phase2>0);
end
end
image((z+2)*15); colormap(jet);
axis equal; data(:,k)=getframe;
end
movie(data,1)
