首先导入数据

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as pltdata = np.array([[30, 4],[47.2, 6.1],[70.2, 9.8],[77.4, 35.2],[36.3, 59.4],[20.6, 41.7],[18.1, 19],[21.4, 13],[22, 8.3],[25.4, 9.1],[27.1, 7.4],[40.3, 8],[57, 12.3],[76.6, 19.5],[52.3, 45.7],[19.5, 51.1],[11.2, 29.7],[7.6, 15.8],[14.6, 9.7],[16.2, 10.1],[24.7, 8.6]
])

设置参数和定义函数

# Set initial parameters to small random values or default values
initial_guess = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]# Define ordinary differential equations
def model(y, t, a, b, c, d, e, f):x, z = ydxdt = a * x - b * x * z - e * x**2dydt = -c * z + d * x * z - f * z**2return [dxdt, dydt]#Define error function
def error(params):a, b, c, d, e, f = paramst = np.linspace(0, 1, len(data))  y0 = [data[0, 0], data[0, 1]]  y_pred = odeint(model, y0, t, args=(a, b, c, d, e, f))x_pred, z_pred = y_pred[:, 0], y_pred[:, 1]error_x = np.sum((x_pred - data[:, 0])**2)error_y = np.sum((z_pred - data[:, 1])**2)total_error = error_x + error_yreturn total_error# Use least squares method to fit parameters
result = minimize(error, initial_guess, method='Nelder-Mead')# Get the fitting parameter values
a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

 

进一步优化

import matplotlib.pyplot as plt# Different optimization methods
methods = ['Nelder-Mead', 'Powell', 'BFGS', 'L-BFGS-B']for method in methods:result = minimize(error, initial_guess, method=method)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectoryt_fit = np.linspace(0, 1, len(data))y0_fit = [data[0, 0], data[0, 1]]y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data pointsplt.figure(figsize=(10, 6))plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')plt.plot(x_fit, z_fit, label=f'Fitting results ({method})', linestyle='-', color='red')plt.xlabel('x')plt.ylabel('y')plt.legend()plt.grid()plt.show()print(f"Parameter values fitted using {method} method：a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

 

Parameter values fitted using Nelder-Mead method：a=1.2173283165346425, b=0.42516102725023064, c=19.726779624261006, d=0.7743814851338301, e=-0.19482192444374966, f=0.37455729849779884

 

Parameter values fitted using Powell method：a=32.49329459442917, b=0.6910719576651617, c=58.98701472032894, d=1.3524516626786816, e=0.47787798383104335, f=-0.5344483269192019

Parameter values fitted using BFGS method：a=1.2171938888848015, b=0.04968374479958104, c=0.9234835772585344, d=0.947268540340848, e=0.010742224447412019, f=0.7985132960108715

 

Parameter values fitted using L-BFGS-B method：a=1.154759061832585, b=0.32168624538800344, c=0.9455699334793284, d=0.9623931795647013, e=0.2936335531513881, f=0.8566315817923148

进一步优化

#Set parameter search range (interval)
bounds = [(-5, 25), (-5, 25), (-5, 25), (-5, 25), (-5, 25), (-5, 25)]# Use interval search to set initial parameter values
result = differential_evolution(error, bounds)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()print(f"Fitting parameter values：a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

Fitting parameter values：a=-0.04946907994199101, b=5.5137169943224755, c=0.6909170053541569, d=10.615879287885402, e=-3.1585499451409937, f=18.4110095977882
发现效果竟然变差了

#Set parameter search range (interval)
bounds = [(-0.1, 10), (-0.1, 10), (-0.1, 10), (-0.1,10), (-0.1, 10), (-0.1, 10)]# Use interval search to set initial parameter values
result = differential_evolution(error, bounds)a_fit, b_fit, c_fit, d_fit, e_fit, f_fit = result.x# Simulate the fitted trajectory
t_fit = np.linspace(0, 1, len(data))
y0_fit = [data[0, 0], data[0, 1]]
y_fit = odeint(model, y0_fit, t_fit, args=(a_fit, b_fit, c_fit, d_fit, e_fit, f_fit))
x_fit, z_fit = y_fit[:, 0], y_fit[:, 1]# Plot the fitting results against the original data points
plt.figure(figsize=(10, 6))
plt.scatter(data[:, 0], data[:, 1], label='Raw data', marker='o', color='blue')
plt.plot(x_fit, z_fit, label='Fitting results', linestyle='-', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()print(f"Fitting parameter values：a={a_fit}, b={b_fit}, c={c_fit}, d={d_fit}, e={e_fit}, f={f_fit}")

 最终优化结果为：

Fitting parameter values：a=10.0, b=0.6320729493793303, c=10.0, d=0.4325244090515547, e=-0.07495645186059174, f=0.18793803443302332

完整代码

创作不易，希望大家多点赞关注评论！！！