我正在尝试使用更大质量的物体不会移动的理想化来解决围绕更大质量物体运行的物体的位置。我正在尝试使用python中的四阶龙格库塔来解决笛卡尔坐标中的位置。
这是我的代码:
dt = .1
t = np.arange(0,10,dt)
vx = np.zeros(len(t))
vy = np.zeros(len(t))
x = np.zeros(len(t))
y = np.zeros(len(t))
vx[0] = 10 #initial x velocity
vy[0] = 10 #initial y velocity
x[0] = 10 #initial x position
y[0] = 0 #initial y position
M = 20
def fx(x,y,t): #x acceleration
return -G*M*x/((x**2+y**2)**(3/2))
def fy(x,y,t): #y acceleration
return -G*M*y/((x**2+y**2)**(3/2))
def rkx(x,y,t,dt): #runge-kutta for x
kx1 = dt * fx(x,y,t)
mx1 = dt * x
kx2 = dt * fx(x + .5*kx1, y + .5*kx1, t + .5*dt)
mx2 = dt * (x + kx1/2)
kx3 = dt * fx(x + .5*kx2, y + .5*kx2, t + .5*dt)
mx3 = dt * (x + kx2/2)
kx4 = dt * fx(x + kx3, y + x3, t + dt)
mx4 = dt * (x + kx3)
return (kx1 + 2*kx2 + 2*kx3 + kx4)/6
return (mx1 + 2*mx2 + 2*mx3 + mx4)/6
def rky(x,y,t,dt): #runge-kutta for y
ky1 = dt * fy(x,y,t)
my1 = dt * y
ky2 = dt * fy(x + .5*ky1, y + .5*ky1, t + .5*dt)
my2 = dt * (y + ky1/2)
ky3 = dt * fy(x + .5*ky2, y + .5*ky2, t + .5*dt)
my3 = dt * (y + ky2/2)
ky4 = dt * fy(x + ky3, y + ky3, t + dt)
my4 = dt * (y + ky3)
return (ky1 + 2*ky2 + 2*ky3 + ky4)/6
return (my1 + 2*my2 + 2*my3 + my4)/6
for n in range(1,len(t)): #solve using RK4 functions
vx[n] = vx[n-1] + fx(x[n-1],y[n-1],t[n-1])*dt
vy[n] = vy[n-1] + fy(x[n-1],y[n-1],t[n-1])*dt
x[n] = x[n-1] + vx[n-1]*dt
y[n] = y[n-1] + vy[n-1]*dt
最初,无论我以哪种方式调整代码,我的 for 循环都会出错,要么是“'float' 类型的对象没有 len()”(我不明白 float python 可能指的是什么),或“使用序列设置数组元素”(我也不明白它的意思是什么)。我设法摆脱了错误,但我的结果是错误的。我得到 10s 的 vx 和 vy 数组,从 10. 到 109. 的整数 x 数组,以及从 0. 到 99. 的整数数组。
我怀疑 fx(x,y,t) 和 fy(x,y,t) 或者我将 runge-kutta 函数编码为与 fx 和 fy 一起使用的方式存在问题,因为我使用了相同的 runge -kutta 代码用于其他功能,它工作正常。
我非常感谢任何帮助找出我的代码为什么不起作用的帮助。谢谢你。