3

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

非常感谢。

4

6 回答 6

6

为 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

首先在一个函数中进行单步 import numpy as np

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()

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

在此处输入图像描述

于 2019-02-03T12:23:39.927 回答
4

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

因此，如果您必须将 ODE 与固定步长集成，只需通过将 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。

于 2019-09-15T08:31:54.827 回答
3

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

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
于 2019-02-02T16:15:41.177 回答
2

如果您对数据修复步长感兴趣，那么我强烈建议您使用该scipy.integrate.solve_ivp函数及其t_eval参数。

此函数将所有 ode 求解器封装scipy.integrate在一个函数中，因此您必须通过为其method参数赋值来选择方法。幸运的是，默认方法是 RK45，因此您不必为此烦恼。

对你来说更有趣的是t_eval参数，你必须给出一个平面数组。该函数在值处对解曲线进行采样，t_eval并仅返回这些点。因此，如果您想要按步长进行统一采样，那么只需给出以下t_eval参数：numpy.linspace(t0, tf, samplingResolution)，其中 t0 是模拟的开始，tf 是模拟的结束。因此，您可以进行统一采样，而不必采用会导致某些 ODE 不稳定的固定步长。

于 2020-10-23T06:49:14.727 回答
1

我建议在py中编写自己的rk4定步程序。有许多互联网示例可以提供帮助。这可以保证您准确地知道每个值是如何计算的。此外，通常不会有 0/0 计算，如果是这样，它们将很容易追踪并提示再次查看正在求解的 ode。

于 2021-02-13T21:53:33.120 回答
1

你说过你想要一个固定时间步的行为，而不仅仅是一个固定的评估时间步。因此，如果您不想自己重新实现求解器，则必须“破解”您的方式。只需将积分公差 atol 和 rtol 设置为 1e90，并将 max_step 和 first_step 设置为您要使用的时间步长的值 dt。这样，估计的积分误差总是很小，从而欺骗求解器不动态地缩小时间步长。

但是，只能将此技巧与 EXPLICIT 算法(RK23,RK45,DOP853) 一起使用！来自“solve_ivp”的隐式算法（Radau、BDF，也可能是 LSODA）根据 atol 和 rtol 调整非线性牛顿求解器的精度，因此您最终可能会得到一个没有任何意义的解决方案......

于 2020-11-10T17:03:38.550 回答