Не удается заставить RK4 определить положение вращающегося тела в Python

Вы нигде не используете функции rkx, rky! В конце определения функции есть два return, которые вы должны использовать return [(kx1 + 2*kx2 + 2*kx3 + kx4)/6, (mx1 + 2*mx2 + 2*mx3 + mx4)/6] (как указано @eapetcho). Также мне непонятна ваша реализация Рунге-Кутты.

У вас есть dv/dt, поэтому вы решаете проблему с v, а затем соответствующим образом обновляете r.

for n in range(1,len(t)): #solve using RK4 functions
    vx[n] = vx[n-1] + rkx(vx[n-1],vy[n-1],t[n-1])*dt
    vy[n] = vy[n-1] + rky(vx[n-1],vy[n-1],t[n-1])*dt
    x[n] = x[n-1] + vx[n-1]*dt
    y[n] = y[n-1] + vy[n-1]*dt

Вот моя версия кода

import numpy as np

#constants
G=1
M=1
h=0.1

#initiating variables
rt = np.arange(0,10,h)
vx = np.zeros(len(rt))
vy = np.zeros(len(rt))
rx = np.zeros(len(rt))
ry = np.zeros(len(rt))

#initial conditions
vx[0] = 10 #initial x velocity
vy[0] = 10 #initial y velocity
rx[0] = 10 #initial x position
ry[0] = 0 #initial y position

def fx(x,y): #x acceleration
     return -G*M*x/((x**2+y**2)**(3/2))

def fy(x,y): #y acceleration
     return -G*M*y/((x**2+y**2)**(3/2))

def rk4(xj, yj):
    k0 = h*fx(xj, yj)
    l0 = h*fx(xj, yj)

    k1 = h*fx(xj + 0.5*k0 , yj + 0.5*l0)
    l1 = h*fy(xj + 0.5*k0 , yj + 0.5*l0)

    k2 = h*fx(xj + 0.5*k1 , yj + 0.5*l1)
    l2 = h*fy(xj + 0.5*k1 , yj + 0.5*l1)

    k3 = h*fx(xj + k2, yj + l2)
    l3 = h*fy(xj + k2, yj + l2)

    xj1 = xj + (1/6)*(k0 + 2*k1 + 2*k2 + k3)
    yj1 = yj + (1/6)*(l0 + 2*l1 + 2*l2 + l3)
    return (xj1, yj1)

for t in range(1,len(rt)):
    nv = rk4(vx[t-1],vy[t-1])
    [vx[t],vy[t]] = nv
    rx[t] = rx[t-1] + vx[t-1]*h
    ry[t] = ry[t-1] + vy[t-1]*h

I suspect there are issues with fx(x,y,t) and fy(x,y,t)

Это так, я только что проверил свой код на fx=3 и fy=y, и я получил хорошую траекторию.

Вот график зависимости ry от rx:

Физика

Закон Ньютона дает ODE u''=F(u) второго порядка с u=[x,y]. Используя v=[x',y'], вы получаете систему первого порядка

u' = v
v' = F(u)

которая является 4-мерной и должна быть решена с использованием 4-х мерного состояния. Единственная доступная редукция - это использование законов Кеплера, которые позволяют уменьшить систему до скалярного порядка одного ОДУ для угла. Но задача не в этом.

Но чтобы получить правильные масштабы, для круговой орбиты радиуса R с угловой скоростью w получается тождество w^2*R^3=G*M, которое подразумевает, что скорость вдоль орбиты равна w*R=sqrt(G*M/R) и периоду T=2*pi*sqrt(R^3/(G*M)). С приведенными данными R ~ 10, w ~ 1, таким образом, G*M ~ 1000 для орбиты, близкой к круговой, поэтому с M=20 для этого потребуется G между 50 и 200 с орбитальным периодом около 2*pi ~ 6. Временной интервал в 10 может соответствовать от половины до 2 или 3 витков.

Метод Эйлера

Вы правильно реализовали метод Эйлера для вычисления значений в последнем цикле вашего кода. То, что это может выглядеть нефизическим, может быть связано с тем, что метод Эйлера непрерывно увеличивает орбиту, поскольку он перемещается за пределы выпуклых траекторий, следующих по касательной. В вашей реализации эту внешнюю спираль можно увидеть для G=100.

Фактически это можно уменьшить, выбрав меньший размер шага, например dt=0.001.

Вы должны выбрать время интегрирования так, чтобы оно составляло хорошую часть полной орбиты, чтобы получить презентабельный результат, с указанными выше параметрами вы получите около 2 петель, что хорошо.

Внедрение РК4

Вы допустили несколько ошибок. Как-то вы потеряли скорости, обновления позиции должны основываться на скоростях.

Тогда вам следовало бы остановиться на 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)
 ky2 = dt * fy(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
 my2 = dt * (vy + 0.5*ky1)

и Т. Д.

Однако, как видите, это уже начинает становиться громоздким. Соберите состояние в вектор и используйте векторную функцию для уравнений системы

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)