Runge Kutta Fehlberg 中公差的验证

Validation of tolerance in Runge Kutta Fehlberg

我编写了一个代码来使用 Runge Kutta Fehlberg 方法求解 ODE 系统:

import numpy as np
import matplotlib.pyplot as plt 
import numba
import time

start_time = time.clock()

@numba.jit()
def S1(x,u):
    alpha = 0.01
    C1 = 0.00096*alpha
    C2 = 865688 * alpha**2
    C3 = 0.71488 * alpha**(3/2)
    C4 = 1715.2 * alpha**(3/2)
    C5 = 0.356
    C6 = 3.18373* 10**8 *alpha
    C7 = 262.9 * 10**20 * alpha**0.5
    C8 = 63* 10**25 * alpha**0.5
    Gam0 = 121
    v0 = 10**(-5)
    b = 2* 10**(-4)
    x0 = 0.00045
    B0 = 0.0005
    BY,nB,neL,neR=u
    B= B0/(b*(2*np.pi)**0.5) * np.exp(-((x-x0)**2)/(2 * b**2))
    v= v0/(b*(2*np.pi)**0.5)* np.exp(-((x-x0)**2)/(2 * b**2))
    nT= neR-neL/2 + (3/8) * nB
    dn= (neR**2-neL**2) * 0.5
    u1= (-C1-C2 * nT)* (BY/(10e+20))**2 * x**3/2
    u2= (C3*B+C4* (dn**2))*v*(BY/(10e+20))**2
    u3= (Gam0 * (1-x)/(x**0.5))*(neR-neL)
    dneL= -(1/4) * (u1+u2) + u3/2
    dneR= u1+u2-u3
    dnB= (3/2)*(u1+u2)
    dBY= (1/(x**0.5))*(-C5-C6*nT)*BY - BY/x + (C7*B+C8* dn**2) * (v/(x**1.5))
    return np.array([dBY,dnB,dneL,dneR])


@numba.jit()
def rkf( f, a, b, x0, tol, hmax, hmin ):

    a2  =   2.500000000000000e-01  #  1/4
    a3  =   3.750000000000000e-01  #  3/8
    a4  =   9.230769230769231e-01  #  12/13
    a5  =   1.000000000000000e+00  #  1
    a6  =   5.000000000000000e-01  #  1/2

    b21 =   2.500000000000000e-01  #  1/4
    b31 =   9.375000000000000e-02  #  3/32
    b32 =   2.812500000000000e-01  #  9/32
    b41 =   8.793809740555303e-01  #  1932/2197
    b42 =  -3.277196176604461e+00  # -7200/2197
    b43 =   3.320892125625853e+00  #  7296/2197
    b51 =   2.032407407407407e+00  #  439/216
    b52 =  -8.000000000000000e+00  # -8
    b53 =   7.173489278752436e+00  #  3680/513
    b54 =  -2.058966861598441e-01  # -845/4104
    b61 =  -2.962962962962963e-01  # -8/27
    b62 =   2.000000000000000e+00  #  2
    b63 =  -1.381676413255361e+00  # -3544/2565
    b64 =   4.529727095516569e-01  #  1859/4104
    b65 =  -2.750000000000000e-01  # -11/40

    r1  =   2.777777777777778e-03  #  1/360
    r3  =  -2.994152046783626e-02  # -128/4275
    r4  =  -2.919989367357789e-02  # -2197/75240
    r5  =   2.000000000000000e-02  #  1/50
    r6  =   3.636363636363636e-02  #  2/55

    c1  =   1.157407407407407e-01  #  25/216
    c3  =   5.489278752436647e-01  #  1408/2565
    c4  =   5.353313840155945e-01  #  2197/4104
    c5  =  -2.000000000000000e-01  # -1/5

    t = a
    x = np.array(x0)
    h = hmax

    T = np.array( [t] )
    X = np.array( [x] )
    
    while t < b:

        if t + h > b:
            h = b - t

        k1 = h * f(t, x)
        k2 = h * f(t + a2 * h, x + b21 * k1 )
        k3 = h * f(t + a3 * h, x + b31 * k1 + b32 * k2)
        k4 = h * f(t + a4 * h, x + b41 * k1 + b42 * k2 + b43 * k3)
        k5 = h * f(t + a5 * h, x + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4)
        k6 = h * f(t + a6 * h, x + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5)
        # print(t)
        r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
        if len( np.shape( r ) ) > 0:
            r = max( r )
        if r <= tol:
            t = t + h
            x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
            T = np.append( T, t )
            X = np.append( X, [x], 0 )

        h = h * min( max( 0.84 * ( tol / r )**0.25, 0.1 ), 4.0 )
        
        if h > hmax:
            h = hmax
        elif h < hmin:
            raise RuntimeError("Error: Could not converge to the required tolerance %e with minimum stepsize  %e." % (tol,hmin))
            break

    return (X, T)
    
S1u0=np.array([0,0,0,0])
np.seterr(all='ignore')

u,t =rkf( f=S1, a=1e-4, b=1, x0=S1u0, tol=1e+10, hmax=1e-1, hmin=1e-40)

print("Execution time:",time.clock() - start_time, "seconds")
BY,nB,neL,neR = u.T
print(BY.shape)
plt.plot(t,BY)
plt.xscale('log')
plt.yscale('log')
plt.show()

运行速度非常快,returns 得到了想要的结果:

(我在 Maple 中解决了完全相同的 ODE,并得到了相同的结果)。但我的问题是关于 公差 。当我将公差设置为 1e+10 时,代码运行得非常快(它只集成了几个点的 ODE)并且与较低的公差相反,它需要更长的时间(即使是 1e+4 的公差)。我的第一个问题是 为什么我会得到像 1e+10 这样的公差的正确答案? 我的第二个问题是如何计算相对和绝对误差?第三,是否可以使此代码更高效?

你在积分中得到这么多点,因为组件有不同的比例,第一个大约在 1e+20,其他三个大约在 1e-10。使用由绝对和相对误差容限构造的每个组件的步长可以获得更好的结果。

def rkf( f, a, b, x0, atol, rtol, hmax, hmin ):
    ...
    while t < b:

        ...

        r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
        r = r / (atol+rtol*(abs(x)+abs(k1)))
        if len( np.shape( r ) ) > 0:
            r = max( r )
        if r <= 1:
            t = t + h
            x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
            T = np.append( T, t )
            X = np.append( X, [x], 0 )

        h = h * min( max( 0.94 * ( 1 / r )**0.25, 0.1 ), 4.0 )
        
        if h > hmax:
            h = hmax
        elif h < hmin or t==t-h:
            raise RuntimeError("Error: Could not converge to the required tolerance %e with minimum stepsize  %e. t=%g h=%e" % (tol,hmin,t,h))
            break

合理公差不同导致点数减少,如

atol=1e+4, rtol=1e-2  ==>  93 points
atol=1e+0, rtol=1e-4  ==> 125 points
atol=1e+10, rtol=1e-6 ==> 245 points

点数在很大程度上对 atol 值不敏感。也许这需要为每个组件设置,设置 atol=np.array([1e+20,1e-10,1e-10,1e-10])*1e-4 在更改最后一个因素时具有一致的效果。