Runge Kutta 4 和 python 中的钟摆模拟答案

【问题标题】：Runge Kutta 4 and pendulum simulation in pythonRunge Kutta 4 和 python 中的钟摆模拟
【发布时间】：2019-07-31 14:51:12
【问题描述】：

我正在尝试制作一个 python 程序，该程序使用 runge 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()

我现在更新了代码。我越来越正弦了。谢谢。。

【问题讨论】：

您需要分享您获得的结果与您期望的结果。还有你尝试过什么来解决你的问题
嗨，预期结果是正弦波。我得到了指数输出。

标签： python numpy scientific-computing runge-kutta

【解决方案1】：

您已经应用了 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()

【讨论】：

我按照建议做了，但仍然得到一个非罪波。请查看更新的代码。
您的函数dydt 应该只返回v，正如我在上面所写的。不涉及余弦或类似。而dvdt 取决于y 的值，而不是v。
我现在尝试更新代码。现在我得到正弦波。但是，也必须使用欧拉方法。我想把代码完全放在 RK4 上。如果我们只返回 v，则表明程序理解该函数是基于 sin 函数的。我对 python 很陌生，对不起，我不明白。
你有一个复杂的方法来实现一个order 1方法，你也可以一直使用欧拉。

【解决方案2】：

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()

【讨论】：