确定求解 R 中的 ode 的正确时间步长 (deSolve)

Question

如中所述，求解 ODE 时的 times 参数改变了结果。

基本上，对于完全相同的 func，seq(1, 10, 0.01) 和 seq(1, 10, 1) 的 times 会产生完全不同的结果。

我的问题是："right" 结果是什么？如何确定"right"时间步长（0.01或1，或任何其他值）?

我会逐渐减少时间步长，解决这个问题，看看结果是否停止变化。但是有更好的方法吗？

谢谢！

Answer 1

link 中的差异是混沌系统在小于 30 的时间跨度内将最小误差放大到混沌吸引子大小的结果。即使比较100 个子区间，其中一个有 101 个子区间。因此，对于这个系统，在大于 20 的时间跨度内计算出可靠的结果是完全不可能的。

python scipy 的 odeint 使用与 R 包相同的 lsode Fortran 代码。使用与 linked question/answer 中的测试示例相同的参数和函数，我通过将错误容限设置得相当低

atol, rtol = 1e-3, 1e-6

for N in [10, 100, 1000]:
    t = np.linspace(0,10,N+1)
    u = odeint(Lorenz,u0,t, atol=atol, rtol=rtol);
    print "N =%6d"%N, ", (t,u) =",zip(t,u)[-1]

输出

N =    10 , (t,u) = (10.0, array([ 15.9506689 ,  -6.49172355, -10.50061322]))
N =   100 , (t,u) = (10.0, array([ 15.86686806,  -6.39874131, -10.3567592 ]))
N =  1000 , (t,u) = (10.0, array([ 15.87163076,  -6.40449548, -10.36581067]))

仍然可以看到分段的影响，但是差异与误差容限成正比。

通过将误差容限设置为更合理的大小，输出中的差异会减少然后消失，因为积分器的内部步长变得比输出时间序列的时间步长小得多：

atol, rtol = 1e-6, 1e-8

N =    10 , (t,u) = (10.0, array([ 16.76057876,  -7.28291262, -11.68849765]))
N =   100 , (t,u) = (10.0, array([ 16.76049974,  -7.28284578, -11.68840157]))
N =  1000 , (t,u) = (10.0, array([ 16.76049991,  -7.28284592, -11.68840176]))

---
atol, rtol = 1e-12, 1e-14

N =    10 , (t,u) = (10.0, array([ 16.76043932,  -7.28277217, -11.68828488]))
N =   100 , (t,u) = (10.0, array([ 16.76043932,  -7.28277217, -11.68828488]))
N =  1000 , (t,u) = (10.0, array([ 16.76043932,  -7.28277217, -11.68828488]))

因此，如果您看到输出时间分段对积分结果的影响，则表明误差容限（给定或默认）不够严格。

确定求解 R 中的 ode 的正确时间步长 (deSolve)

Determine the right time step for solving ode in R (deSolve)

r

ode

differential-equations