二自由度摩托车运动微分方程python

下面是用 Python 实现的二自由度摩托车运动微分方程：

import numpy as np
from scipy.integrate import odeint

def motorbike_equations(state, t, m1, m2, l1, l2, g):
    theta1, theta2, dtheta1, dtheta2 = state
    
    # Compute second derivatives of theta1 and theta2
    c = np.cos(theta1 - theta2)
    s = np.sin(theta1 - theta2)
    denominator = m1 * l1**2 + m2 * l2**2 + 2 * m2 * l1 * l2 * c
    d2theta1 = (m2 * l2**2 * (dtheta2**2) * s - 2 * m2 * l1 * l2 * dtheta1 * dtheta2 * s + \
                (m1 + m2) * g * l1 * np.sin(theta1) - m2 * g * l2 * np.sin(theta2)) / denominator
    d2theta2 = ((m1 + m2) * l1**2 * (dtheta1**2) * s + 2 * m2 * l1 * l2 * dtheta1 * dtheta2 * s + \
                (m1 + m2) * g * l1 * np.cos(theta1) * s + (m1 + m2) * g * l2 * np.sin(theta2)) / denominator
    
    # Return state derivatives
    return [dtheta1, dtheta2, d2theta1, d2theta2]

# Initial conditions
theta1_0 = np.pi / 2
theta2_0 = np.pi / 2
dtheta1_0 = 0
dtheta2_0 = 0
state0 = [theta1_0, theta2_0, dtheta1_0, dtheta2_0]

# Parameters
m1 = 100  # mass of body
m2 = 20  # mass of wheel
l1 = 1.5  # distance from body center of mass to front wheel center of mass
l2 = 1  # distance from front wheel center of mass to rear wheel center of mass
g = 9.81  # acceleration due to gravity

# Time points to solve for
t = np.linspace(0, 10, 100)

# Solve the equations
sol = odeint(motorbike_equations, state0, t, args=(m1, m2, l1, l2, g))

# Plot the results
import matplotlib.pyplot as plt
plt.plot(t, sol[:, 0], 'b', label='theta1(t)')
plt.plot(t, sol[:, 1], 'g', label='theta2(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

这段代码中，我们使用了 `odeint` 函数来求解微分方程，给出了初始状态和时间点等参数。最后，我们使用 Matplotlib 库来将结果可视化。