Verlet 算法的简单示例（Matlab）

贡献者： addis; ACertainUser

本文需要更多参考文献。

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

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

\begin{equation} \begin{aligned} \boldsymbol{\mathbf{r}} (t) &= 2 \boldsymbol{\mathbf{r}} (t-\Delta t) - \boldsymbol{\mathbf{r}} (t-2\Delta t) + \frac{ \boldsymbol{\mathbf{f}} (t-\Delta t)}{m}(\Delta t)^2 \qquad t>0~,\\ \boldsymbol{\mathbf{r}} (0) &= \boldsymbol{\mathbf{r}} _0~,\\ \boldsymbol{\mathbf{r}} (- \Delta t) &= \boldsymbol{\mathbf{r}} _0 - \boldsymbol{\mathbf{v}} _0 \Delta t + \frac{ \boldsymbol{\mathbf{f}} (0)}{2m}(\Delta t)^2~.\\ \end{aligned} \end{equation}

　　其中，$ \boldsymbol{\mathbf{r}} _0, \boldsymbol{\mathbf{v}} _0$ 分别为初始时刻 $t=0$ 粒子的位置与速度，是需要给定的已知量；$ \boldsymbol{\mathbf{f}} = \boldsymbol{\mathbf{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