Validation of tolerance in Runge Kutta Fehlberg
I've written a code to solve a system of ODEs using Runge Kutta Fehlberg method:
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()
Which works very fast, and returns the desired result:
( I have solved the exact same ODEs in Maple, and got the same result). But my question is regarding the tolerance. When I set the tolerance to 1e+10 , code works really fast (it integrates ODEs with only few points) and in contrary with lower tolerance, well, it takes much much longer (even for tolerance of 1e+4). My first question is why I'm getting the right answer with a tolerance like 1e+10 ? My second question is how can I calculate relative and absolute error? And third, is it possible to make this code more efficient?
You get that many points in the integration because the components have different scales, the first at about 1e+20, the other three at about 1e-10. One gets better results with a per-component step-size constructed from absolute and relative error tolerances.
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
Different reasonable tolerances result in a reduced number of points, such as
atol=1e+4, rtol=1e-2 ==> 93 points
atol=1e+0, rtol=1e-4 ==> 125 points
atol=1e+10, rtol=1e-6 ==> 245 points
The number of points is largely insensitive to the atol value. Perhaps that needs to be set per component, setting atol=np.array([1e+20,1e-10,1e-10,1e-10])*1e-4 has a consistent effect when changing the last factor.