logo标签列表

如何在scipy.integrate中设置固定步长大小？

pythonscipyintegraterunge-kutta
4
我正在寻找一种方法，在Python中使用Runge-Kutta方法解决初始值问题时设置固定的步长。因此，我该如何告诉scipy.integrate.RK45保持其积分过程的恒定更新（步长）？非常感谢。
2
你想要实现什么目标？你希望算法像固定步长方法一样运行，还是希望得到等间距采样的解以便于绘图或进一步计算？ - Lutz Lehmann
感谢您的留言。实际上，我正在寻找前一部分，即固定步长法。 - behzad baghapour
由于这似乎是为教育目的而完成的，因此直接使用Dormand-Prince表格实现固定步长循环可能要容易得多。你可以使用RK45类来伪造它，在执行下一步之前不断重置步长，但没有明确保证只执行一步。 - Lutz Lehmann
1
我不知道你的动机是什么，但是作为一个集成模块的作者，在过去的几年里，我收到了十几个像你这样的请求。对于每一个请求，最终都证明这是由于完全误解步长自适应或者是一个畸形的XY问题，即固定步长是另一个问题的过于复杂或者糟糕的解决方案。因此，除非你真的知道你在做什么，否则我强烈建议你质疑自己是否真的需要这个（这反过来可能是一个好问题[scicomp.se]）。 - Wrzlprmft
6个回答
7

Scipy.integrate通常通过控制误差TOL（一步误差）来使用可变步长方法进行数值积分。通常通过与另一种数值方法进行比较来计算TOL。例如，RK45使用第5阶Runge-Kutta方法来检查第4阶Runge-Kutta方法的TOL以确定积分步骤。

因此，如果必须使用固定步长积分ODEs，则只需通过将atol，rtol设置为相当大的常量来关闭TOL检查。例如，如下形式：

solve_ivp(您的函数，t_span = [0, 10]，y0 = ...，method =“RK45”，max_step = 0.01，atol = 1，rtol = 1)

TOL检查被设置为非常大，以便积分步长将是您选择的max_step。

7

编写Dormand-Prince RK45方法的Butcher表格非常容易。

0
1/5   |  1/5
3/10  |  3/40        9/40
4/5   |  44/45       −56/15        32/9
8/9   |  19372/6561  −25360/2187   64448/6561   −212/729
1     |  9017/3168   −355/33       46732/5247   49/176     −5103/18656
1     |  35/384           0        500/1113     125/192    −2187/6784     11/84     
-----------------------------------------------------------------------------------------
      |  35/384           0        500/1113     125/192    −2187/6784     11/84     0
      |  5179/57600       0        7571/16695   393/640    −92097/339200  187/2100  1/40

首先，在单步函数中导入numpy库。

def DoPri45Step(f,t,x,h):

    k1 = f(t,x)
    k2 = f(t + 1./5*h, x + h*(1./5*k1) )
    k3 = f(t + 3./10*h, x + h*(3./40*k1 + 9./40*k2) )
    k4 = f(t + 4./5*h, x + h*(44./45*k1 - 56./15*k2 + 32./9*k3) )
    k5 = f(t + 8./9*h, x + h*(19372./6561*k1 - 25360./2187*k2 + 64448./6561*k3 - 212./729*k4) )
    k6 = f(t + h, x + h*(9017./3168*k1 - 355./33*k2 + 46732./5247*k3 + 49./176*k4 - 5103./18656*k5) )

    v5 = 35./384*k1 + 500./1113*k3 + 125./192*k4 - 2187./6784*k5 + 11./84*k6
    k7 = f(t + h, x + h*v5)
    v4 = 5179./57600*k1 + 7571./16695*k3 + 393./640*k4 - 92097./339200*k5 + 187./2100*k6 + 1./40*k7;

    return v4,v5

然后在一个标准的固定步长循环中执行

def DoPri45integrate(f, t, x0):
    N = len(t)
    x = [x0]
    for k in range(N-1):
        v4, v5 = DoPri45Step(f,t[k],x[k],t[k+1]-t[k])
        x.append(x[k] + (t[k+1]-t[k])*v5)
    return np.array(x)

然后，使用已知的精确解y(t)=sin(t)对其进行一些玩具示例的测试。

def mms_ode(t,y): return np.array([ y[1], sin(sin(t))-sin(t)-sin(y[0]) ])
mms_x0 = [0.0, 1.0]

并绘制由 h^5 缩放的误差图

for h in [0.2, 0.1, 0.08, 0.05, 0.01][::-1]:
    t = np.arange(0,20,h);
    y = DoPri45integrate(mms_ode,t,mms_x0)
    plt.plot(t, (y[:,0]-np.sin(t))/h**5, 'o', ms=3, label = "h=%.4f"%h);
plt.grid(); plt.legend(); plt.show()

为了确认这确实是一个五阶方法，需要注意误差系数的图形是否趋近。请参考下图：enter image description here
感谢您的出色工作，这节省了我很多时间。尽管在DoPri45integrate中设置数组切片似乎不可能，并且会产生非常恶心的错误。我将代码编辑为更健壮的列表附加。 - AlexNe
我不确定为什么你会在那里遇到错误，在我的 Python 2 和 Python 3 上都可以正常运行。 - Lutz Lehmann
我也不确定。可能是新的numpy版本或者维度不匹配。我在github上发布了代码：https://github.com/AlexanderNenninger/DynamicalSystems1-SS2020/blob/master/Van-der-PolOscillatorEx6Problem21.ipynb。我还使用调试器在jupyter笔记本之外测试了代码，切片值分配似乎是问题所在。 - AlexNe
1
你知道你使用的不是标准的Van der Pol振荡器x''-mu*(1-x^2)x'+x=0吗？你实现了方程eps*x'' - x'(1-x')^2+x=0。在它的导数y=x'中，它变成了eps*y''-y'(1-4y+3y^2)+y=0，经过时间重新调整后，它与VdP有些相似。 - Lutz Lehmann
感谢您的侦查工作，我之前并不知道！这是一项家庭作业任务，它们经常被错误地规定或包含错误。现在我会坚持按照所写的去做，并在备注中提到这一点。 - AlexNe
VdP振荡器来自于一个物理模型，因此你的公式可能仍更接近物理学，而不是通常使用的归一化方程。但是由于该方程不再对称，所以极限环也将不对称。 - Lutz Lehmann
3
如果您对数据固定步长感兴趣，我强烈建议您使用scipy.integrate.solve_ivp函数及其t_eval参数。
此函数包装了所有scipy.integrateode求解器的功能，因此您必须通过给其method参数赋值来选择方法。幸运的是，默认方法是RK45，所以您不必担心这个问题。
对您更有趣的是t_eval参数，在其中您必须提供一个平坦的数组。该函数在t_eval值处采样解曲线并仅返回这些点。因此，如果您想通过步长进行均匀采样，只需将t_eval参数设置为：numpy.linspace（t0，tf，samplingResolution），其中t0是模拟开始时间，tf是结束时间。 这样，您就可以获得均匀的采样，而无需使用会导致某些ODE不稳定的固定步长。
3

通过查看步骤的实现，您会发现最好的方法就是通过在调用RK45.step之前设置属性h_abs来控制初始步长（在最小和最大步长范围内）：

In [27]: rk = RK45(lambda t, y: t, 0, [0], 1e6)

In [28]: rk.h_abs = 30

In [29]: rk.step()

In [30]: rk.step_size
Out[30]: 30.0
1
您说过您想要一个固定时间步长的行为，而不仅仅是一个固定的评估时间步长。因此，如果您不想重新实现求解器，您必须通过“黑客”方式来解决这个问题。只需将积分公差atol和rtol设置为1e90，将max_step和first_step设置为所需时间步长dt的值。这样，估计的积分误差将始终非常小，从而欺骗求解器不会动态缩小时间步长。

然而，只有在显式算法(RK23、RK45、DOP853)中才使用这个技巧！"solve_ivp"中的隐式算法(Radau、BDF，可能还包括LSODA)根据atol和rtol调整非线性牛顿求解器的精度，因此您可能会得到一个毫无意义的解...

1
我建议你在Python中编写自己的rk4固定步长程序。有很多互联网上的例子可以帮助你。这可以确保你知道每个值是如何计算的。此外，通常不会出现0/0计算，如果出现，很容易追踪并提示重新审视正在解决的ODE。
也许分享一个入门的网站？ - Seth B
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接