ODE数值求解问题：y'=2x在x=10处结果偏差排查

问题澄清 #

首先纠正一个关键误解：给定初始条件x=0时y=1，ODE y'=2x的解析解应为y(x)=x²+1，因此x=10处的正确结果是101，而非你认为的100。你得到的数值结果101.000000001241属于数值积分过程中的极小浮点误差，是正常现象，以下是优化精度的具体方法：

优化方案 #

1. 避免冗余插值，直接获取结果 #

你的代码中使用np.interp(10, x, y)完全多余——np.linspace(0,10)生成的最后一个点就是x=10，直接取y[-1]即可得到对应数值解，省去插值带来的额外误差。

2. 调整求解器精度参数 #

odeint支持通过rtol（相对误差容忍度）和atol（绝对误差容忍度）控制求解精度，默认值为1.49012e-8，设置更小的阈值可以进一步缩小偏差：

solution = odeint(deriv, y0, x, rtol=1e-12, atol=1e-12)

3. 确认导数函数参数顺序 #

odeint要求导数函数的第一个参数是状态变量y，第二个是自变量x，你的代码def deriv(y,x)符合要求，但若顺序颠倒会导致严重错误，需注意保持正确格式。

修正后的完整代码 #

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

def deriv(y, x):
    return 2 * x

# 生成自变量点
x = np.linspace(0, 10)
# 初始条件
y0 = 1

# 高精度求解ODE
solution = odeint(deriv, y0, x, rtol=1e-12, atol=1e-12)
y = solution[:, 0]

# 绘制结果
plt.plot(x, y)
plt.title("ODE y' = 2x 的数值解")
plt.ylabel("y值")
plt.xlabel("x")
plt.show()

# 直接获取x=10处的数值解
print("x=10处的数值解：", y[-1])
# 输出解析解作为对比
print("x=10处的解析解：", 10**2 + 1)

运行上述代码后，数值解与解析解的偏差会进一步缩小，可忽略不计。

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