Verlet에서 Python으로 입자가 누설 됨

Question

두 개의 몸체에 대한 운동 방정식을 해결할 수있는 휴 용 코드가 있지만 그 결과 입자가 돌아가고 어디에서 오류가 발생하는지 찾을 수 없습니다.

import numpy as np 
import matplotlib.pyplot as plt 

plt.style.use('ggplot') 

DIM = 2 
N = 2 
ITER = 1000 

def acc(r, v, a): 
    for i in range(N - 1): 
     for j in range(i+1, N): 
      r2 = 0.0 
      rij = [0.0, 0.0] 
      for k in range(DIM): 
       rij[k] = r[i][k] - r[j][k] 
       r2 += rij[k] * rij[k] 
      fg = -1.0 /np.power(np.sqrt(r2), 3) 
      for k in range(DIM): 
       a[i][k] += fg * rij[k] 
       a[j][k] -= fg * rij[k] 
    return a 

def verlet(r, v, a, dt): 
    for i in range(N): 
     for k in range(DIM): 
      v[i][k] += 0.5 * a[i][k] * dt 
      r[i][k] += v[i][k] * dt 
    a = acc(r, v, a) 
    for i in range(N): 
     for k in range(DIM): 
      v[i][k] += 0.5 * a[i][k] * dt 
    return [r,v] 


r = np.zeros((N, DIM)) 
v = np.zeros((N ,DIM)) 
a = np.zeros((N, DIM)) 

r[0] = [0.5,0.0] 
v[0] = [0.0,1.0] 

r[1] = [-0.5,0.0] 
v[1] = [0.0,-1.0] 

dt = 0.01 
plt.ion() 
for i in range(ITER): 
    r,v = verlet(r, v, a, dt) 
    plt.scatter(r[0][0], r[0][1]) 
    plt.scatter(r[1][0], r[1][1],color='yellow') 
    plt.pause(0.00005) 

그리고 시간이 지남에 축적되지 않습니다 velocity Verlet

Answer 1

가속에 설명 된 알고리즘을 사용하는 방식의 속도와 위치 할 : 그것은 프로세스의 각 단계에서 처음부터 계산해야한다. 그래서,

1. 는 acc의 시작 부분에 제로에 의해 모두 acc 및
2. verlet 초기화 a의 인수 목록에서 a를 제거합니다.
3. a = acc(r, v) 호를 verlet의 처음으로 이동하십시오.

예상대로 서로에 대해 회전하는 기사가 표시됩니다. 문제에 관련없는하지만 NumPy와의 효과적인 사용을위한 중요한

가 : 벡터를 추가하고 뺄 NumPy와 그것을두고, 벡터 좌표를 통해 루프를 제거. 나는이 방법으로 acc 방법을 다시 썼다 :

def acc(r, v): 
    a = np.zeros((N, DIM)) 
    for i in range(N - 1): 
     for j in range(i+1, N): 
      rij = r[i, :] - r[j, :] 
      r2 = np.dot(rij, rij) 
      fg = -1.0 /np.power(np.sqrt(r2), 3) 
      a[i, :] += fg * rij 
      a[j, :] -= fg * rij 
    return a