1 | r""" |
---|
2 | Show numerical precision of $3 j_1(x)/x$ vs. its Taylor series expansion. |
---|
3 | |
---|
4 | The choice of the number of terms in the series and the cutoff value for |
---|
5 | switching between series and direct calculation depends on the numeric |
---|
6 | precision. |
---|
7 | |
---|
8 | Point 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 | |
---|
13 | Point where Taylor series reaches machine precision (eps), where taylor |
---|
14 | series 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 | |
---|
24 | Note: relative error on single precision starts increase on the direct |
---|
25 | method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3. This should be |
---|
26 | safe for the sans range, with objects of 100 nm supported to a q of 0.1 |
---|
27 | while maintaining 5 digits of precision. For usans/sesans, the objects |
---|
28 | are larger but the q is smaller, so again it should be fine. |
---|
29 | """ |
---|
30 | |
---|
31 | import numpy as np |
---|
32 | from sympy.mpmath import mp |
---|
33 | #import matplotlib; matplotlib.use('TkAgg') |
---|
34 | import pylab |
---|
35 | |
---|
36 | mp.dps = 150 # number of digits to use in estimating true value |
---|
37 | |
---|
38 | SHOW_DIFF = True # True if show diff rather than function value |
---|
39 | LINEAR_X = False # True if q is linearly spaced instead of log spaced |
---|
40 | |
---|
41 | def mp_j1c(vec): |
---|
42 | """ |
---|
43 | Direct calculation using sympy multiprecision library. |
---|
44 | """ |
---|
45 | return [_mp_j1c(mp.mpf(x)) for x in vec] |
---|
46 | |
---|
47 | def _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 | |
---|
53 | def 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 | |
---|
60 | def 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 | |
---|
78 | def 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 | |
---|
88 | def 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 | |
---|
101 | def 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 | |
---|
125 | def 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 | |
---|
142 | if __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() |
---|