Verlet 算法的简单示例（Matlab）

贡献者： addis; ACertainUser

本文需要更多参考文献。

　　 Verlet 算法常用于数值计算一些经典动力学问题。一些论坛探讨可以参考 https://www.zhihu.com/question/22531466/answer/27695021 与 https://www.zhihu.com/question/22531466/answer/763079874。

　　简而言之，Verlet 算法认为，某时刻粒子的位置可由前两时刻粒子的位置与受力递推得到。假定时间步长为 $Δ t$ .

\begin{matrix} (1) & \begin{aligned} r (t) & = 2 r (t - Δ t) - r (t - 2 Δ t) + \frac{f (t - Δ t)}{m} (Δ t)^{2} t > 0, \\ r (0) & = r_{0}, \\ r (- Δ t) & = r_{0} - v_{0} Δ t + \frac{f (0)}{2 m} (Δ t)^{2} . \end{aligned} \end{matrix}

　　其中， $r_{0}, v_{0}$ 分别为初始时刻 $t = 0$ 粒子的位置与速度，是需要给定的已知量； $f = f (t)$ 是 $t$ 时刻粒子的受力。

　　以下是使用 Verlet 算法模拟单粒子圆周运动的 Matlab 代码。

clc
clear

function f = centriforce(r) %有心力，例如引力、点电荷的静电力
  mag = 1/(norm(r)^2);
  f = -mag*r/norm(r);
end

r0=[10;0]; %点粒子的初始位置
v0=[0;2]; %点粒子的初始速度。v0=[0;sqrt(10)]时为匀速率圆周运动
dt = 0.01;

R0 = r0 - v0*dt + 100*centriforce(r0)/2*dt^2;
R1 = r0;
R2 = [0;0];
i=0;

while 1
  R2 = 2*R1 - R0 + 100*centriforce(R1)*dt^2;
  if mod(i,5) == 0
    clf
    hold on
    scatter(0,0);
    scatter(R2(1),R2(2));
    axis equal
    axis([-10 10 -10 10])
    hold off
    pause(0.05);
    i = 0;
  endif

  R0 = R1;
  R1 = R2;
  R2 = [0;0];
  i=i+1;
end

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。