Scipy 对一组二阶微分方程进行积分

Question

我有一组二阶微分方程：

我想在 python 中使用 odeint 解决。这可能吗？

我尝试这样做，但我的代码没有给我预期的结果。

import numpy as np, matplotlib.pyplot as plt, matplotlib.font_manager as fm,os
from scipy.integrate import odeint
from mpl_toolkits.mplot3d.axes3d import Axes3D

initial_state = [0.1, 0, 0]
a = 4
time_points = np.linspace(0, 100, 10000)

def my_system(current_state, t):
    theta1, theta2, theta3 = current_state

    d2theta1_dt2 = -a*(2*theta1-theta2-theta3)
    d2theta2_dt2 = -a*(2*theta2-theta1-theta3)
    d2theta3_dt2 = -a*(2*theta3-theta1-theta2)

    return [d2theta1_dt2, d2theta2_dt2, d2theta3_dt2]

xyz = odeint(my_system, initial_state, time_points)

theta1 = xyz[:, 0]
theta2 = xyz[:, 1]
theta3 = xyz[:, 2]

fig = plt.figure(figsize=(12, 9))
ax = fig.gca(projection='3d')
ax.xaxis.set_pane_color((1,1,1,1))
ax.yaxis.set_pane_color((1,1,1,1))
ax.zaxis.set_pane_color((1,1,1,1))
ax.plot(theta1, theta2, theta3, color='g', alpha=0.7, linewidth=0.6)
ax.set_title('plot', fontproperties=title_font)
plt.show()

Answer 1

实际上，插入一阶导数的修复相当快，

def my_system(current_state, t):
    theta1, theta2, theta3, omega1, omega2, omega3 = current_state

    d2theta1_dt2 = -a*(2*theta1-theta2-theta3)
    d2theta2_dt2 = -a*(2*theta2-theta1-theta3)
    d2theta3_dt2 = -a*(2*theta3-theta1-theta2)

    return [ omega1, omega2, omega3, d2theta1_dt2, d2theta2_dt2, d2theta3_dt2]

状态的前三个条目也是角度，后三个条目是它们的速度。

您还需要将 angular 速度的初始值添加到初始状态向量。

Scipy 对一组二阶微分方程进行积分

Integrating a set of second order differential equations with Scipy

python

integrate

scipy

ode