from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # # 全局参数与上一轮一致 # alpha_k 0.5 beta_m 0.3 m_crit 5.0 m_eff_fixed 0.0 delta 0.1 epsilon 0.05 zeta_base 0.02 alpha_xi 0.3 beta_xi 0.1 gamma_xi 0.01 def theta(xi_L): return np.exp(-0.5 * xi_L) def phi_t(t, phi_00.1): pulses [ {A: 1.0, t: 41, sigma: 10}, {A: 0.8, t: 72, sigma: 8}, {A: 0.5, t: 99, sigma: 5}, {A: 0.3, t: 152, sigma: 3}, ] phi phi_0 for p in pulses: phi p[A] * np.exp(-((t - p[t])**2) / (2 * p[sigma]**2)) if t 169: phi - 0.7 elif t 166: phi - 0.4 if t 184: phi 0.2 return max(phi, 0.01) def Gamma(t): shocks [ {t: 107, I: 3.0}, {t: 140, I: 2.0}, {t: 184, I: 5.0}, ] g 0.0 for s in shocks: g s[I] * np.exp(-((t - s[t])**2) / 2) return g def compute_sigma_star(sigma_0, kernel_ratio, k_eff, m_eff): suppression 1.0 / (1.0 alpha_k * k_eff beta_m * m_eff) sigma_star sigma_0 * kernel_ratio * suppression return np.clip(sigma_star, 0, sigma_0) def system(t, y, sigma_star, m_eff): Delta_true, xi_L y phi phi_t(t) zeta zeta_base * 0.1 if m_eff m_crit else zeta_base dDelta_dt delta * xi_L - epsilon * Delta_true * phi zeta * Gamma(t) dxi_dt alpha_xi * Delta_true * (1 - theta(xi_L)) - beta_xi * xi_L * phi gamma_xi * xi_L**2 return [dDelta_dt, dxi_dt] # # 灵帝坐标k_eff0.3, kernel_ratio0.58, sigma_01.0, year168 # lingdi_ke 0.3 lingdi_kr 0.58 lingdi_sigma_0 1.0 lingdi_year 168 t_end min(lingdi_year 30, 220) t_span (25, t_end) y0 [0.05, 0.0] t_eval np.linspace(*t_span, 400) # # k1 扫描 # k1_values np.arange(0.50, 0.701, 0.01) results [] for k1 in k1_values: # 原版sigma_star 用 k1 作为 kernel_ratio 的替代即 k_eff0 时的 sigma_star sigma_0 * k1 # 注意这里 k1 实际对应的是 kernel_ratio 参数保持与热力图一致 sigma_star_original compute_sigma_star(lingdi_sigma_0, k1, lingdi_ke, m_eff_fixed) sigma_star_variant compute_sigma_star(lingdi_sigma_0, k1, lingdi_ke, m_eff_fixed) # 原版theta exp(-0.5 * xi_L)但生成项中无观测失真耦合 # 变体在 system 中通过 theta 引入观测失真 # 为区分原版/变体我们定义两个不同的 system def system_original(t, y, sigma_star, m_eff): Delta_true, xi_L y phi phi_t(t) zeta zeta_base * 0.1 if m_eff m_crit else zeta_base # 原版生成项中 theta 不影响 Delta_true 的反馈 dDelta_dt delta * xi_L - epsilon * Delta_true * phi zeta * Gamma(t) dxi_dt alpha_xi * Delta_true * (1 - theta(xi_L)) - beta_xi * xi_L * phi gamma_xi * xi_L**2 return [dDelta_dt, dxi_dt] def system_variant(t, y, sigma_star, m_eff): Delta_true, xi_L y phi phi_t(t) zeta zeta_base * 0.1 if m_eff m_crit else zeta_base # 变体观测失真通过 theta 放大 Delta_true 的生成项 # 这里用 Delta_true * (1 theta(xi_L)) 模拟观测失真对截流差的放大 dDelta_dt delta * xi_L * (1 theta(xi_L)) - epsilon * Delta_true * phi zeta * Gamma(t) dxi_dt alpha_xi * Delta_true * (1 - theta(xi_L)) - beta_xi * xi_L * phi gamma_xi * xi_L**2 return [dDelta_dt, dxi_dt] sol_orig solve_ivp(system_original, t_span, y0, t_evalt_eval, args(sigma_star_original, m_eff_fixed)) sol_vari solve_ivp(system_variant, t_span, y0, t_evalt_eval, args(sigma_star_variant, m_eff_fixed)) xi_orig sol_orig.y[1][-1] xi_vari sol_vari.y[1][-1] delta_xi xi_vari - xi_orig accel_pct (delta_xi / xi_orig * 100) if xi_orig 1e-6 else float(inf) results.append({ k1: k1, xi_orig: xi_orig, xi_vari: xi_vari, delta_xi: delta_xi, accel_pct: accel_pct }) # # 控制台输出 # print(*90) print(k1 扫描结果灵帝thetaexp(-0.5*xi_L)) print(*90) print(f{k1:8} {原版ξ_L:12} {变体ξ_L:12} {Δξ_L:12} {加速效应%:12}) print(-*90) for r in results: accel_str f{r[accel_pct]:.2f} if r[accel_pct] ! float(inf) else inf print(f{r[k1]:8.2f} {r[xi_orig]:12.4f} {r[xi_vari]:12.4f} {r[delta_xi]:12.4f} {accel_str:12}) print(*90) # 找临界点 k1_crit_orig None k1_crit_vari None for i, r in enumerate(results): if r[xi_orig] 10 and k1_crit_orig is None: # 线性插值找精确临界点 if i 0: r_prev results[i-1] k1_crit_orig r_prev[k1] (10 - r_prev[xi_orig]) / (r[xi_orig] - r_prev[xi_orig]) * (r[k1] - r_prev[k1]) else: k1_crit_orig r[k1] if r[xi_vari] 10 and k1_crit_vari is None: if i 0: r_prev results[i-1] k1_crit_vari r_prev[k1] (10 - r_prev[xi_vari]) / (r[xi_vari] - r_prev[xi_vari]) * (r[k1] - r_prev[k1]) else: k1_crit_vari r[k1] print(f 【临界点分析】) print(f原版临界 k1^原版 {k1_crit_orig:.4f}原版ξ_L10) print(f变体临界 k1^变体 {k1_crit_vari:.4f}变体ξ_L10) if k1_crit_orig and k1_crit_vari: delta_k1 k1_crit_orig - k1_crit_vari print(fΔk1 k1^原版 - k1^变体 {delta_k1:.4f}) print(f观测失真降低相变阈值 {delta_k1:.4f}相对降幅 {delta_k1/k1_crit_orig*100:.2f}%) print(*90) # # 绘图k1 - ξ_L 曲线 # fig, ax plt.subplots(figsize(12, 7)) k1_arr [r[k1] for r in results] xi_orig_arr [r[xi_orig] for r in results] xi_vari_arr [r[xi_vari] for r in results] ax.plot(k1_arr, xi_orig_arr, b-o, linewidth2, markersize5, label原版 ξ_L) ax.plot(k1_arr, xi_vari_arr, r-s, linewidth2, markersize5, label变体 ξ_L) ax.axhline(y10, colorgreen, linestyle--, linewidth2, label临界阈值 ξ_L10) if k1_crit_orig: ax.axvline(xk1_crit_orig, colorblue, linestyle:, linewidth1.5, alpha0.7) ax.annotate(fk1^原版{k1_crit_orig:.3f}, xy(k1_crit_orig, 10), xytext(k1_crit_orig0.005, 12), fontsize10, colorblue, bboxdict(boxstyleround,pad0.3, facecolorlightblue, alpha0.8)) if k1_crit_vari: ax.axvline(xk1_crit_vari, colorred, linestyle:, linewidth1.5, alpha0.7) ax.annotate(fk1^变体{k1_crit_vari:.3f}, xy(k1_crit_vari, 10), xytext(k1_crit_vari0.005, 8), fontsize10, colorred, bboxdict(boxstyleround,pad0.3, facecolorlightcoral, alpha0.8)) ax.fill_between(k1_arr, 0, 10, alpha0.05, colorgreen, label稳态区) ax.set_xlabel($k_1$kernel_ratio / 耦合强度, fontsize12) ax.set_ylabel(灵帝 ξ_L 终值, fontsize12) ax.set_title(k1 扫描原版 vs 变体 ξ_L 相变边界, fontsize14) ax.legend() ax.grid(alpha0.3) plt.tight_layout() plt.savefig(k1_scan_phase_boundary.png, dpi150, bbox_inchestight) plt.show() print( 图表已保存k1_scan_phase_boundary.png) print(控制台表格 临界点分析可直接填入报告。)项目描述目标分析观测失真对系统相变阈值的影响通过对比原版和变体模型的仿真结果。方法对k1kernel_ratio进行扫描计算原版和变体模型的ξ_L终值并比较其差异。关键函数system_original和system_variant分别表示原版和变体模型其中变体模型通过theta引入观测失真。结果指标ξ_L终值、Δξ_L终值差异、加速效应%相对变化百分比。临界点分析通过线性插值法确定原版和变体模型的相变阈值k1_crit并计算其差异。结论观测失真显著降低了系统的相变阈值表明观测误差可能使系统更容易进入发散状态。参考来源numpy、scipy、pandas、matplotlib的读书报告手把手教你使用Numpy、Matplotlib、Scipy等5个Python库python如何安装Numpy、SciPy、MatPlotLib简述Python的Numpy,SciPy和Pandas,Matplotlib的区别简述 Python 的 Numpy、SciPy、Pandas、Matplotlib 的区别