试图模拟一维波

Question

我正在尝试使用一系列谐波振荡器对一维波进行简化模拟.根据教科书的说法x[i],第i个振荡器位置的微分方程(假设均衡x[i]=0)证明了这个:

m*x[i]''=-k(2*x[i]-x[i-1]-x[i+1])

(导数是时间)所以我试图用以下算法数值计算动力学.这osc[i]是第i个振荡器作为具有属性loc(绝对位置),vel(速度),acc(加速度)和bloc(平衡位置)的对象,dt是时间增量:

for every i:
    osc[i].vel+=dt*osc[i].acc;
    osc[i].loc+=dt*osc[i].vel;
    osc[i].vel*=0.99;    
    osc[i].acc=-k*2*(osc[i].loc-osc[i].bloc);

    if(i!=0){
      osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc);
    }
    if(i!=N-1){
      osc[i].acc+=+k*(osc[i+1].loc-osc[i+1].bloc);
    }

您可以在此处查看算法,只需单击即可向第6个振荡器发出脉冲.你可以看到它不仅不会产生波浪,而且还会产生一个总能量增加的运动(即使我加了阻尼!).我究竟做错了什么？

Answer 1

所有给出的答案都提供了大量有用的信息。\n你必须做的是：

\n\n

\n
• 使用 ja72\xc2\xb4s 观察来执行更新 - 您需要保持一致并根据同一迭代批次的位置计算节点处的加速度（即不要混合 $x\xe2\x81\xbdk\xe2\x81 \xbe$ 和 $x\xe2\x81\xbdk+1\xe2\x81\xbe$ 在相同的加速度表达式中，因为这些是属于不同迭代步骤的量）
\n
• 忘记你原来的位置，正如 Tom10 所说 - 你只需要考虑相对位置 - 这个波动方程的行为就像应用于折线的拉普拉斯平滑滤波器 - 使用它的两个直接连接的邻居来松弛一个点，被拉向由这些邻居确定的段的中间
\n
• 最后，如果您希望能量守恒（即不让水面停止振动），那么值得使用辛积分器。辛/半隐式欧拉就可以了。对于这种模拟，您不需要更高阶的积分器，但如果您有时间，请考虑使用 Verlet 或 Ruth-Forest，因为它们都是辛的，并且在精度方面为 2 或 3 阶。龙格-库塔将像隐式欧拉一样，从系统中吸取能量（慢得多），因此您将获得固有阻尼。显式的欧拉最终会给你带来厄运 - 这是一个糟糕的选择 - 总是（特别是对于刚性或无阻尼系统）。还有大量其他集成商，所以请继续尝试。
\n

\n\n

PS，您没有实现真正的隐式欧拉，因为它需要同时影响所有粒子。为此，您必须使用投影共轭梯度法，导出雅可比行列式并求解粒子串的稀疏线性系统。显式或半隐式方法将使您获得动画所期望的“跟随领导者”行为。

\n\n

现在，如果您想要更多，我在 SciLab 中实现并测试了这个答案（懒得用 C++ 对其进行编程）：

\n\n

n=50;\nfor i=1:n\n    x(i) = 1;\nend\n\ndt = 0.02;\n\nk = 0.05;\nx(20) = 1.1;\nxold = x;\nv = 0 * x;\nplot(x, \'r\');\nplot(x*0,\'r\');\nfor iteration=1:10000\n    for i = 1:n\n        if (i>1 & i < n) then\n            acc = k*(0.5*(xold(i-1) + xold(i+1)) - xold(i));\n            v(i) = v(i) + dt * acc;\n            x(i) = xold(i) + v(i) *dt;\n        end\n    end\n    if (iteration/500  == round(iteration / 500) ) then\n        plot(x - iteration/10000,\'g\');\n\n    end\n    xold = x;\nend\nplot(x,\'b\');\n

\n\n

波的演变可以在下面的堆积图中看到：\n

\n