地球轨道-数值方法跳跃

Question

编辑:我现在从地球轨道上得到这些地块。然而，总能量似乎不正确。它不应该在0之间振荡吗？

再次感谢！

​

​

我正在尝试这种数值方法，作业练习b部分)。地球轨道

我的目标是修改a部分的代码，从而返回地球轨道的动能和势能，但我得到了这样的错误：

(8.12).py:31: RuntimeWarning: overflow encountered in multiply
  z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])

我的计算地球轨道动能和势能的程序显然是错误的。我的问题是，我不知道如何修改我的代码，在a)找到势能和动能。

这是我的代码：

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as spi


"""a"""


G = 6.6738*10**-11
M = 1.9891*10**30
h = 3600
y = 1.4710*10**11
vx = 3.0287*10**4

def LeapFrog(f, t_start, t_stop, z0, h):

    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    d = len(z0)  
    z_vec = np.zeros((n, d))
    z_vec[0,:] = z0
    z_half_step=z_vec[0 , :] + (1/2)*h*f(z0,t_vec[0]) 
    
    
    for i in range(0, n - 1):
        z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])
        z_half_step += h*f(z_vec[i+1,:], t_vec[i])   


    return t_vec, z_vec, 


def f(z,t):   #dz/dt
    
    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)
    
    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*x/r**3
    dz[3]=-G*M*y/r**3

    return dz

t_start = 0
t_stop = 24*365*h*5 #5 years
z0 = np.array([0,y,vx,0])
t_vec, z_vec = LeapFrog(f, t_start, t_stop, z0, h)



figure, ax = plt.subplots()

plt.title("Orbit of the Earth")
plt.plot(0,0,'yo', label = 'Sun positon')  
plt.plot(z_vec[:,0],z_vec[:,1], 'g', markersize=1, label='Earth trajectory')
ax.set_aspect('equal')

ax.set(xlim=(-1.55*10**11, 1.55*10**11), ylim = (-1.55*10**11, 1.55*10**11))
a_circle = plt.Circle((0, 0), 1.4710*10**11, fill=False,color='r')
ax.add_artist(a_circle)
ax.set_aspect('equal')

legend = plt.legend(['Sun position','Earths 5-year trajectory','Perfect circle'],loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()



"""

We can see that the trajectory of the earth for 5 years is very slightly non-circular.

"""

对于第二部分，我修改了"dz2“和"dz3”。很明显，这不管用：

"""b"""

    m = 5.9722*10**24
    def f(z,t):   #dz/dt
    
        x=z[0]
        y=z[1] 
        vx=z[2] 
        vy=z[3] 
        r=np.sqrt(x**2+y**2)
    
        dz=np.zeros(4)
    
        dz[0]=vx 
        dz[1]=vy
        dz[2]=-G*M*m/r
        dz[3]=0.5*m*y**2
    
        return dz
    
    t_start = 0
    t_stop = 24*365*h*5 #5 years
    z0 = np.array([0,y,vx,0])
    LeapFrog(f, t_start, t_stop, z0, h)

谢谢

Answer 1

在b部分中不应该有积分)，您使用a部分的结果来计算每个点的能量。

    x,y,vx,vy = z_vec.T
    r = np.hypot(x,y)
    E_kin = 0.5*m*(vx**2+vy**2)
    E_pot = -G*M*m/r
    E = E_kin+E_pot
    plt.plot(t,E)