source: sasmodels/explore/sph_j1c.py @ e6f1410

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e6f1410 was e6f1410, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

examine numerical precision of J1c, sinc, j1c and log gamma

  • Property mode set to 100644
File size: 4.4 KB
Line 
1r"""
2Show numerical precision of $3 j_1(x)/x$ vs. its Taylor series expansion.
3
4The choice of the number of terms in the series and the cutoff value for
5switching between series and direct calculation depends on the numeric
6precision.
7
8Point where direct calculation reaches machine precision::
9
10    single machine precision eps 3e-8 at qr=1.1 (see below)
11    double machine precision eps 4e-16 at qr=1.1
12
13Point where Taylor series reaches machine precision (eps), where taylor
14series matches direct calculation (cross) and the error at that point::
15
16    prec   n eps  cross  error
17    single 3 0.28  0.4   6.2e-7
18    single 4 0.68  0.7   2.3e-7
19    single 5 1.18  1.2   7.5e-8
20    double 3 0.01  0.03  2.3e-13
21    double 4 0.06  0.1   3.1e-14
22    double 5 0.16  0.2   5.0e-15
23
24Note: relative error on single precision starts increase on the direct
25method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
26safe for the sans range, with objects of 100 nm supported to a q of 0.1
27while maintaining 5 digits of precision.  For usans/sesans, the objects
28are larger but the q is smaller, so again it should be fine.
29"""
30
31import numpy as np
32from sympy.mpmath import mp
33#import matplotlib; matplotlib.use('TkAgg')
34import pylab
35
36mp.dps = 150 # number of digits to use in estimating true value
37
38SHOW_DIFF = True  # True if show diff rather than function value
39LINEAR_X = False  # True if q is linearly spaced instead of log spaced
40
41def mp_j1c(vec):
42    """
43    Direct calculation using sympy multiprecision library.
44    """
45    return [_mp_j1c(mp.mpf(x)) for x in vec]
46
47def _mp_j1c(x):
48    """
49    Helper funciton for mp_j1c
50    """
51    return mp.mpf(3)*(mp.sin(x)/x - mp.cos(x))/(x*x)
52
53def np_j1c(x, dtype):
54    """
55    Direct calculation using numpy.
56    """
57    x = np.asarray(x, dtype)
58    return np.asarray(3, dtype)*(np.sin(x) - x*np.cos(x))/(x*x*x)
59
60def lin_j1c(x, dtype):
61    r"""
62    Direct calculation using linear combination of sin/cos.
63
64    Use the following trig identity:
65
66    .. math::
67
68        a \sin(x) + b \cos(x) = c \sin(x + \phi)
69
70    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
71    numerator $\sin(x) - x\cos(x)$.
72    """
73    x = np.asarray(x, dtype)
74    c = np.sqrt(np.asarray(1,dtype) + x*x)
75    phi = np.arctan2(-np.asarray(x,dtype),np.asarray(1,dtype))
76    return np.asarray(3, dtype)*(c*np.sin(x+phi))/(x*x*x)
77
78def taylor_j1c(x, dtype, n):
79    """
80    Calculation using taylor series.
81    """
82    # Generate coefficients using the precision of the target value.
83    cinv = [3991680, -45360, 840, -30, 3]
84    x = np.asarray(x, dtype)
85    p = np.asarray(3, dtype)/np.array(cinv, dtype)
86    return np.polyval(p[-n:], x*x)
87
88def plotdiff(x, target, actual, label):
89    """
90    Plot the computed value.
91
92    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
93    """
94    if SHOW_DIFF:
95        err = np.clip(abs((target-actual)/target), 0, 1)
96        pylab.loglog(x, err, '-', label=label)
97    else:
98        limits = np.min(target), np.max(target)
99        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
100
101def compare(x, precision):
102    r"""
103    Compare the different computation methods using the given precision.
104    """
105    target = np.asarray(mp_j1c(x), 'double')
106    direct = np_j1c(x, precision)
107    comb = lin_j1c(x, precision)
108    taylor3 = taylor_j1c(x, precision, 3)
109    taylor4 = taylor_j1c(x, precision, 4)
110    taylor5 = taylor_j1c(x, precision, 5)
111    plotdiff(x, target, direct, 'direct '+precision)
112    #plotdiff(x, target, comb, 'c sin(x+phi) '+precision)
113    plotdiff(x, target, taylor3, 'taylor 3 '+precision)
114    plotdiff(x, target, taylor4, 'taylor 4 '+precision)
115    plotdiff(x, target, taylor5, 'taylor 5 '+precision)
116    pylab.xlabel("qr (1/Ang)")
117    if SHOW_DIFF:
118        pylab.ylabel("relative error")
119    else:
120        pylab.ylabel("3 j1(x)/x")
121        pylab.semilogx(x, target,  '-', label="true value")
122    if LINEAR_X:
123        pylab.xscale('linear')
124
125def main():
126    r"""
127    Compare accuracy of different methods for computing $3 j_1(x)/x$.
128    :return:
129    """
130    if LINEAR_X:
131        qr = np.linspace(1,1000,2000)
132    else:
133        qr = np.logspace(-3,5,400)
134    pylab.subplot(121)
135    compare(qr, 'single')
136    pylab.legend(loc='best')
137    pylab.subplot(122)
138    compare(qr, 'double')
139    pylab.legend(loc='best')
140    pylab.suptitle('3 j1(x)/x')
141
142if __name__ == "__main__":
143    #print "\n".join(str(x) for x in mp_j1c([1e-6,1e-5,1e-4,1e-3]))
144    main()
145    pylab.show()
Note: See TracBrowser for help on using the repository browser.