| 1 | r""" |
|---|
| 2 | Show numerical precision of $sin(x)/x$. |
|---|
| 3 | """ |
|---|
| 4 | |
|---|
| 5 | import numpy as np |
|---|
| 6 | from sympy.mpmath import mp |
|---|
| 7 | #import matplotlib; matplotlib.use('TkAgg') |
|---|
| 8 | import pylab |
|---|
| 9 | |
|---|
| 10 | mp.dps = 150 # number of digits to use in estimating true value |
|---|
| 11 | |
|---|
| 12 | SHOW_DIFF = True # True if show diff rather than function value |
|---|
| 13 | LINEAR_X = False # True if q is linearly spaced instead of log spaced |
|---|
| 14 | |
|---|
| 15 | def mp_sinc(vec): |
|---|
| 16 | """ |
|---|
| 17 | Direct calculation using sympy multiprecision library. |
|---|
| 18 | """ |
|---|
| 19 | return [_mp_sinc(mp.mpf(x)) for x in vec] |
|---|
| 20 | |
|---|
| 21 | def _mp_sinc(x): |
|---|
| 22 | """ |
|---|
| 23 | Helper funciton for mp_j1c |
|---|
| 24 | """ |
|---|
| 25 | return mp.sin(x)/x |
|---|
| 26 | |
|---|
| 27 | def np_sinc(x, dtype): |
|---|
| 28 | """ |
|---|
| 29 | Direct calculation using scipy. |
|---|
| 30 | """ |
|---|
| 31 | x = np.asarray(x, dtype) |
|---|
| 32 | return np.sin(x)/x |
|---|
| 33 | #return np.asarray(np.sin(np.double(x))/np.double(x),dtype) |
|---|
| 34 | |
|---|
| 35 | def plotdiff(x, target, actual, label): |
|---|
| 36 | """ |
|---|
| 37 | Plot the computed value. |
|---|
| 38 | |
|---|
| 39 | Use relative error if SHOW_DIFF, otherwise just plot the value directly. |
|---|
| 40 | """ |
|---|
| 41 | if SHOW_DIFF: |
|---|
| 42 | err = np.clip(abs((target-actual)/target), 0, 1) |
|---|
| 43 | pylab.loglog(x, err, '-', label=label) |
|---|
| 44 | else: |
|---|
| 45 | limits = np.min(target), np.max(target) |
|---|
| 46 | pylab.semilogx(x, np.clip(actual,*limits), '-', label=label) |
|---|
| 47 | |
|---|
| 48 | def compare(x, precision): |
|---|
| 49 | r""" |
|---|
| 50 | Compare the different computation methods using the given precision. |
|---|
| 51 | """ |
|---|
| 52 | target = np.asarray(mp_sinc(x), 'double') |
|---|
| 53 | direct = np_sinc(x, precision) |
|---|
| 54 | plotdiff(x, target, direct, 'direct '+precision) |
|---|
| 55 | pylab.xlabel("qr (1/Ang)") |
|---|
| 56 | if SHOW_DIFF: |
|---|
| 57 | pylab.ylabel("relative error") |
|---|
| 58 | else: |
|---|
| 59 | pylab.ylabel("sin(x)/x") |
|---|
| 60 | pylab.semilogx(x, target, '-', label="true value") |
|---|
| 61 | if LINEAR_X: |
|---|
| 62 | pylab.xscale('linear') |
|---|
| 63 | |
|---|
| 64 | def main(): |
|---|
| 65 | r""" |
|---|
| 66 | Compare accuracy of different methods for computing $3 j_1(x)/x$. |
|---|
| 67 | :return: |
|---|
| 68 | """ |
|---|
| 69 | if LINEAR_X: |
|---|
| 70 | qr = np.linspace(1,1000,2000) |
|---|
| 71 | else: |
|---|
| 72 | qr = np.logspace(-3,5,400) |
|---|
| 73 | pylab.subplot(121) |
|---|
| 74 | compare(qr, 'single') |
|---|
| 75 | pylab.legend(loc='best') |
|---|
| 76 | pylab.subplot(122) |
|---|
| 77 | compare(qr, 'double') |
|---|
| 78 | pylab.legend(loc='best') |
|---|
| 79 | pylab.suptitle('sin(x)/x') |
|---|
| 80 | |
|---|
| 81 | if __name__ == "__main__": |
|---|
| 82 | print "\n".join(str(x) for x in mp_sinc([1e-6,1e-5,1e-4,1e-3])) |
|---|
| 83 | main() |
|---|
| 84 | pylab.show() |
|---|
