1

Python の float の精度に問題がありました。明示的に記述された球面ベッセル関数 J_n (x) を使用したいので、高精度が必要です。これは、numpyfloat が使用されている場合 (特に n>5 の場合)、低い x 値で理論値から逸脱します (正確な 15 桁)。

より正確な数値を維持するために、特にmpmathとから多くのオプションを試しました。function があることがわかるまで、関数内の精度を配列sympyと組み合わせるときに問題がありました。最後に、最初の問題に対するこの解決策を得ました。mpmathnumpynumpy.vectorize

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()

私が今抱えている問題は、明示的で正確なものをプロットするだけで、かなり長い時間がかかることです ( の scipy 関数よりも約 31 倍遅いmp.dps=100)。小さくdpsしても、これらのプロセスはそれほど速くはなりませんmp.dps=15。. これをより速くする方法はありますか？

4

1 に答える 1