使用scipy的solve_ivp求解微分方程后绘图报错，求解决方法

我看你在用solve_ivp求解状态空间形式的微分方程组时，绘图环节出了问题，咱们来一步步搞定它：

先说说你原来代码里的核心问题 #

你最后写的plt.plot(t,X)里的X是你定义的初始全零状态数组，根本不是solve_ivp计算出来的解！而且X是4元素的一维数组，t是10000个点的数组，维度不匹配，这肯定会报错。

正确的姿势：利用solve_ivp的返回结果 #

solve_ivp返回的sol对象里包含了关键的求解结果：

• sol.t：求解过程中自适应选取的时间点数组
• sol.y：对应每个时间点的状态变量，形状是**(状态数, 时间点数)**，你的例子里就是(4, n)
• sol.sol：一个插值函数，可以用来获取任意时间点的状态值（适合自定义时间点的场景）

完整修正后的代码 #

这里给你两种绘图方案，按需选择：

方案1：直接用solve_ivp自带的时间点绘图（简单快捷）

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

vin = 12; vdon = .7;
U = np.array([[vin], [vdon]]);
A1 = np.array([[(-.01-10)/1e-5,-1/1e-5,10/1e-5,0] ,
               [1/1e-8,-1/(.05*1e-8),0,0],
               [10/1e-8,0,(-10-1)/1e-8,-1/1e-8],
               [0,0,1/1e-3,0]]);
B1 = np.array([[1/1e-5,0],
               [0,1/(.05*1e-8)],
               [0,0],
               [0,0]]);

def conv(t, X):
    # 简化计算：矩阵乘法后转成一维数组返回
    xdot = A1.dot(X) + B1.dot(U).flatten()
    return xdot

tspan = [0, .001]
X0 = np.array([0, 0, 0, 0])
sol = solve_ivp(conv, tspan, X0)

# 绘制每个状态变量随时间变化的曲线
plt.figure(figsize=(10,6))
for state_idx in range(sol.y.shape[0]):
    plt.plot(sol.t, sol.y[state_idx, :], label=f'State X{state_idx+1}')
plt.xlabel('Time (s)')
plt.ylabel('State Value')
plt.title('State Variables vs Time')
plt.legend()
plt.grid(True)
plt.show()

方案2：用自定义时间点绘图（需要插值）

如果你一定要用自己生成的t（比如你代码里的10000个点），需要用sol.sol()做插值，因为solve_ivp是自适应步长，不会刚好和你自定义的时间点重合：

# 接上面的代码，在求解得到sol之后：
t_custom = np.linspace(0, .001, 10000)
# 插值获取自定义时间点对应的状态值
X_custom = sol.sol(t_custom)

plt.figure(figsize=(10,6))
for state_idx in range(X_custom.shape[0]):
    plt.plot(t_custom, X_custom[state_idx, :], label=f'State X{state_idx+1}')
plt.xlabel('Time (s)')
plt.ylabel('State Value')
plt.title('State Variables vs Time (Custom Time Points)')
plt.legend()
plt.grid(True)
plt.show()

额外的小优化 #

我还简化了你的conv函数：原来的xdot计算里，A1.dot(X)已经是一维数组（因为X是一维的初始状态），B1.dot(U)是(4,1)的二维数组，用flatten()转成一维后直接相加即可，不需要多余的reshape和asarray操作，代码更简洁。

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