如何使用 scipy.integrate 设置固定步长？

Question

我正在寻找一种方法来设置固定步长，以通过 Python 中的 Runge-Kutta 方法解决我的初始值问题。因此，我如何告诉scipy.integrate.RK45其集成过程保持不断更新（步长）？

非常感谢。

Answer 1

为 Dormand-Prince RK45 方法编写 Butcher 画面非常容易。

\n

\n 0     |\n 1/5   |  1/5\n 3/10  |  3/40        9/40\n 4/5   |  44/45       \xe2\x88\x9256/15        32/9\n 8/9   |  19372/6561  \xe2\x88\x9225360/2187   64448/6561   \xe2\x88\x92212/729\n 1     |  9017/3168   \xe2\x88\x92355/33       46732/5247   49/176     \xe2\x88\x925103/18656\n 1     |  35/384           0        500/1113     125/192    \xe2\x88\x922187/6784     11/84    \n-----------------------------------------------------------------------------------------\n dy4   |  35/384           0        500/1113     125/192    \xe2\x88\x922187/6784     11/84     0\n dy5   |  5179/57600       0        7571/16695   393/640    \xe2\x88\x9292097/339200  187/2100     1/40\n

\n

函数中的第一个单步

\n

import numpy as np\n\ndef DoPri45Step(f,t,x,h):\n    \n    k1 = f(t,x)\n    k2 = f(t + 1./5*h, x + h*(1./5*k1) )\n    k3 = f(t + 3./10*h, x + h*(3./40*k1 + 9./40*k2) )\n    k4 = f(t + 4./5*h, x + h*(44./45*k1 - 56./15*k2 + 32./9*k3) )\n    k5 = f(t + 8./9*h, x + h*(19372./6561*k1 - 25360./2187*k2 + 64448./6561*k3 - 212./729*k4) )\n    k6 = f(t + h, x + h*(9017./3168*k1 - 355./33*k2 + 46732./5247*k3 + 49./176*k4 - 5103./18656*k5) )\n\n    v5 = 35./384*k1 + 500./1113*k3 + 125./192*k4 - 2187./6784*k5 + 11./84*k6\n    k7 = f(t + h, x + h*v5)\n    v4 = 5179./57600*k1 + 7571./16695*k3 + 393./640*k4 - 92097./339200*k5 + 187./2100*k6 + 1./40*k7;\n    \n    return v4,v5\n

\n

然后在标准的固定步长循环中

\n

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

\n

然后用已知的精确解决方案测试一些玩具示例y(t)=sin(t)

\n

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

\n

并绘制按比例缩放的误差h^5

\n

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

\n

确认这确实是 5 阶方法，因为误差系数的图形非常接近。

\n

\n

Answer 2

Scipy.integrate 通常与可变步长方法一起使用，通过在数值积分时控制 TOL（一步误差）。TOL 通常是通过使用另一种数值方法进行检查来计算的。例如RK45使用5阶龙格-库塔法检查4阶龙格-库塔法的TOL来确定积分步长。

因此，如果必须对固定步长的常微分方程进行积分，只需通过将 atol、rtol 设置为相当大的常数来关闭 TOL 检查即可。例如，像这样的形式：

solve_ivp(your function, t_span=[0, 10], y0=..., method="RK45", max_step=0.01, atol = 1, rtol = 1)

TOL 检查设置得很大，以致于积分步长将是您选择的 max_step。