在你这里也贴一下,看看我的方法。
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from matplotlib import style
style.use('ggplot')
fig = plt.figure()
Ax = fig.gca()
epsilon = 0.0001
n = 100000
position = np.zeros((n,2))
velocity = np.zeros((2))
acceleration = np.zeros((2))
position[0,0] = float(input("输入x初始坐标值:"))
position[0,1] = float(input("输入y初始坐标值:"))
velocity[0] = float(input("输入x分量初始速度:"))
velocity[1] = float(input("输入y分量初始速度:"))
r = norm(position[0])
r3 = 1/(r**3)
acceleration = position[0] * -r3
velocity += acceleration * epsilon * 0.5
star = 1
if(norm(velocity) > np.sqrt(2/r)):
star = 0
else:
i = 1
t = 0
while(i<n):
position[i] = position[i-1]+velocity*epsilon
if(position[i-1,1]*position[i,1]< 0):
t += 1
if(t==2):
break
r = norm(position[i])
r3 = 1/(r**3)
acceleration = position[i]*-r3
velocity += acceleration*epsilon
if(norm(velocity) > np.sqrt(2/r)):
star = 0
break
i += 1
if(star == 1):
print('行星的公转周期为:%f' %(i*epsilon))
plt.plot(position[:i,0], position[:i,1],'b-')
plt.plot(0,0,"ro")
plt.show()
else:
print("start is gone")
沒錯,跟我的想法差不多!我是把x、y座標互乘,計算變號次數,要算四次。而脫離速度,我是以速度的點積及 2/r 來比較,自認可以降低浮點數誤差。我們方法不同,但結果相同!