Физика
Закон Ньютона дает вам ОДУ второго порядка u''=F(u) с u=[x,y].Используя v=[x',y'], вы получите систему первого порядка
u' = v
v' = F(u)
, которая является 4-мерной и должна быть решена с использованием 4-мерного состояния.Единственное доступное сокращение - это использование законов Кеплера, которые позволяют привести систему к скалярному порядку на один ОДУ для угла.Но это здесь не задача.
Метод Эйлера
Вы правильно реализовали метод Эйлера для вычисления значений в последнем цикле вашего кода.То, что он может выглядеть нефизическим, может быть связано с тем, что метод Эйлера непрерывно увеличивает орбиту, поскольку он движется за пределы выпуклых траекторий, следующих за касательной.В вашей реализации эта внешняя спираль видна для G=100.
Эффект можно уменьшить, выбрав меньший размер шага, например dt=0.001.
Вы должны выбрать время интегрирования, чтобы быть хорошей частью полной орбиты, чтобы получить презентабельный результат, с указанными выше параметрами вы получите около 2 циклов, что хорошо.
Реализация RK4
Вы допустили несколько ошибок.Каким-то образом вы потеряли скорости, обновления позиций должны быть основаны на скоростях.
Тогда вам следует остановиться на fx(x + .5*kx1, y + .5*kx1, t + .5*dt), чтобы пересмотреть ваш подход, поскольку это несовместимо с любым соглашением об именах.Последовательный, правильный вариант -
fx(x + .5*kx1, y + .5*ky1, t + .5*dt)
, который показывает, что вы не можете отделить интеграцию связанной системы, так как вам нужны обновления y вместе с обновлениями x.Кроме того, значения функции являются ускорениями, поэтому обновляются скорости.Обновления положения используют скорости текущего состояния.Таким образом, шаг должен начинаться с
kx1 = dt * fx(x,y,t) # vx update
mx1 = dt * vx # x update
ky1 = dt * fy(x,y,t) # vy update
my1 = dt * vy # y update
kx2 = dt * fx(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
mx2 = dt * (vx + 0.5*kx1/2)
ky2 = dt * fy(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
my2 = dt * (vy + 0.5*ky1/2)
и т. Д.
Однако, как вы видите, это уже начинает становиться громоздким.Соберите состояние в вектор и используйте векторную функцию для системных уравнений
M, G = 20, 100
def orbitsys(u):
x,y,vx,vy = u
r = np.hypot(x,y)
f = G*M/r**3
return np.array([vx, vy, -f*x, -f*y]);
Затем вы можете использовать поваренную реализацию шага Эйлера или Рунге-Кутты
def Eulerstep(f,u,dt): return u+dt*f(u)
def RK4step(f,u,dt):
k1 = dt*f(u)
k2 = dt*f(u+0.5*k1)
k3 = dt*f(u+0.5*k2)
k4 = dt*f(u+k3)
return u + (k1+2*k2+2*k3+k4)/6
и объединить их в интеграционный цикл
def Eulerintegrate(f, y0, tspan):
y = np.zeros([len(tspan),len(y0)])
y[0,:]=y0
for k in range(1, len(tspan)):
y[k,:] = Eulerstep(f, y[k-1], tspan[k]-tspan[k-1])
return y
def RK4integrate(f, y0, tspan):
y = np.zeros([len(tspan),len(y0)])
y[0,:]=y0
for k in range(1, len(tspan)):
y[k,:] = RK4step(f, y[k-1], tspan[k]-tspan[k-1])
return y
и вызвать их с вашей заданной проблемой
dt = .1
t = np.arange(0,10,dt)
y0 = np.array([10, 0.0, 10, 10])
sol_euler = Eulerintegrate(orbitsys, y0, t)
x,y,vx,vy = sol_euler.T
plt.plot(x,y)
sol_RK4 = RK4integrate(orbitsys, y0, t)
x,y,vx,vy = sol_RK4.T
plt.plot(x,y)