Lotka-Volterra模型是一种经典的非线性微分方程组，用于描述捕食者与猎物之间的种群动态关系。在实际应用中，研究人员常常需要根据观测数据对模型参数进行估计，以提高模型的预测能力。非线性最小二乘法是一种常用的参数估计方法，适用于这种非线性系统。 ### 数学原理 非线性最小二乘法的目标是通过最小化模型输出与实际观测值之间的残差平方和来估计模型参数。设观测数据为 $ (x_i, y_i) $，其中 $ i = 1, 2, \dots, n $，模型的预测值为 $ \hat{x}(t; \theta), \hat{y}(t; \theta) $，其中 $ \theta $ 表示待估计的参数向量。目标函数定义为： $$ S(\theta) = \sum_{i=1}^{n} \left[ (x_i - \hat{x}(t_i; \theta))^2 + (y_i - \hat{y}(t_i; \theta))^2 \right] $$ 通过优化算法（如Levenberg-Marquardt算法）可以求解使 $ S(\theta) $ 最小的参数 $ \theta $ [^4]。 ### Python 实现 以下是一个使用 `scipy.optimize.curve_fit` 实现 Lotka-Volterra 模型参数估计的完整代码示例： ```python import numpy as np from scipy.integrate import solve_ivp from scipy.optimize import curve_fit import matplotlib.pyplot as plt # Lotka-Volterra 模型定义 def lotka_volterra(t, z, a, b, c, d): x, y = z # x: prey, y: predator dxdt = x * (a - b * y) dydt = y * (c * x - d) return [dxdt, dydt] # 解微分方程 def solve_lotka_volterra(params, t_span, z0): sol = solve_ivp(lambda t, z: lotka_volterra(t, z, *params), t_span, z0, dense_output=True) return sol.sol # 生成模拟数据（用于演示） def generate_data(true_params, t_eval, noise_level=0.1): z0 = [10, 5] # 初始种群 sol = solve_lotka_volterra(true_params, [t_eval[0], t_eval[-1]], z0) x_true, y_true = sol(t_eval) x_data = x_true + np.random.normal(0, noise_level, len(t_eval)) y_data = y_true + np.random.normal(0, noise_level, len(t_eval)) return x_data, y_data # 非线性拟合函数 def model_func(t, a, b, c, d): z0 = [10, 5] sol = solve_lotka_volterra([a, b, c, d], [t[0], t[-1]], z0) x_sol, y_sol = sol(t) return np.concatenate((x_sol, y_sol)) # 主程序 if __name__ == "__main__": true_params = [1.5, 0.1, 0.2, 0.5] # 真实参数 t_eval = np.linspace(0, 20, 100) x_data, y_data = generate_data(true_params, t_eval) # 构造输入输出数据 data_input = np.concatenate((x_data, y_data)) # 初始猜测 initial_guess = [1.0, 0.1, 0.1, 0.5] # 参数估计 params_opt, params_cov = curve_fit(model_func, t_eval, data_input, p0=initial_guess) print("Estimated parameters:", params_opt) print("True parameters:", true_params) # 可视化结果 sol_opt = solve_lotka_volterra(params_opt, [t_eval[0], t_eval[-1]], [10, 5]) x_opt, y_opt = sol_opt(t_eval) plt.figure(figsize=(12, 6)) plt.plot(t_eval, x_data, 'b.', label='Prey data') plt.plot(t_eval, y_data, 'r.', label='Predator data') plt.plot(t_eval, x_opt, 'b-', label='Prey fit') plt.plot(t_eval, y_opt, 'r-', label='Predator fit') plt.legend() plt.xlabel('Time') plt.ylabel('Population') plt.title('Lotka-Volterra Model Fitting with Nonlinear Least Squares') plt.show() ``` ### 说明 - 使用 `solve_ivp` 对 Lotka-Volterra 微分方程进行数值积分。 - 使用 `curve_fit` 进行非线性最小二乘拟合，自动调用 Levenberg-Marquardt 算法。 - 通过构造目标函数 `model_func`，将模型输出与观测数据进行比较。 - 最终输出估计的参数并绘制拟合曲线与原始数据对比图。 这种方法可以有效地估计 Lotka-Volterra 模型中的参数，并且适用于其他类似的非线性动力系统模型。 ---