我想在对象坐标中双重积分 2D 加速度数据以获得世界坐标中的 2D 位置。物体总是指向速度方向（假设例如一列火车）。
因此，我尝试将加速度值与速度 verlet积分进行数值积分，将每一步的方向更改为世界坐标中的前一个速度，由速度 verlet 算法提供：
import numpy as np
from math import sqrt
from matplotlib import pyplot as plt
def rotate(a, newXAxis):
r = newXAxis
normX = r / sqrt(np.dot(r.T,r))
normY = [-normX[1], normX[0]]
b = np.dot(np.array([normX, normY]).T, a)
return(b)
"""return true if v > 1 km/h or any speed given"""
def isMoving(deltaXPosition, deltaYPosition, deltaTime, fasterThankmh=1.0):
x = deltaXPosition
y = deltaYPosition
t = deltaTime
if t*t == 0.:
return False
if hasattr(x, "__len__"):
x = x[0]
if hasattr(y, "__len__"):
y = y[0]
if hasattr(t, "__len__"):
t = t[0]
speed = float(fasterThankmh)
return((x*x + y*y) / (t*t) > 0.077160*speed*speed)
def velocity_verlet_integration(Xacc, Yacc,
x0=0., y0=0.,
vx_0=0, vy_0=0,
forward=np.array([1.0, 0.0])):
vx = np.zeros(len(Xacc))
vy = np.zeros(len(Xacc))
x = np.zeros(len(Xacc))
y = np.zeros(len(Xacc))
x[0] = x0
y[0] = y0
vx[0] = vx_0
vy[0] = vy_0
for i in range(len(Xacc)-1):
dt = Xacc[i+1]-Xacc[i]
a = rotate(Yacc[i,:], forward)
x[i+1] = x[i] + vx[i]*dt + 1.0/2.0*a[0]*dt*dt
y[i+1] = y[i] + vy[i]*dt + 1.0/2.0*a[1]*dt*dt
if isMoving(x[i+1]-x[i], y[i+1]-y[i], dt):
forward = np.array([x[i+1]-x[i], y[i+1]-y[i]])
aNext = rotate(Yacc[i+1,:], forward)
vx[i+1] = vx[i] + dt*(a[0] + aNext[0])/2
vy[i+1] = vy[i] + dt*(a[1] + aNext[1])/2
return x, y
用一个简单的圆周运动测试这个：
"""test circle"""
centripetal=-0.2
N = 0.01
xCircle = np.array(range(int(100*10**N)))/float(10**N)
yCircle = np.array([[0.0, centripetal] for i in xCircle])
xvvi, yvvi = velocity_verlet_integration(xCircle, yCircle, 0., 0., 1., 0.)
#plot it
plt.plot(xvvi, yvvi, ".-", label='position with "velocity verlet" integration')
这会导致向外漂移，因为当前方向是基于最后一个速度的，这显然是一个不好的近似值。
谁能指出我更好的解决方案？
一些想法：
- 最佳情况下，需要世界坐标中的最后一个和下一个速度，通过对它们进行平均（例如相加）来做出更好的近似。但在我的方法中，下一个速度取决于世界坐标中的下一个加速度，而后者又需要下一个方向（追逐它自己的尾巴）。
- 如果我使用我的方法来获得下一个速度的第一个近似值，从而获得下一个方向，我可以通过上面的想法使用它来改进当前方向。现在我可以对下一个速度和下一个方向做出更好的近似，并再次使用它来改进当前方向。这可能是可能的解决方案，尽管它看起来真的很难看。
