4-й порядок Рунге-Кутты для решения системы ODE 2-го порядка с использованием Python

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

И у меня очень странный график решений ... У меня есть: введите описание изображения здесь

Правильный график: введите описание изображения здесь

Я не могу найти проблем в своей рунге-кутте. Помогите, пожалуйста.

Мой код здесь:

dt = 0.04

#initial conditions
t.append(0)
zdot.append(0)
z.append(A)
thetadot.append(0)
theta.append(B)

#derrive functions
def zdotdot(z_cur, theta_cur):
   return -omega_z * z_cur - epsilon / 2 / m * theta_cur
def thetadotdot(z_cur, theta_cur):
   return -omega_theta * theta_cur - epsilon / 2 / I * z_cur 
i = 0
while True:
    # runge_kutta
    k1_zdot = zdotdot(z[i], theta[i])
    k1_thetadot = thetadotdot(z[i], theta[i])

    k2_zdot = zdotdot(z[i] + dt/2 * k1_zdot, theta[i])
    k2_thetadot = thetadotdot(z[i], theta[i]  + dt/2 * k1_thetadot)

    k3_zdot = zdotdot(z[i] + dt/2 * k2_zdot, theta[i])
    k3_thetadot = thetadotdot(z[i], theta[i]  + dt/2 * k2_thetadot)

    k4_zdot = zdotdot(z[i] + dt * k3_zdot, theta[i])
    k4_thetadot = thetadotdot(z[i], theta[i]  + dt * k3_thetadot)

    zdot.append (zdot[i] + (k1_zdot + 2*k2_zdot + 2*k3_zdot + k4_zdot) * dt / 6)
    thetadot.append (thetadot[i] + (k1_thetadot + 2*k2_thetadot + 2*k3_thetadot + k4_thetadot) * dt / 6)

    z.append (z[i] + zdot[i + 1] * dt)
    theta.append (theta[i] + thetadot[i + 1] * dt)
    i += 1

python numerical-methods ode runge-kutta
person Nikolay    schedule 05.05.2019    source источник
comment
В вашем вопросе нет вопросов ;-) Пожалуйста, укажите, что вы хотите знать.   -  person Kaadzia    schedule 05.05.2019

Ответы (1)


arrow_upward
0
arrow_downward

Ваше состояние состоит из 4 компонентов, поэтому вам нужно по 4 уклона на каждом этапе. Должно быть очевидно, что наклон / обновление для z не может исходить из k1_zdot, он должен быть k1_z, который должен быть вычислен ранее как

k1_z = zdot

и на следующем этапе

k2_z = zdot + dt/2*k1_zdot

и т.п.

Но лучше использовать векторизованный интерфейс.

def derivs(t,u):
    z, theta, dz, dtheta = u
    ddz = -omega_z * z - epsilon / 2 / m * theta
    ddtheta = -omega_theta * theta - epsilon / 2 / I * z 
    return np.array([dz, dtheta, ddz, ddtheta]);

а затем использовать стандартные формулы для RK4

i = 0
while True:
    # runge_kutta
    k1 = derivs(t[i], u[i])
    k2 = derivs(t[i] + dt/2, u[i] + dt/2 * k1)
    k3 = derivs(t[i] + dt/2, u[i] + dt/2 * k2)
    k4 = derivs(t[i] + dt, u[i] + dt * k3)

    u.append (u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
    i += 1

а затем распаковать как

z, theta, dz, dtheta = np.asarray(u).T
person Lutz Lehmann    schedule 05.05.2019