如何使用scipy.integrate设置固定步长? [英] How to set fixed step size with scipy.integrate?

问题描述

我正在寻找一种设置固定步长的方法，以通过Python中的Runge-Kutta方法解决我的初值问题.因此，如何告诉scipy.integrate.RK45对其集成过程进行恒定的更新(步长)?

I am looking for a way to set a fixed step size for solving my initial value problem by Runge-Kutta method in Python. Accordingly, how I can tell the scipy.integrate.RK45 to keep a constant update (step size) for its integration procedure?

非常感谢您.

推荐答案

为Dormand-Prince RK45方法编写Butcher表很容易.

It is quite easy to code the Butcher tableau for the Dormand-Prince RK45 method.

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导入为np

first in a function for a single step 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

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

and then in a standard fixed-step loop

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)

Then test it for some toy example with known exact solution 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阶方法，因为误差系数的图形非常接近.

to get the confirmation that this is indeed an order 5 method, as the graphs of the error coefficients come close together.

这篇关于如何使用scipy.integrate设置固定步长?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！