1

我正在使用四阶 Runge-Kutta 方法求解两个耦合的常微分方程。由于应用此方法，我无法打印 z 的值。源代码如下供参考。请通过打印 z 和 theta 的更新值来帮助我修复此代码。谢谢你。

#Import neeeded modules
import numpy as np
import matplotlib.pyplot as plt

#Input parameters
k = 5 #longitudinal torsional constant
delta = 10**-3 #longitudinal torsional constant
I = 10**-4 #Rotational Inertia
eps = 10**-2 #Coupling constant
m = 0.5


#Time Step
#Setting time array for graph visualization
dt = 0.002 #Time Step
tStop = 0.30 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+dt, dt) #Array of time
z = np.zeros(len(t))
dz = np.zeros(len(t))
theta = np.zeros(len(t))
dtheta = np.zeros(len(t))

#Functions that include the equations of motion
def dYdt(t,u):
    z, dz, theta, dtheta = u
    ddz = -(k*z+0.5*eps*theta)/m
    ddtheta = -(delta*theta+0.5*eps*z)/I
    return np.array([dz, ddz, dtheta, ddtheta])

def rk4(t, u, dt):
    for i in range(len(t)-1):
     # runge_kutta
        k1 = dYdt(t[i], u[i])
        k2 = dYdt(t[i] + dt/2, u[i] + dt/2 * k1)
        k3 = dYdt(t[i] + dt/2, u[i] + dt/2 * k2)
        k4 = dYdt(t[i] + dt, u[i] + dt * k3)
        u.append(u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
        #Unpacking
        z, dz, theta, dtheta = np.asarray(u).T

print(z)

这是我用来提供源代码的运动方程。耦合 ODE

4

1 回答 1

0

这就是我认为你想要的，但我不知道要传递什么参数rk4。另外，u是一维向量，所以我不知道您期望z, dz, theta, dtheta = np.asarray(u).T做什么。结果，此代码不会运行，但它会向您展示一个潜在的设计。

import numpy as np
import matplotlib.pyplot as plt

#Input parameters
k = 5 #longitudinal torsional constant
delta = 10**-3 #longitudinal torsional constant
I = 10**-4 #Rotational Inertia
eps = 10**-2 #Coupling constant
m = 0.5


#Time Step
#Setting time array for graph visualization
dt = 0.002 #Time Step
tStop = 0.30 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+dt, dt) #Array of time

#Functions that include the equations of motion
def dYdt(t,u):
    z, dz, theta, dtheta = u
    ddz = -(k*z+0.5*eps*theta)/m
    ddtheta = -(delta*theta+0.5*eps*z)/I
    return np.array([dz, ddz, dtheta, ddtheta])

def rk4(t, u, dt):
    for i in range(len(t)-1):
     # runge_kutta
        k1 = dYdt(t[i], u[i])
        k2 = dYdt(t[i] + dt/2, u[i] + dt/2 * k1)
        k3 = dYdt(t[i] + dt/2, u[i] + dt/2 * k2)
        k4 = dYdt(t[i] + dt, u[i] + dt * k3)
        u.append(u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
    #Unpacking
    return np.asarray(u).T

z, dz, theta, dtheta = rk4(t, u, dt)
print(z)
于 2021-04-17T03:44:44.697 回答