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

mbyulnm0  于 5个月前  发布在  其他
关注(0)|答案(6)|浏览(88)

我正在寻找一种方法,在Python中设置一个固定步长来解决我的龙格-库塔方法的初值问题。因此,我如何告诉scipy.integrate.RK45保持一个恒定的更新(步长),为它的积分过程?
非常感谢

drnojrws

drnojrws1#

Scipy.integrate通常与变步长方法一起使用,通过控制数值积分时的TOL(一步误差)来计算。TOL通常通过与另一种数值方法进行检查来计算。例如RK 45使用五阶Runge-Kutta来检查四阶Runge-Kutta方法的TOL以确定积分步长。
因此,如果你必须集成固定步长的ODE,只需通过设置一个相当大的常数atol,rtol来关闭TOL检查。例如,像这样的形式:
第一个月
TOL检查设置得很大,以至于积分步长将是您选择的max_step。

jmo0nnb3

jmo0nnb32#

为Dormand-Prince RK 45方法编写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    
-----------------------------------------------------------------------------------------
 dy4   |  35/384           0        500/1113     125/192    −2187/6784     11/84     0
 dy5   |  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阶方法,因为误差系数的曲线图非常接近。
x1c 0d1x的数据

gopyfrb3

gopyfrb33#

通过查看步骤的实现,您会发现您可以做的最好的事情是通过在调用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

字符串

yizd12fk

yizd12fk4#

如果您对数据方面的修复步长感兴趣,那么我强烈建议您使用scipy.integrate.solve_ivp函数及其t_eval参数。
这个函数将所有的scipy.integrate ode求解器打包在一个函数中,因此你必须通过给它的method参数赋值来选择方法。幸运的是,默认方法是RK45,所以你不必为此烦恼。
对你来说更有趣的是t_eval参数,在这里你必须给予一个平面数组。该函数在t_eval值处对解曲线进行采样,并只返回这些点。所以如果你想按步长均匀采样,那么只需给予t_eval参数如下:numpy.linspace(t0, tf, samplingResolution),其中t0是模拟的开始,tf是模拟的结束。因此,您可以进行均匀采样,而不必采取固定步长,这会导致某些常微分方程的不稳定性。

g6ll5ycj

g6ll5ycj5#

你说过你想要一个固定的时间步长行为,而不仅仅是一个固定的评估时间步长。因此,如果你不想自己重新实现求解器,你必须“破解”你的方式。只需将积分容差atol和rtol设置为1 e90,将max_step和first_step设置为要使用的时间步长的值dt。这样,估计的积分误差将始终非常小,从而欺骗求解器不动态地收缩时间步长。
然而,只使用显式算法(RK 23,RK 45,DOP 853)!来自“solve_ivp”的隐式算法(Radau,BDF,也许还有LSODA)根据atol和rtol调整非线性牛顿求解器的精度,因此您可能最终得到一个没有任何意义的解决方案。

ehxuflar

ehxuflar6#

我建议你用py编写自己的rk 4固定步长程序。网上有很多例子可以帮助你。这保证了你精确地知道每个值是如何计算的。此外,通常不会有0/0计算,如果是这样的话,它们将很容易跟踪并提示你再次查看正在求解的ode。

相关问题