Рунге Кутта 4 и симуляция маятника в питоне

Question

Я пытаюсь создать программу на Python, которая отображает колебания маятника, используя ступеньку kutta 4. У меня есть уравнение углового увеличения = -(m*g*r/I) * np.sin(y).

Пожалуйста, найдите мой код.Я новичок в Python.

import numpy as np 
import matplotlib.pyplot as plt 
m = 3.0
g = 9.8
r = 2.0
I = 12.0
h = 0.0025 
l=2.0
cycle = 10.0
t = np.arange(0, cycle, h)
n = (int)((cycle)/h)  
initial_angle = 90.0
y=np.zeros(n)
v=np.zeros(n)
def accel(theta):
  return -(m*g*r/I)*np.sin(theta)
y[0] = np.radians(initial_angle) 
v[0] = np.radians(0.0)
for i in range(0, n-1): 
    k1y = h*v[i]
    k1v = h*accel(y[i])
    k2y = h*(v[i]+0.5*k1v)
    k2v = h*accel(y[i]+0.5*k1y)
    k3y = h*(v[i]+0.5*k2v)
    k3v = h*accel(y[i]+0.5*k2y)
    k4y = h*(v[i]+k3v)
    k4v = h*accel(y[i]+k3y)
    y[i+1] = y[i] + (k1y + 2 * k2y + 2 * k3y + k4y) / 6.0 
    v[i+1] = v[i] + (k1v + 2 * k2v + 2 * k3v + k4v) / 6.0

plt.plot(t, y)
plt.title('Pendulum Motion:')
plt.xlabel('time (s)')
plt.ylabel('angle (rad)')
plt.grid(True)
plt.show()

Я обновил код сейчас.Я получаю синусоидальный путь.Спасибо ..

Answer 1

import numpy as np 
import matplotlib.pyplot as plt 
m = 3.0
g = 9.8
r = 2.0
I = 12.0
h = 0.0025 
l=2.0
cycle = 10.0
t = np.arange(0, cycle, h)
# step height h 
n = len(t) 
initial_angle = 90.0
y=np.zeros(n)
v=np.zeros(n)
def accel(theta): return -(m*g*r/I)*np.sin(theta)
y[0] = np.radians(initial_angle) 
v[0] = np.radians(0.0)

for i in range(0, n-1): 
    k1y = h*v[i]
    k1v = h*accel(y[i])

    k2y = h*(v[i]+0.5*k1v)
    k2v = h*accel(y[i]+0.5*k1y)

    k3y = h*(v[i]+0.5*k2v)
    k3v = h*accel(y[i]+0.5*k2y)

    k4y = h*(v[i]+k3v)
    k4v = h*accel(y[i]+k3y)

    # Update next value of y 
    y[i+1] = y[i] + (k1y + 2 * k2y + 2 * k3y + k4y) / 6.0 
    v[i+1] = v[i] + (k1v + 2 * k2v + 2 * k3v + k4v) / 6.0

plt.plot(t, y)
plt.title('Pendulum Motion:')
plt.xlabel('time (s)')
plt.ylabel('angle (rad)')
plt.grid(True)
plt.show()

Answer 2

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

v = dy/dt
acceleration = dv/dt

, чтобы каждый шаг в RK4 имел два компонента

k1y = h*v
k1v = h*accel(y)

k2y = h*(v+0.5*k1v)
k2v = h*accel(y+0.5*k1y)

и т. Д.

---

Полный код дает хорошую синусоидальную волну

import numpy as np 
import matplotlib.pyplot as plt 
m = 3.0
g = 9.8
r = 2.0
I = 12.0
h = 0.0025 
l=2.0
cycle = 10.0
t = np.arange(0, cycle, h)
# step height h 
n = len(t) 
initial_angle = 90.0
y=np.zeros(n)
v=np.zeros(n)
def accel(theta): return -(m*g*r/I)*np.sin(theta)
y[0] = np.radians(initial_angle) 
v[0] = np.radians(0.0)

for i in range(0, n-1): 
    k1y = h*v[i]
    k1v = h*accel(y[i])

    k2y = h*(v[i]+0.5*k1v)
    k2v = h*accel(y[i]+0.5*k1y)

    k3y = h*(v[i]+0.5*k2v)
    k3v = h*accel(y[i]+0.5*k2y)

    k4y = h*(v[i]+k3v)
    k4v = h*accel(y[i]+k3y)

    # Update next value of y 
    y[i+1] = y[i] + (k1y + 2 * k2y + 2 * k3y + k4y) / 6.0 
    v[i+1] = v[i] + (k1v + 2 * k2v + 2 * k3v + k4v) / 6.0

plt.plot(t, y)
plt.title('Pendulum Motion:')
plt.xlabel('time (s)')
plt.ylabel('angle (rad)')
plt.grid(True)
plt.show()