import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # 定义微分方程组 def phosphorus model(y, t, params): x1, x2, x3 = y k01, k10, k12, k21, k23, k31, k30 = params dx1dt = k01 + k21∗x2 + k31∗x3 − (k10 + k12)∗x1 dx2dt = k12∗x1 − (k21 + k23)∗x2 dx3dt = k23∗x2 − (k31 + k30)∗x3 return [dx1dt, dx2dt, dx3dt] # 蒙特卡洛模拟参数设置 num simulations = 1000 # 模拟次数 simulation time = 30 # 模拟时长(天) t = np.linspace(0, simulation time, 100) # 时间点 # 实验数据的均值(示例数据,需根据实际情况调整) mean params = { ’k01’: 200, # 外部输入速率均值 ’k10’: 27.3, # 自然流失速率均值 ’k12’: 247.4, # 植物吸收速率均值 ’k21’: 9.1, # 植物释放速率均值 ’k23’: 238.4, # 动物摄食速率均值 ’k31’: 65.7, # 动物排放速率均值 ’k30’: 172.7 # 动物带走速率均值 } # 参数标准差(根据数据波动性设定) std dev = { ’k01’: 40, 附录 9 ’k10’: 5, ’k12’: 50, ’k21’: 2, ’k23’: 40, ’k31’: 10, ’k30’: 30 } # 初始化存储矩阵 results = np.zeros((num simulations, 3)) # 存储稳态磷浓度 # 蒙特卡洛模拟主循环 for i in range(num simulations): # 生成随机参数(正态分布,确保正数) params = [ np.abs(np.random.normal(mean params[’k01’], std dev[’k01’])), np.abs(np.random.normal(mean params[’k10’], std dev[’k10’])), np.abs(np.random.normal(mean params[’k12’], std dev[’k12’])), np.abs(np.random.normal(mean params[’k21’], std dev[’k21’])), np.abs(np.random.normal(mean params[’k23’], std dev[’k23’])), np.abs(np.random.normal(mean params[’k31’], std dev[’k31’])), np.abs(np.random.normal(mean params[’k30’], std dev[’k30’])) ] # 初始磷含量(可根据实际情况随机化) initial = [15.0, 1.8, 13.1] # x1, x2, x3 # 求解微分方程 solution = odeint(phosphorus model, initial, t, args=(params,)) # 记录稳态结果(取最后时间点的平均值)10% steady state = np.mean(solution[−len(t)//10:], axis=0) results[i] = steady state # 结果分析 print(”\模拟结果统计:n”) print(f”湖水磷浓度均值: 10 {np.mean(results[:,0]):.2f}±{np.std(results[:,0]):.2f}mg”) print(f”植物磷浓度均值: {np.mean(results[:,1]):.2f}±{np.std(results[:,1]):.2f}mg”) print(f”动物磷浓度均值: {np.mean(results[:,2]):.2f}±{np.std(results[:,2]):.2f}mg”) # 可视化分布 plt.figure(figsize=(12,6)) plt.subplot(131) plt.hist(results[:,0], bins=30, color=’blue’, alpha=0.7) plt.title(’湖水磷浓度分布’) plt.xl
最新发布