Commit 4b8d3549 authored by Francesco Montanari's avatar Francesco Montanari

Fix Simpson's rule

* cosmo.el (cosmo--trapzd, cosmo-qsimp): Check Simpson's rule for
  convergence and stop if the desired accuracy is reached. The
  previous algorithm did't perform any convergence check. At low
  redshift z~1 numerical integrals are now about one order of
  magnitude faster, at large redshift z~1 one order of magnitude
  slower (but it is guaranteed to reach the desired accuracy).

* cosmo.el (cosmo, cosmo-int-prec, cosmo-int-maxsteps): Define
  customization group to allow users to set numerical integrals
  precision variable through the `customize' command.

* README.org: Explain how to customize variables.
parent 3da5fec4
......@@ -138,6 +138,10 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
- cosmo-transverse-comoving-distance :: Display transverse
comoving distance in mini-buffer.
Use the command =M-x customize= to set variables related to
internal implementation details (e.g., the precision of
numerical integrals), under the =Cosmo= group.
* Support
Bugs and issues are tracked through the [[https://gitlab.com/montanari/cosmo-el][git repository page]]. Please
......
......@@ -61,26 +61,39 @@
;; In priority order:
;;
;; - Allow users to load a customized file with the cosmo-pedia
;; command (or describe how to do it from ~/.emacs).
;;
;; - Simpson's rule performs well for standard cosmologies and for z <
;; 1000. If non-standard cosmologies or very large redshifts are
;; required, more steps may be required. They should be set
;; automatically (the user should only pass a max number of steps)
;; or an error should be shown if the result does not converge.
;;
;; - Refactor tests, now the code is too cumbersome and repetitive.
;; - Refactor tests, now they are too cumbersome and repetitive.
;;
;; - Add all quantities from Hogg 1999.
;;
;; - Suggest default parameters when reading them with the related
;; command; set the to default values if none is entered.
;;
;; - Allow users to load a customized file with the cosmo-pedia
;; command (or describe how to do it from ~/.emacs).
;;
;; - (Consider the following only if performance becomes critical.) At
;; a fixed redshift, only cosmo-get-los-comoving-distance perform
;; the integral. Other distances are defined in terms of this
;; one. However, each time that other distances are required this
;; distance is called once again. This is not necessary, it can be
;; called just once for a given redshift (or use memoization).
;;; Code:
;;; Global variables.
(defgroup cosmo nil
"Cosmological calculator."
:group 'convenience)
(defcustom cosmo-int-prec 1e-3
"Fractional accuracy for numerical integrals."
:group 'cosmo)
(defcustom cosmo-int-maxsteps 20
"Maximum number of steps in the numerical integral algorithm."
:group 'cosmo)
(defvar cosmo-pedia-text
;; Text appearing in the *Cosmopedia* buffer.
"(`q` to quite)
......@@ -168,84 +181,59 @@ Important constants
"Hyperbolic sine of real arguments X."
(* 0.5 (- (exp x) (exp (- x)))))
(defun cosmo-simps (func a b &optional n)
"Simpson rule.
Integrate a function FUNC of one argument from A to B in N steps.
The values A and B will be considered as float.
Example:
\(cosmo-simps #'\(lambda \(x\) x\) 0.0 1.0\)"
(let ((a (float a)) ; Extremes must be floats.
(b (float b))
(n (or n 50)))
(loop with h = (/ (- b a) n)
with sum1 = (funcall func (+ a (/ h 2)))
with sum2 = 0
for i from 1 below n
do (incf sum1 (funcall func (+ a (* h i) (/ h 2))))
do (incf sum2 (funcall func (+ a (* h i))))
finally (return (* (/ h 6)
(+ (funcall func a)
(funcall func b)
(* 4 sum1)
(* 2 sum2)))))))
;; (defun cosmo-trapzd (func lo hi sum niter)
;; "Extended trapezoidal rule to integrate FUNC from LO to HI.
;; Compute the refinement to an old SUM result, previous to the
;; current NITER stage of refinement. The function must be called
;; for NITER = 1, 2, 3, ... in order. Boundaries LO and HI are
;; always converted to floats."
;; (let ((lo (float lo)) ; Always assume float boundaries.
;; (hi (float hi)))
;; ;; NITER=1 corresponds simply to integrate(func(x), x, lo, hi)
;; (if (= niter 1)
;; (* 0.5 (- hi lo) (+ (funcall func lo) (funcall func hi)))
;; ;; Subsequent calls NITER = 2, 3, ... (in order) improve the
;; ;; accuracy by adding 2^(NITER-2) interior points.
;; (let* ((iter (expt 2 (- niter 2)))
;; (step (/ (- hi lo) iter))
;; (xval (+ lo (* 0.5 step)))
;; (sumf 0.0)
;; (j 0))
;; (while (< j iter)
;; (setq sumf (+ sumf (funcall func xval)))
;; (setq xval (+ xval step))
;; (setq j (1+ j)))
;; (* 0.5 (+ sum (* (- hi lo) (/ sumf iter))))))))
;; (defun cosmo-qsimp (func lo hi &optional eps jmax)
;; "Simpson's integration.
;; Argument FUNC integrand.
;; Argument LO lower bound.
;; Argument HI upper bound.
;; Optional argument EPS fractional accuracy.
;; Optional argument JMAX maximum number of steps.
;; "
;; (let ((eps (or eps 1e-6)) ; Set defaults.
;; (jmax (or jmax 20))
;; (nsum 0.0) ; N-th iteration sum.
;; (nsumt 0.0)
;; (osum 0.0) ; Old sum.
;; (osumt 0.0)
;; (converged nil)
;; (j 1))
;; (while (and (<= j jmax) (not converged))
;; (setq nsumt (cosmo-trapzd func lo hi nsum j))
;; (setq nsum (/ (- (* 4.0 nsumt) osumt) 3.0))
;; (if (and (> j 5) ; Avoid spurious early convergence.
;; (or (< (abs (- nsum osum)) (* eps (abs osum)))
;; (and (= nsum 0.0) (= osum 0.0))))
;; (setq converged t))
;; (setq osum nsum)
;; (setq osumt nsumt)
;; (setq j (1+ j)))
;; (if converged
;; nsum
;; (error "Error: The integral did not converge. Try to
;; increase the number of steps"))))
(defun cosmo--trapzd (func lo hi sum niter)
"Extended trapezoidal rule to integrate FUNC from LO to HI.
Compute the refinement to an old SUM result, previous to the
current NITER stage of refinement. The function must be called
for NITER = 1, 2, 3, ... in order. Boundaries LO and HI are
always converted to floats."
(let ((lo (float lo)) ; Always assume float boundaries.
(hi (float hi)))
;; NITER=1 corresponds simply to integrate(func(x), x, lo, hi)
(if (= niter 1)
(* 0.5 (- hi lo) (+ (funcall func lo) (funcall func hi)))
;; Subsequent calls NITER = 2, 3, ... (in order) improve the
;; accuracy by adding 2^(NITER-2) interior points.
(let* ((iter (expt 2 (- niter 2)))
(step (/ (- hi lo) iter))
(xval (+ lo (* 0.5 step)))
(sumf 0.0)
(j 0))
(while (< j iter)
(setq sumf (+ sumf (funcall func xval)))
(setq xval (+ xval step))
(setq j (1+ j)))
(* 0.5 (+ sum (* (- hi lo) (/ sumf iter))))))))
(defun cosmo-qsimp (func lo hi &optional eps jmax)
"Simpson's integration.
Argument FUNC integrand.
Argument LO lower bound.
Argument HI upper bound.
Optional argument EPS fractional accuracy.
Optional argument JMAX maximum number of steps."
(let ((eps (or eps 1e-6)) ; Set defaults.
(jmax (or jmax 20))
(nsum 0.0) ; N-th iteration sum.
(nsumt 0.0)
(osum 0.0) ; Old sum.
(osumt 0.0)
(converged nil)
(j 1))
(while (and (<= j jmax) (not converged))
(setq nsumt (cosmo--trapzd func lo hi nsum j))
(setq nsum (/ (- (* 4.0 nsumt) osumt) 3.0))
(if (and (> j 5) ; Avoid spurious early convergence.
(or (< (abs (- nsum osum)) (* eps (abs osum)))
(and (= nsum 0.0) (= osum 0.0))))
(setq converged t))
(setq osum nsum)
(setq osumt nsumt)
(setq j (1+ j)))
(if converged
nsum
(error "Error: the integral did not converge, try to
increase the number of steps"))))
;;; Compute cosmological functions.
......@@ -306,8 +294,8 @@ Example:
(defun cosmo-get-los-comoving-distance (redshift)
"Line-of-sight comoving distance [Mpc] for Lambda-CDM at a given REDSHIFT."
(let ((DH (cosmo-get-hubble-distance))
(int (cosmo-simps #'cosmo-inv-efunc 0.0 redshift 1000)))
;; (int (cosmo-qsimp #'cosmo-inv-efunc 0.0 redshift 1e-3 20)))
(int (cosmo-qsimp #'cosmo-inv-efunc 0.0 redshift
cosmo-int-prec cosmo-int-maxsteps)))
(* DH int)))
(defun cosmo-los-comoving-distance ()
......@@ -444,9 +432,6 @@ Argument ANGULAR-DIST angular diameter distance at given redshift."
los-dist transverse-dist luminosity-dist
angular-dist))))
;;; Cosmopedia
(defun cosmo-pedia ()
"Display a reference to basic cosmological definitions."
(interactive)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment