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 在更改最后一个因素时具有一致的效果。
我编写了一个代码来使用 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 在更改最后一个因素时具有一致的效果。