基于ODE解的含噪声数据拟合：ODE系统曲线拟合与绘图

核心思路 #

我们需要将带有5%高斯噪声的ODE解数据，拟合回原始ODE系统——本质是通过优化算法找到能最小化预测值与噪声数据误差的ODE参数，最后绘制拟合曲线与原始数据对比。

完整实现代码 #

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

# 替换为你实际的ODE系统（示例为双参数线性ODE）
def dose(y, D, params):
    k1, k2 = params
    # 根据你的实际ODE逻辑修改微分方程
    dy_dt = -k1 * y + k2 * (D - y)
    return dy_dt

# 加载你的D和mat2数据（这里模拟示例数据，实际中直接替换为你的数组即可）
true_params = [0.5, 0.2]  # 仅用于模拟真实数据，实际无需此变量
D = np.linspace(0, 10, 7)  # 你的7个D值
y0 = 0.0  # 替换为你实际的ODE初始条件

# 定义从ODE提取特征值的函数（按需修改：比如取稳态值、特定时间点值等）
def get_response(D_vals, params):
    responses = []
    for d in D_vals:
        # 求解ODE，这里取t=10时的稳态值作为特征值
        y_sol = odeint(dose, y0, [0, 10], args=(d, params))[-1]
        responses.append(y_sol[0])
    return np.array(responses)

# 生成带5%噪声的模拟数据（实际中直接用你的mat2变量）
true_mat = get_response(D, true_params)
noise_level = 0.05
mat2 = true_mat + noise_level * np.std(true_mat) * np.random.randn(len(true_mat))

# 拟合误差函数：计算预测值与噪声数据的残差
def residual_func(params, D_data, y_data):
    pred_vals = get_response(D_data, params)
    return pred_vals - y_data

# 初始化参数猜测值（根据ODE物理意义设置合理初始值）
initial_guess = [0.4, 0.1]
# 执行非线性拟合
fit_result = least_squares(residual_func, initial_guess, args=(D, mat2))
fitted_params = fit_result.x

# 打印拟合结果
print(f"拟合得到的参数：{fitted_params.round(4)}")
print(f"真实参数（模拟用）：{true_params}")

# 绘制结果对比图
D_dense = np.linspace(D.min(), D.max(), 100)
true_curve = get_response(D_dense, true_params)
fitted_curve = get_response(D_dense, fitted_params)

plt.figure(figsize=(10, 6))
plt.scatter(D, mat2, color='#ff4c4c', s=60, label='带噪声数据')
plt.plot(D_dense, true_curve, color='#2e8b57', linestyle='--', label='真实ODE曲线')
plt.plot(D_dense, fitted_curve, color='#1e90ff', label='拟合ODE曲线')
plt.xlabel('D值')
plt.ylabel('特征值')
plt.legend(frameon=True, loc='best')
plt.title('ODE系统拟合结果对比')
plt.grid(alpha=0.3)
plt.show()

关键注意事项 #

1. ODE系统替换：将示例中的dose函数替换为你实际的ODE定义，确保参数params包含所有需要拟合的未知变量，微分方程逻辑与原始系统完全一致。
2. 特征值提取逻辑：如果你的特征值不是稳态值（比如某个特定时间点的解、积分值等），修改get_response函数里的ODE求解部分，比如通过t_eval参数指定提取对应时间点的解。
3. 初始参数选择：初始猜测值initial_guess尽量贴近真实参数范围（可通过物理意义或初步试算获取），避免优化算法陷入局部最优解。
4. 优化器选择：示例用least_squares处理非线性拟合，比curve_fit更适合多参数、残差非线性的场景；如果你的ODE是线性的，也可以改用线性回归类方法提升效率。

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