本文试图用龙格库塔四阶法数值求解两个常微分方程组。 初始系统: 要解决的系统:
我在我的龙格库塔找不到麻烦。请帮帮我。在
我的代码在这里:
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
您的状态有4个组件,因此每个阶段需要4个斜坡。应该很明显,},这是之前要计算的
z
的斜率/更新不能来自k1_zdot
,它必须是{在下一阶段
^{pr2}$等等
但最好是使用矢量化接口
然后使用RK4的标准公式
然后将其解包为
相关问题 更多 >
编程相关推荐