问题描述
我编写了一个代码来使用 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])
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,neR = u.T
print(BY.shape)
plt.plot(t,BY)
plt.xscale('log')
plt.yscale('log')
plt.show()
运行速度非常快，并返回所需的结果：
（我已经在 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-4 在更改最后一个因素时具有一致的效果。