python - 静电力 - 使用 Runge Kutta 模拟测试粒子的轨迹 - 力总是排斥

Question

在 2D 平面的中心，正静电荷 Q 放置在r_prime位置。这种电荷会产生一个静电场E。

现在我想在这个静态电场中放置一个带有电荷 Q 和位置矢量r的测试粒子，并使用 4 阶 Runge-Kutta 方法计算它的轨迹。

对于初始条件

• Q = 1, r_prime=(0,0)
• q = -1, r = (9, 0), v = (0,0)

可以预期，带负电的测试粒子应该向中心的正电荷移动。相反，对于测试粒子 x 分量的时间演化，我得到以下结果：

[9.0,
 9.0,
 8.999876557604697,
 8.99839964155741,
 8.992891977287334,
 8.979313669171093,
 8.95243555913327,
 8.906134626441052,
 8.83385018027209,
 8.729257993736123,
 8.587258805422984,
 8.405449606446608,
 8.186368339303788,
 7.940995661159361,
 7.694260386250479,
 7.493501689700884,
 7.420415546859942,
 7.604287312065716,
 8.226733652039988,
 9.498656905483394,
 11.60015461031076,
 14.621662121713964,
 18.56593806599109,
 ....

第一个迭代步骤的结果显示了正确的行为，但随后粒子奇怪地被排斥到无穷大。我的龙格-库塔方法的实现一定有一个重大缺陷，但我检查了几次，我找不到任何...

有人可以快速查看我的实现，看看他们是否能找到错误。

"""
Computes the trajectory of a test particle with Charge q with position vector r = R[:2] in a 
electrostatic field that is generated by charge Q with fixed position r_prime
"""

import numpy as np
import matplotlib.pyplot as plt

def distance(r, r_prime, n):
    
    """
    computes the euclidean distance to the power n between position x and x_prime
    """
    return np.linalg.norm(r - r_prime)**n


def f(R, r_prime, q, Q):
    
    """ 
    The equations of motion for the particle is given by: 
    
    d^2/dt^2 r(t) = F = constants * q * Q * (r - r_prime)/||r - r_prime||^3   
        
    To apply the Runge-Kutta-Method we transform the above (constants are set to 1) 
    two dimensional second order ODE into four one dimensional ODEs:
        
    d/dt r_x = v_x
    d/dt r_y = v_y
    d/dt v_x = q * Q * (r_x - r_prime_x)/||r - r_prime||^3
    d/dt v_y = q * Q * (r_y - r_prime_y)/||r - r_prime||^3 '''  
    """
    
    r_x, r_y = R[0], R[1]
    v_x, v_y = R[2], R[3]
    
    dist = 1 / distance(np.array([r_x, r_y]), r_prime, 3)

    # Right Hand Side of the 1D Odes
    f1 = v_x
    f2 = v_y
    f3 = q * Q * dist * (r_x - r_prime[0]) 
    f4 = q * Q * dist * (r_y - r_prime[1])
        
    return np.array([f1, f2, f3, f4], float)

# Constants for the Simulation
a = 0.0       # t_0
b = 10.0       # t_end
N = 100    # number of iterations
h = (b-a) / N # time step

tpoints = np.arange(a,b+h,h)

# Create lists to store the computed values
xpoints, ypoints = [], []
vxpoints, vypoints = [], []


# Initial Conditions

Q, r_prime = 1, np.array([0,0], float)  # charge and position of particle that creates the static E-Field

q, R = -1, np.array([9,0,0,0], float)    # charge and its initial position + velocity r=R[:2], v=[2:] 


for dt in tpoints:
    xpoints.append(R[0])
    ypoints.append(R[1])
    vxpoints.append(R[2])
    vypoints.append(R[3])
    
    # Runge-Kutta-4th Order Method
    k1 = dt * f(R, r_prime, q, Q)
    k2 = dt * f(R + 0.5 * k1, r_prime, q, Q)
    k3 = dt * f(R + 0.5 * k2, r_prime, q, Q)
    k4 = dt * f(R + k3, r_prime, q, Q)
    R += (k1 + 2*k2 * 2*k3 + k4)
    
plt.plot(tpoints, xpoints) # should converge to 0

标签： pythonsimulationphysics