如何在Python中求解并绘制三个关联二阶ODE的3D曲线

原代码存在几个核心问题：

• 不必要混用SymPy符号运算和数值求解逻辑，手动转换一阶方程组更直接高效
• 未定义关键变量（如dy/dz/dx）且未正确导入SymPy别名sp
• 初始条件维度错误：二阶转一阶后需要6个初始值（x, x', y, y', z, z'），而非1个
• 参数传递格式不符合solve_ivp要求

正确实现步骤 #

我们直接手动将二阶方程组转换为一阶，再用SciPy求解，最后用Matplotlib绘制3D轨迹：

1. 定义一阶状态方程组

设状态变量数组s = [x, x', y, y', z, z']，则一阶导数为：

• s[0]' = s[1]（x的一阶导数就是x'）
• s[1]' = b*s[1] + c*s[3] + d*s[5] + e*s[4] + f*s[2] + g*s[0]（原x''的表达式）
• s[2]' = s[3]（y的一阶导数就是y'）
• s[3]' = q*s[1] + h*s[3] + i*s[5] + p*s[4] + l*s[2] + m*s[0]（原y''的表达式）
• s[4]' = s[5]（z的一阶导数就是z'）
• s[5]' = a*s[1] + w*s[3] + v*s[5] + u*s[4] + o*s[2] + n*s[0]（原z''的表达式）

2. 完整可运行代码

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def ode_system(t, s, params):
    # 解包参数：b,c,d,e,f,g,q,h,i,p,l,m,a,w,v,u,o,n
    b, c, d, e, f, g, q, h, i, p, l, m, a, w, v, u, o, n = params
    # 解包状态变量：x, x', y, y', z, z'
    x, dx, y, dy, z, dz = s
    
    # 计算各状态变量的一阶导数
    dsdt = [
        dx,
        b*dx + c*dy + d*dz + e*z + f*y + g*x,
        dy,
        q*dx + h*dy + i*dz + p*z + l*y + m*x,
        dz,
        a*dx + w*dy + v*dz + u*z + o*y + n*x
    ]
    return dsdt

# 设置参数：按顺序b,c,d,e,f,g,q,h,i,p,l,m,a,w,v,u,o,n
params = [0, 0, 0, 1, 9, 0, 1, 2, 4, 0.5, -1, 0.3, 0.1, -0.2, 0.4, -0.6, 0.8, -0.5]

# 设置初始条件：x(0), x'(0), y(0), y'(0), z(0), z'(0)
initial_conditions = [1.0, 0.0, 0.0, 1.0, 0.0, 0.5]

# 时间区间和采样点
t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 500)

# 求解方程组
sol = solve_ivp(ode_system, t_span, initial_conditions, t_eval=t_eval, args=(params,))

# 提取结果
x_sol = sol.y[0]
y_sol = sol.y[2]
z_sol = sol.y[4]

# 绘制3D曲线
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(x_sol, y_sol, z_sol, label='3D Trajectory')

ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')
ax.set_title('3D Trajectory of the ODE System')
ax.legend()
plt.show()

代码说明 #

• ode_system函数接收时间t、状态变量s和参数组params，返回状态变量的一阶导数数组
• 参数params按原方程的参数顺序传入，确保对应正确
• 初始条件需包含每个变量的初始值和初始导数值，共6个元素
• 求解后提取sol.y中对应x、y、z的行（索引0、2、4），绘制3D轨迹

内容的提问来源于stack exchange，提问作者Mo711