# 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@…>, 8 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.