我正在寻找一种方法,在Python中设置一个固定步长来解决我的龙格-库塔方法的初值问题。因此,我如何告诉scipy.integrate.RK45保持一个恒定的更新(步长),为它的积分过程?非常感谢
scipy.integrate.RK45
drnojrws1#
Scipy.integrate通常与变步长方法一起使用,通过控制数值积分时的TOL(一步误差)来计算。TOL通常通过与另一种数值方法进行检查来计算。例如RK 45使用五阶Runge-Kutta来检查四阶Runge-Kutta方法的TOL以确定积分步长。因此,如果你必须集成固定步长的ODE,只需通过设置一个相当大的常数atol,rtol来关闭TOL检查。例如,像这样的形式:第一个月TOL检查设置得很大,以至于积分步长将是您选择的max_step。
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
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 npdef 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
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)
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)来测试它
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]
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
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()
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的数据
gopyfrb33#
通过查看步骤的实现,您会发现您可以做的最好的事情是通过在调用RK45.step之前设置属性h_abs来控制初始步长(在最小和最大步长设置的范围内):
RK45.step
h_abs
In [27]: rk = RK45(lambda t, y: t, 0, [0], 1e6)In [28]: rk.h_abs = 30In [29]: rk.step()In [30]: rk.step_sizeOut[30]: 30.0
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
字符串
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是模拟的结束。因此,您可以进行均匀采样,而不必采取固定步长,这会导致某些常微分方程的不稳定性。
scipy.integrate.solve_ivp
t_eval
scipy.integrate
method
numpy.linspace(t0, tf, samplingResolution)
g6ll5ycj5#
你说过你想要一个固定的时间步长行为,而不仅仅是一个固定的评估时间步长。因此,如果你不想自己重新实现求解器,你必须“破解”你的方式。只需将积分容差atol和rtol设置为1 e90,将max_step和first_step设置为要使用的时间步长的值dt。这样,估计的积分误差将始终非常小,从而欺骗求解器不动态地收缩时间步长。然而,只使用显式算法(RK 23,RK 45,DOP 853)!来自“solve_ivp”的隐式算法(Radau,BDF,也许还有LSODA)根据atol和rtol调整非线性牛顿求解器的精度,因此您可能最终得到一个没有任何意义的解决方案。
ehxuflar6#
我建议你用py编写自己的rk 4固定步长程序。网上有很多例子可以帮助你。这保证了你精确地知道每个值是如何计算的。此外,通常不会有0/0计算,如果是这样的话,它们将很容易跟踪并提示你再次查看正在求解的ode。
6条答案
按热度按时间drnojrws1#
Scipy.integrate通常与变步长方法一起使用,通过控制数值积分时的TOL(一步误差)来计算。TOL通常通过与另一种数值方法进行检查来计算。例如RK 45使用五阶Runge-Kutta来检查四阶Runge-Kutta方法的TOL以确定积分步长。
因此,如果你必须集成固定步长的ODE,只需通过设置一个相当大的常数atol,rtol来关闭TOL检查。例如,像这样的形式:
第一个月
TOL检查设置得很大,以至于积分步长将是您选择的max_step。
jmo0nnb32#
为Dormand-Prince RK 45方法编写Butcher表非常容易。
字符串
在单个步骤的函数中的第一个
型
然后在一个标准的固定步长循环中
型
然后用已知的精确解
y(t)=sin(t)
来测试它型
并将误差标绘为
h^5
型
以确认这确实是5阶方法,因为误差系数的曲线图非常接近。
x1c 0d1x的数据
gopyfrb33#
通过查看步骤的实现,您会发现您可以做的最好的事情是通过在调用
RK45.step
之前设置属性h_abs
来控制初始步长(在最小和最大步长设置的范围内):字符串
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是模拟的结束。因此,您可以进行均匀采样,而不必采取固定步长,这会导致某些常微分方程的不稳定性。g6ll5ycj5#
你说过你想要一个固定的时间步长行为,而不仅仅是一个固定的评估时间步长。因此,如果你不想自己重新实现求解器,你必须“破解”你的方式。只需将积分容差atol和rtol设置为1 e90,将max_step和first_step设置为要使用的时间步长的值dt。这样,估计的积分误差将始终非常小,从而欺骗求解器不动态地收缩时间步长。
然而,只使用显式算法(RK 23,RK 45,DOP 853)!来自“solve_ivp”的隐式算法(Radau,BDF,也许还有LSODA)根据atol和rtol调整非线性牛顿求解器的精度,因此您可能最终得到一个没有任何意义的解决方案。
ehxuflar6#
我建议你用py编写自己的rk 4固定步长程序。网上有很多例子可以帮助你。这保证了你精确地知道每个值是如何计算的。此外,通常不会有0/0计算,如果是这样的话,它们将很容易跟踪并提示你再次查看正在求解的ode。