Python求解ODE后按时间点取分位数并绘制曲线问题

问题原因 #

你的核心错误是没有将100次仿真的结果按时间点对齐聚合，两次尝试的循环都只遍历单条仿真数据，每次仅基于单次仿真的所有时间点计算分位数，得到的是单个数值，无法和长度为99的时间序列t[1:]匹配，因此触发维度错误。

正确实现逻辑 #

1. 先将存储100条仿真J差分的列表J_diffs转换为二维numpy数组，形状为(100次仿真, 99个时间点)
2. 沿**仿真次数维度（轴0）**计算三个目标分位数，得到的每个分位数结果都是长度为99的数组，和t[1:]一一对应
3. 直接绘制三个分位数曲线即可

修正后的完整代码 #

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import random

N = 1000
I0, R0 = 1, 0
S0 = N - I0 - R0
J0 = I0
beta, gamma = 2/7, 1/7
t = np.linspace(0, 100,100)

empty = []
for i in range(100):
    empty.append(random.uniform(1.5, 2.5)*gamma)

def deriv(y, t, N, beta, gamma):
    S, I, R, J = y
    dS = ((-beta * S * I) / N)
    dI = ((beta * S * I) / N) - (gamma * I)
    dR = (gamma * I)
    dJ = ((beta * S * I) / N)
    return dS, dI, dR, dJ

solns =  []
for empt in empty:
    ces = odeint(deriv, (S0, I0, R0, J0), t, args=(N, empt, gamma))
    solns.append(ces)

J_diffs = []
for sol in solns:
    S, I, R, J = sol.T
    J_diffs.append(np.diff(J))

# 核心改动：转换为二维数组并按时间点计算分位数
J_diffs_arr = np.array(J_diffs)
p2_5 = np.percentile(J_diffs_arr, 2.5, axis=0)
p50 = np.percentile(J_diffs_arr, 50, axis=0)
p97_5 = np.percentile(J_diffs_arr, 97.5, axis=0)

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.set_xlabel('Time in days')
ax.set_ylabel('Number')
ax.grid(visible=True, which='major', c='w', lw=2, ls='-')

# 原始100条曲线调低透明度，突出分位数
for J_diff in J_diffs:
    ax.plot(t[1:], J_diff, 'blue', alpha=0.1, lw=1)

# 绘制分位数曲线
ax.plot(t[1:], p2_5, 'green', label='2.5%分位数', lw=2)
ax.plot(t[1:], p50, 'red', label='50%中位数', lw=2)
ax.plot(t[1:], p97_5, 'green', label='97.5%分位数', lw=2)
# 可选填充上下分位区间，更直观
ax.fill_between(t[1:], p2_5, p97_5, color='green', alpha=0.2)

ax.legend()
plt.show()

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