I was having problems with the accuracy of floats in Python. I need high accuracy because I want to use explicitly written spherical bessel functions J_n (x), which deviate (especially for n>5) from their theoretical values at low x values if numpy floats are used (15 precise digits).
I have tried many options, especially from mpmath and sympy, in order to keep more precise numbers. I had problems when combining the accuracy of mpmath inside the functions with numpy arrays, until I knew there was the function numpy.vectorize. Finally I got this solution to my initial problem:
import time
% matplotlib qt
import scipy
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from sympy import *
from mpmath import *
mp.dps=100
#explicit inaccurate
def bessel6_expi(z):
return -((z**6-210*z**4+4725*z**2-10395)*np.sin(z)+(21*z**5-1260*z**3+10395*z)*np.cos(z))/z**7
#explicit inaccurate 1, computation time increases, a bit less inaccuracy
def bessel6_exp1(z):
def bv(z):
return -((z**6-210*z**4+4725*z**2-10395)*mp.sin(z)+(21*z**5-1260*z**3+10395*z)*mp.cos(z))/z**7
bvec=np.vectorize(bv)
return bvec(z)
#explicit accurate 2, computation time increases markedly, accurate
def bessel6_exp2(z):
def bv(z):
return -((mpf(z)**mpf(6)-mpf(210)*mpf(z)**mpf(4)+mpf(4725)*mpf(z)**mpf(2)-mpf(10395))*mp.sin(mpf(z))+(mpf(21)*mpf(z)**mpf(5)-mpf(1260)*mpf(z)**mpf(3)+mpf(10395)*mpf(z))*mp.cos(mpf(z)))/mpf(z)**mpf(7)
bvec=np.vectorize(bv)
return bvec(z)
#explicit accurate 3, computation time increases markedly, accurate
def bessel6_exp3(z):
def bv(z):
return -((mpf(z)**6-210*mpf(z)**4+4725*mpf(z)**2-10395)*mp.sin(mpf(z))+(21*mpf(z)**5-1260*mpf(z)**3+10395*mpf(z))*mp.cos(mpf(z)))/mpf(z)**7
bvec=np.vectorize(bv)
return bvec(z)
#implemented in scipy, accurate, fast
def bessel6_imp(z):
def bv(z):
return scipy.special.sph_jn(6,(z))[0][6]
bvec=np.vectorize(bv)
return bvec(z)
a=np.arange(0.0001,17,0.0001)
plt.figure()
start = time.time()
plt.plot(a,bessel6_expi(a),'b',lw=1,label='expi')
end = time.time()
print(end - start)
start = time.time()
plt.plot(a,bessel6_exp1(a),'m',lw=1,label='exp1')
end = time.time()
print(end - start)
start = time.time()
plt.plot(a,bessel6_exp2(a),'c',lw=3,label='exp2')
end = time.time()
print(end - start)
start = time.time()
plt.plot(a,bessel6_exp2(a),'y',lw=5,linestyle='--',label='exp3')
end = time.time()
print(end - start)
start = time.time()
plt.plot(a,bessel6_imp(a),'r',lw=1,label='imp')
end = time.time()
print(end - start)
plt.ylim(-0.5/10**7,2.5/10**7)
plt.xlim(0,2.0)
plt.legend()
plt.show()
The problem I have now is that just for plotting the explicit, accurate ones, it takes quite a long time (about 31 times slower than the scipy function for mp.dps=100). Smaller dps do not make these processes much faster, even with mp.dps=15, they are still 26 times slower. Is there a way to make this faster?
Note that the loss of accuracy you observe near zero comes from the fact that you are subtracting two nearly equal terms both of the form 10395 z^-6 + O(z^-4). As the true value is 1/135135 z^6 + O(z^8) you will lose a factor of ~1.4 x 10^9 z^-12 in accuracy. So if you want to calculate the value at z=0.01 to, say, 7 decimals you need to start with >40 decimals precision.
The solution is of course to avoid this cancellation. A straight-forward way of achieving this is to compute the power series around 0.
You could use sympy to obtain the power series:
>>> z = sympy.Symbol('z')
>>> f = -((z**6-210*z**4+4725*z**2-10395)*sympy.sin(z)+(21*z**5-1260*z**3+10395*z)*sympy.cos(z))/z**7
>>> f.nseries(n=20)
z**6/135135 - z**8/4054050 + z**10/275675400 - z**12/31426995600 + z**14/5279735260800 - z**16/1214339109984000 + z**18/364301732995200000 + O(z**20)
For small z a small number of terms appear to be enough for good accuracy.
>>> ply = f.nseries(n=20).removeO().as_poly()
>>> float(ply.subs(z, 0.1))
7.397541093587708e-12
You can export the coefficients for use with numpy.
>>> monoms = np.array(ply.monoms(), dtype=int).ravel()
>>> coeffs = np.array(ply.coeffs(), dtype=float)
>>>
>>> (np.linspace(-0.1, 0.1, 21)[:, None]**monoms * coeffs).sum(axis=1)
array([7.39754109e-12, 3.93160564e-12, 1.93945374e-12, 8.70461282e-13,
3.45213317e-13, 1.15615481e-13, 3.03088138e-14, 5.39444356e-15,
4.73594159e-16, 7.39998273e-18, 0.00000000e+00, 7.39998273e-18,
4.73594159e-16, 5.39444356e-15, 3.03088138e-14, 1.15615481e-13,
3.45213317e-13, 8.70461282e-13, 1.93945374e-12, 3.93160564e-12,
7.39754109e-12])