;; -*- Mode:Lisp; Syntax:ANSI-Common-LISP; Coding:us-ascii-unix; fill-column:158 -*- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;; @file use-intg.lisp ;; @author Mitch Richling ;; @brief Numerical integration of univariate real functions.@EOL ;; @std Common Lisp ;; @see tst-intg.lisp ;; @copyright ;; @parblock ;; Copyright (c) 1997,1998,2004,2008,2013,2015, Mitchell Jay Richling All rights reserved. ;; ;; Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: ;; ;; 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer. ;; ;; 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation ;; and/or other materials provided with the distribution. ;; ;; 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software ;; without specific prior written permission. ;; ;; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ;; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE ;; LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS ;; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT ;; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH ;; DAMAGE. ;; @endparblock ;; @todo Integrate data instead of a function.@EOL@EOL ;; @todo Arc length of f:R->R^n on an interval [a,b] via various base integration rules.@EOL@EOL ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defpackage :MJR_INTG (:USE :COMMON-LISP :MJR_NUMU :MJR_VVEC :MJR_UTIL :MJR_PRNG :MJR_EPS :MJR_MXP) (:DOCUMENTATION "Brief: Numerical integration of univariate real functions.;") (:EXPORT #:mjr_intg_help #:mjr_intg_simple-rect-left #:mjr_intg_simple-rect-right ;; fixed-order simple Riemann rules #:mjr_intg_simple-rect-mid #:mjr_intg_simple-milne ;; fixed-order simple Open NC rules #:mjr_intg_simple-trap #:mjr_intg_simple-simpson ;; fixed-order simple Closed NC rules #:mjr_intg_simple-simpson-3/8 #:mjr_intg_simple-boole ;; fixed-order simple Closed NC rules #:mjr_intg_simple-newton-cotes ;; variable-order simple All NC rules #:mjr_intg_simple-gauss-legendre ;; variable-order simple Gauss Legendre #:mjr_intg_simple-monte-carlo ;; variable-order simple Monte Carlo #:mjr_intg_simple-gauss-kronrod ;; variable-order simple embedded GK #:mjr_intg_composite-trap #:mjr_intg_composite-simpson ;; variable-order Composite NC rules #:mjr_intg_glb-adp-composite-romberg #:mjr_intg_glb-adp-composite-trapezoidal ;; Global Adaptive Composite Closed NC rules #:mjr_intg_loc-adp-dnc-trapezoidal #:mjr_intg_loc-adp-dnc-simpson ;; Local Adaptive Div&Conq Closed NC rules ;; Experimental #:mjr_intg_gbl-adp-factory-order ;; Global Adaptive Order various (adapter) #:mjr_intg_divide-and-conquer ;; Local Adaptive Div&Conq various (adapter) ;; Not exported (backend use) ;;#:mjr_intg_glb-adp-composite-trap-or-romberg )) (in-package :MJR_INTG) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_help () "Help for MJR_INTG: Single variable definite integrals computed over an interval The functions: mjr_intg_simple-newton-cotes .............. Good for smooth functions well approximated by a polynomial of degree < order of method mjr_intg_simple-gauss-legendre ............ Fastest for smooth functions well approximated by a polynomial. mjr_intg_glb-adp-composite-romberg ........ Good choice for well behaved functions -- $C^\infty$ mjr_intg_glb-adp-composite-trapezoidal .... Good choice for moderately ill-behaved functions. For example if the derivative fails to exist at a few points, but the variation is generally consistent across interval of integration. mjr_intg_loc-adp-dnc-trapezoidal .......... Good choice for functions with a widely changing second derivative or strong variation that break mjr_intg_glb-adp-composite-trapezoidal. Also good when a particular mesh must be refined. mjr_intg_loc-adp-dnc-simpson .............. Like mjr_intg_loc-adp-dnc-trapezoidal, but uses simpson's rule. mjr_intg_composite-trap ................... Fixed step size composite trapezoidal rule. mjr_intg_composite-simpson ................ Fixed step size composite simpson rule. This is a pathological example CL> (mjr_intrp_poly #(0 1/4 1/2 3/4 1) #(0 1/2 1/2 0 1)) #(64/3 -32 32/3 1 0) CL> (float (mjr_poly_integrate #(64/3 -32 32/3 1 0) 0 1)) 0.32222223 CL> (float (mjr_intg_glb-adp-composite-trapezoidal (lambda (x) (mjr_poly_eval #(64/3 -32 32/3 1 0) x)) :start 0 :end 1 :good-err-stop 1)) 0.5 CL> (float (mjr_intg_glb-adp-composite-trapezoidal (lambda (x) (mjr_poly_eval #(64/3 -32 32/3 1 0) x)) :start 0 :end 1 :good-err-stop 2)) 0.3222256 CL> (float (mjr_intg_loc-adp-dnc-trapezoidal (lambda (x) (mjr_poly_eval #(64/3 -32 32/3 1 0) x)) :start 0 :end 1 :len 2)) 0.5 CL> (float (mjr_intg_loc-adp-dnc-trapezoidal (lambda (x) (mjr_poly_eval #(64/3 -32 32/3 1 0) x)) :start 0 :end 1 :len 4)) 0.3222203" (documentation 'mjr_intg_help 'function)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-monte-carlo (fun &key start end (order 10000) show-progress) "Compute the definite integral of FUN on the given interval using the Monte Carlo method." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (let ((wid (- end start)) (fsum 0) (spi (if (> order 30) (truncate order 30) 1))) (loop for ect from 1 upto order for x = (mjr_prng_float-cc start end) finally (return (* wid fsum (/ order))) do (incf fsum (funcall fun x)) do (if (and show-progress (zerop (mod ect spi))) (format 't "~7d :: ~15f ~%" ect (* wid fsum (/ ect)))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_composite-trap (fun &key start end step len) "Compute the definite integral of FUN using the composite Trapezoidal Rule. START, END, STEP, LEN are processed by MJR_VVEC_NORMALIZE-KW-VVT-ASEQ" (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (multiple-value-bind (start end step len) (mjr_util_get-kwarg-vals '(:start :end :step :len) (mjr_vvec_normalize-kw-vvt-aseq (mjr_util_strip-nil-val-kwarg (list :start start :end end :step step :len len)))) (cond ((< len 2) (error "mjr_intg_trap: :LEN must be greater than 1!")) ((> start end) (error "mjr_intg_trap: :START must be numerically less than :END!"))) (* 1/2 (* step (+ (funcall fun start) (funcall fun end) (* 2 (mjr_numu_sum :seq-fun fun :start (+ start step) :end (- end step) :step step :len (- len 2))))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_composite-simpson (fun &key start end step len) "Compute the definite integral of FUN using the composite Simpson Rule. START, END, STEP, LEN are processed by MJR_VVEC_NORMALIZE-ASEQ" (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (multiple-value-bind (start end step len) (mjr_util_get-kwarg-vals '(:start :end :step :len) (mjr_vvec_normalize-kw-vvt-aseq (mjr_util_strip-nil-val-kwarg (list :start start :end end :step step :len len)))) (cond ((< len 2) (error "mjr_intg_simpson: :LEN must be greater than 5!")) ((evenp len) (error "mjr_intg_simpson: :LEN must be odd!")) ((> start end) (error "mjr_intg_simpson: :START must be numerically less than :END!"))) (* 1/3 (* step (+ (funcall fun start) (funcall fun end) (loop for x = (+ start step) then (+ x step) for i from 1 upto (- len 2) sum (* (if (oddp i) 4 2) (funcall fun x))))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_loc-adp-dnc-simpson (fun &key points start end step len seq-next min-width (the-err 1e-5) max-evals show-progress) "Compute the definite integral of FUN between A and B using an adaptive Simpson's rule. Each sub-interval of the given partition is recursively bisected until the error estimate on each interval is below THE-ERR, the number of function evaluations grows beyond MAX-EVALS, or the intervals become smaller than MIN-WIDTH. If only :START and :END are given of the partition, then :LEN will be set to 2. If provided, :POINTS must be a vector or a non-NIL list. If the error goal is satisfied, then the return is the integral approximation and the number of function evaluations. If any of the limits are violated on any part of the interval, then the return will be nil, the integral approximation, the number of function evaluations required, the total length of the sub-intervals upon which no limits were violated during the computation, and the number of intervals that had a violated limit. References: Kythe & Schaferkotter (2005); Handbook of Computational Methods for Integration; pp 89-94" (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (flet ((simpson (a b fa fm fb) (* (/ (- b a) 6) (+ fa (* 4 fm) fb)))) (let ((ect 0) (len (if (and start end (not (or points step len seq-next))) 2 len))) (let ((ilist (loop for x0 = nil then x1 for x1 in (mjr_vvec_to-list (mjr_util_strip-nil-val-kwarg (list :points points :start start :end end :step step :len len :map-fun seq-next))) for f0 = nil then f1 for f1 = (funcall fun x1) do (incf ect) when x0 collect (let* ((m (/ (+ x0 x1) 2)) (fm (funcall fun m))) (incf ect) (list x0 m x1 f0 fm f1 (simpson x0 f1 f0 fm f1) (/ (* the-err (- x1 x0)) (- end start)))))) (badt 0) (cerw 0) (civ 0)) (loop for (a m b fa fm fb vo eps) = (pop ilist) for iwid = (and a (- b a)) while iwid finally (return (if (zerop badt) (values civ ect) (values nil civ ect cerw badt))) do (if (and min-width (< iwid min-width)) (progn (if show-progress (format 't "~4d :: ~15f TLW ~15f ~%" ect civ iwid)) (incf badt) (incf civ vo)) (if (and max-evals (> ect max-evals)) (progn (if show-progress (format 't "~4d :: ~15f TLF ~15f~%" ect civ iwid)) (incf badt) (incf civ vo)) (let* ((ml (/ (+ a m) 2)) (mr (/ (+ m b) 2)) (fml (funcall fun ml)) (fmr (funcall fun mr)) (eps2 (/ eps 2)) (vl (simpson a m fa fml fm)) (vr (simpson m b fm fmr fb)) (err (abs (- vo vl vr)))) (incf ect 2) (if (< err eps2) (progn (if show-progress (format 't "~4d :: ~15f TER ~15f ~%" ect civ iwid)) (incf cerw iwid) (incf civ (+ vl vr))) (progn (if show-progress (format 't "~4d :: ~15f CER ~15f ~15f ~%" ect civ iwid err)) (setq ilist (append ilist (list (list a ml m fa fml fm vl eps2) (list m mr b fm fmr fb vr eps2))))))))))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_glb-adp-composite-trap-or-romberg (fun &key start end max-evals min-width the-err good-err-stop show-progress do-romberg) "Back end function to do adaptive trapezoidal or adaptive Romberg. See: MJR_INTG_TRAP-ADP and/or MJR_INTG_ROMBERG-ADP" (cond ((not (numberp start)) (error "mjr_intg_glb-adp-composite-trap-or-romberg: Lower limit (START) must be a number!")) ((complexp start) (error "mjr_intg_glb-adp-composite-trap-or-romberg: Lower limit (START) must be a real number!")) ((not (numberp end)) (error "mjr_intg_glb-adp-composite-trap-or-romberg: Upper limit (END) must be a number!")) ((complexp end) (error "mjr_intg_glb-adp-composite-trap-or-romberg: Upper limit (END) must be a real number!")) ((> start end) (error "mjr_intg_glb-adp-composite-trap-or-romberg: Lower limit (START) must be numerically less than upper limit (END)!")) ((and max-evals (not (integerp max-evals))) (error "mjr_intg_glb-adp-composite-trap-or-romberg: MAX-EVALS must be an integer!")) ((and max-evals (< max-evals 3)) (error "mjr_intg_glb-adp-composite-trap-or-romberg: MAX-EVALS must be greater than 2!")) ((not (integerp good-err-stop)) (error "mjr_intg_glb-adp-composite-trap-or-romberg: GOOD-ERR-STOP must be an integer!")) ((< good-err-stop 1) (error "mjr_intg_glb-adp-composite-trap-or-romberg: GOOD-ERR-STOP must be 1 or greater!")) ((not (numberp the-err)) (error "mjr_intg_glb-adp-composite-trap-or-romberg: THE-ERR must be a number!")) ((complexp the-err) (error "mjr_intg_glb-adp-composite-trap-or-romberg: THE-ERR must be a real number!")) ((<= the-err 0) (error "mjr_intg_glb-adp-composite-trap-or-romberg: THE-ERR must be a positive number!"))) (loop with good-c = 0 with lim-del = (- end start) for n from 2 for np = 1 then (* np 2) for ect = 2 then (+ ect np) for hn = (/ lim-del 2 np) for old-trap = (/ (* lim-del (+ (funcall fun start) (funcall fun end))) 2) then new-trap for new-trap = (+ (/ old-trap 2) (* hn (loop for j from 1 upto np sum (funcall fun (+ start (* hn (- (* 2 j) 1))))))) for old-romb = old-trap then new-romb for rom-list = (list old-trap) then (if do-romberg (loop for m = 1 then (* 4 m) for j from 2 upto n for or in (append (list nil) rom-list) for r = old-trap then (+ r (/ (- r or) (- m 1))) collect r)) for new-romb = (car (last rom-list)) for new-int = (if do-romberg new-romb new-trap) for old-int = (if do-romberg old-romb old-trap) ;;do (format 't "~a~%" rom-list) do (if show-progress (format 't "~4d :: ~15f ~%" ect new-int)) do (if (< (mjr_numu_absdif old-int new-int) the-err) (incf good-c) (setq good-c 0)) do (if (and (< 4 ect) (>= good-c good-err-stop)) (return-from mjr_intg_glb-adp-composite-trap-or-romberg (values new-int ect))) do (if (and min-width (> min-width hn)) (return-from mjr_intg_glb-adp-composite-trap-or-romberg (values nil new-int ect))) do (if (and max-evals (> ect max-evals)) (return-from mjr_intg_glb-adp-composite-trap-or-romberg (values nil new-int ect))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_glb-adp-composite-romberg (fun &key start end max-evals min-width (the-err 1e-5) (good-err-stop 2) show-progress) "Compute the definite integral of FUN between A and B using an adaptive Romberg scheme. The number of steps starts at 2 and is doubled (roughly) until it grows beyond MAX-N or we see GOOD-ERR-STOP consecutive iterations with an integral approximation that changes by less than THE-ERR. If the error goal is satisfied, then the return is the integral approximation and the number of function evaluations. If MAX-N is violated, then nil, the integral approximation, and the number of function evaluations returned. min-width sets the smallest intervals that will be used, and if it is violated then nil, the integral approximation, and the number of function evaluations are returned." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (mjr_intg_glb-adp-composite-trap-or-romberg fun :start start :end end :max-evals max-evals :min-width min-width :the-err the-err :good-err-stop good-err-stop :show-progress show-progress :do-romberg 't))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_glb-adp-composite-trapezoidal (fun &key start end max-evals min-width (the-err 1e-5) (good-err-stop 2) show-progress) "Compute the definite integral of FUN between A and B using adaptive the trapezoidal rule. The number of steps starts at 2 and is doubled (roughly) until it grows beyond MAX-N or we see GOOD-ERR-STOP consecutive iterations with an integral approximation that changes by less than THE-ERR. If the error goal is satisfied, then the return is the integral approximation and the number of function evaluations. If MAX-N is violated, then nil, the integral approximation, and the number of function evaluations returned. min-width sets the smallest intervals that will be used, and if it is violated then nil, the integral approximation, and the number of function evaluations are returned." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (mjr_intg_glb-adp-composite-trap-or-romberg fun :start start :end end :max-evals max-evals :min-width min-width :the-err the-err :good-err-stop good-err-stop :show-progress show-progress :do-romberg nil))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_divide-and-conquer (fun simple-integration-func vvec &rest int-args) "Compute the definite integral of FUN on the given numeric sequence using the given interval integration rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (mjr_vvec_map-sum vvec :pair-fun (lambda (x0 x1 fx0 fx1 i) (declare (ignore fx0 fx1 i)) (apply simple-integration-func fun :start x0 :end x1 int-args))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_loc-adp-dnc-trapezoidal (fun vvec &key min-width max-evals (good-err-stop 2) (the-err 1e-5) suppress-warnings show-progress) "Compute the definite integral of FUN between A and B using a doubly adaptive trapezoidal rule. The initial intervals specified by the virtual vector passed in vvec are recursively bisected until intervals are obtained that are the result of :GOOD-ERR-STOP consecutive iterations with an integral approximation that changes by less than :THE-ERR scaled for the interval's size. If the error goal is satisfied, then the return is the integral approximation and the number of times the function was evaluated. If :MAX-EVALS is violated, the return is NIL and the number of times the function was evaluated. If an interval with an unsatisfactory integral approximation is less than :MIN-WIDTH, it will not be subdivided -- the unsatisfactory approximation will be used and a warning will be printed unless :SUPPRESS-WARNINGS is non-NIL." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (cond ((and max-evals (not (integerp max-evals))) (error "mjr_intg_loc-adp-dnc-trapezoidal: MAX-EVALS must be an integer!")) ((and max-evals (< max-evals 3)) (error "mjr_intg_loc-adp-dnc-trapezoidal: MAX-EVALS must be greater than 2!")) ((not (integerp good-err-stop)) (error "mjr_intg_loc-adp-dnc-trapezoidal: GOOD-ERR-STOP must be an integer!")) ((< good-err-stop 1) (error "mjr_intg_loc-adp-dnc-trapezoidal: GOOD-ERR-STOP must be 1 or greater!")) ((and min-width (not (numberp min-width))) (error "mjr_intg_loc-adp-dnc-trapezoidal: MIN-WIDTH must be a number!")) ((and min-width (complexp min-width)) (error "mjr_intg_loc-adp-dnc-trapezoidal: MIN-WIDTH must be a real number!")) ((and min-width (<= min-width 0)) (error "mjr_intg_loc-adp-dnc-trapezoidal: MIN-WIDTH must be a positive number!")) ((not (numberp the-err)) (error "mjr_intg_loc-adp-dnc-trapezoidal: THE-ERR must be a number!")) ((complexp the-err) (error "mjr_intg_loc-adp-dnc-trapezoidal: THE-ERR must be a real number!")) ((<= the-err 0) (error "mjr_intg_loc-adp-dnc-trapezoidal: THE-ERR must be a positive number!"))) (let* ((points (mjr_vvec_to-vec vvec)) (len (length points))) (cond ((< len 2) (error "mjr_intg_loc-adp-dnc-trapezoidal: :LEN must be greater than 1!"))) (let* ((intervals (mjr_vvec_map-filter-reduce 'list points :point-fun (lambda (x i) (declare (ignore i)) (funcall fun x)) :pair-fun (lambda (x0 x1 fx0 fx1 i) (declare (ignore i)) (list x0 x1 fx0 fx1 0)))) (int-width (- (aref points (1- len)) (aref points 0)))) ;; Refine intervals till we get a good value. (loop with intv = 0 for ect from (1+ len) ;; func evaluation count for (cur-x0 cur-x1 cur-fx0 cur-fx1 cur-gct) = (pop intervals) for cur-sum = (+ cur-fx1 cur-fx0) for cur-wid = (- cur-x1 cur-x0) for cur-itg = (/ (* cur-sum cur-wid) 2) for new-x = (/ (+ cur-x0 cur-x1) 2) for new-mfx = (funcall fun new-x) for new-itg = (/ (* (+ cur-sum (* 2 new-mfx)) cur-wid) 4) do (if show-progress (format 't "~4d :: [~15f,~15f] :: ~15f :: ~15f ~%" ect cur-x0 cur-x1 cur-itg new-itg)) do (if (< (mjr_numu_absdif cur-itg new-itg) (/ (* the-err cur-wid) int-width)) (if (>= (1+ cur-gct) good-err-stop) (incf intv new-itg) (push (list cur-x0 cur-x1 cur-fx0 cur-fx1 (1+ cur-gct)) intervals)) (if (and min-width (> min-width (- new-x cur-x0))) (progn (incf intv new-itg) (or suppress-warnings (warn "mjr_intg_loc-adp-dnc-trapezoidal: Violated MIN-WIDTH constraint!"))) (progn (push (list new-x cur-x1 new-mfx cur-fx1 0) intervals) (push (list cur-x0 new-x cur-fx0 new-mfx 0) intervals)))) do (if (null intervals) (return-from mjr_intg_loc-adp-dnc-trapezoidal (values intv ect))) do (if (and max-evals (> ect max-evals)) (return-from mjr_intg_loc-adp-dnc-trapezoidal (values nil ect)))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-trap (fun &key start end) "Compute the definite integral of FUN between A and B using the Trapezoidal rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (- end start) (/ (+ (funcall fun start) (funcall fun end)) 2)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-rect-left (fun &key start end) "Compute the definite integral of FUN between A and B using the Left Rectangle rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (- end start) (funcall fun start)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-rect-right (fun &key start end) "Compute the definite integral of FUN between A and B using the Right Rectangle rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (- end start) (funcall fun end)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-rect-mid (fun &key start end) "Compute the definite integral of FUN between A and B using the Midpoint Rectangle rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (- end start) (funcall fun (/ (+ start end) 2))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-simpson (fun &key start end) "Compute the definite integral of FUN between A and B using Simpson's rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (/ (- end start) 6) (+ (funcall fun start) (* 4 (funcall fun (/ (+ end start) 2))) (funcall fun end))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-simpson-3/8 (fun &key start end) "Compute the definite integral of FUN between A and B using Simpson's 3/8 rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (* (/ (- end start) 8) (+ (funcall fun start) (* 3 (funcall fun (/ (+ end (* 2 start)) 3))) (* 3 (funcall fun (/ (+ (* 2 end) start) 3))) (funcall fun end))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-boole (fun &key start end) "Compute the definite integral of FUN between A and B using Boole's rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x")) (h (/ (- end start) 4))) (* (/ (- end start) 90) (+ (* 7 (funcall fun start)) (* 32 (funcall fun (+ start h))) (* 12 (funcall fun (+ start (* 2 h)))) (* 32 (funcall fun (+ start (* 3 h)))) (* 7 (funcall fun end)))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-milne (fun &key start end) "Compute the definite integral of FUN between A and B using Milne's rule." (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x")) (h (/ (- end start) 4))) (* (/ (- end start) 3) (+ (* 2 (funcall fun (+ start h))) (- (funcall fun (+ start (* 2 h)))) (* 2 (funcall fun (+ start (* 3 h)))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-newton-cotes (fun &key start end (order 3) (closed 't)) "Newton Cotes Integration (both open and closed) When closed is non-NIL, the order must be in [1,20]. When closed is NIL, the order must be in [2,20]. |-------+-------------+-----------------------------+-------------| | order | open/closed | LISP function | Name | |-------+-------------+-----------------------------+-------------| | 1 | closed | mjr_intg_simple-trap | Trapezoidal | | 2 | closed | mjr_intg_simple-simpson | Simpson | | 3 | closed | mjr_intg_simple-simpson-3/8 | Simpson 3/8 | | 4 | closed | mjr_intg_simple-boole | Boole | |-------+-------------+-----------------------------+-------------| | 2 | open | mjr_intg_simple-mid | Midpoint | | 4 | open | mjr_intg_simple-milne | Milne | |-------+-------------+-----------------------------+-------------|" ;; Accuracy generally improves as order increases, but only up to a point. For example, we can integrate 1/x from 1 to 2 with various orders to illustrate ;; the point: ;; ---------------------------------------------------------------------------------------------------- ;; (loop for o from 1 upto 20 ;; for ic = (if (> o 0) (mjr_intg_simple-newton-cotes (lambda (x) (/ (1+ x))) :start 0d0 :end 1d0 :order o :closed 't)) ;; for io = (if (> o 1) (mjr_intg_simple-newton-cotes (lambda (x) (/ (1+ x))) :start 0d0 :end 1d0 :order o :closed nil)) ;; do (if io ;; (format 't "~5d : ~15,10f ~15,10f ~15,10f ~15,10f~%" o ic (abs (- ic (log 2.0d0))) io (abs (- io (log 2.0d0)))) ;; (format 't "~5d : ~15,10f ~15,10f ~15a ~15a~%" o ic (abs (- ic (log 2.0d0))) "-" "-"))) ;; ---------------------------------------------------------------------------------------------------- ;; 1 : 0.7500000000 0.0568528194 - - ;; 2 : 0.6944444444 0.0012972639 0.6666666667 0.0264805139 ;; 3 : 0.6937500000 0.0006028194 0.6750000000 0.0181471806 ;; 4 : 0.6931746032 0.0000274226 0.6920634921 0.0010836885 ;; 5 : 0.6931630291 0.0000158485 0.6923776455 0.0007695351 ;; 6 : 0.6931480623 0.0000008817 0.6930952381 0.0000519425 ;; 7 : 0.6931477333 0.0000005528 0.6931097433 0.0000374373 ;; 8 : 0.6931472145 0.0000000340 0.6931445039 0.0000026766 ;; 9 : 0.6931472028 0.0000000222 0.6931452353 0.0000019452 ;; 10 : 0.6931471820 0.0000000015 0.6931470368 0.0000001438 ;; 11 : 0.6931471815 0.0000000010 0.6931470755 0.0000001051 ;; 12 : 0.6931471806 0.0000000001 0.6931471726 0.0000000079 ;; 13 : 0.6931471806 0.0000000000 0.6931471747 0.0000000058 ;; 14 : 0.6931471806 0.0000000000 0.6931471801 0.0000000004 ;; 15 : 0.6931471806 0.0000000000 0.6931471802 0.0000000003 ;; 16 : 0.6931471806 0.0000000000 0.6931471805 0.0000000000 ;; 17 : 0.6931471806 0.0000000000 0.6931471805 0.0000000000 ;; 18 : 0.6931471806 0.0000000000 0.6931471806 0.0000000000 ;; 19 : 0.6931471806 0.0000000000 0.6931471806 0.0000000000 ;; 20 : 0.6931471806 0.0000000000 0.6931471806 0.0000000000 ;; ---------------------------------------------------------------------------------------------------- ;; The closed weights have the following formula: ;; $$\\frac{-1^{n-r}}{r!(n-r)!}\\int_0^n\\left(\\prod_{k=0}^{r-1}(t-k)\\right)\\left(\\prod_{k=r+1}^n(t-k)\\right) dx$$ ;; ---------------------------------------------------------------------------------------------------- ;; The closed weights were computed like so: ;; (loop for n from 1 upto 20 ;; collect (loop for r from 0 upto n ;; collect (* (if (evenp (- n r)) 1 -1) ;; (/ (* (mjr_combe_! r) ;; (mjr_combe_! (- n r)))) ;; (mjr_poly_intg (mjr_poly_make-from-roots (loop for i from 0 upto n ;; when (not (= i r)) ;; collect i)) ;; 0 n)))) ;; ---------------------------------------------------------------------------------------------------- ;; The open weights have the following formula: ;; $$\\frac{-1^{n-r}}{(r-1)!(1+n-r)!}\\int_0^n\\left(\\prod_{k=0}^{r-1}(t-k)\\right)\\left(\\prod_{k=r+1}^n(t-k)\\right) dx$$ ;; ---------------------------------------------------------------------------------------------------- ;; The closed weights were computed like so: ;; (loop for n from 3 upto 20 ;; collect (loop for r from 1 upto (1- n) ;; collect (* (if (evenp (1+ (- n r))) 1 -1) ;; (/ (* (mjr_combe_! (1- r)) ;; (mjr_combe_! (- n (1+ r))))) ;; (mjr_poly_intg (mjr_poly_make-from-roots (loop for i from 1 upto (1- n) ;; when (not (= i r)) ;; collect i)) ;; 0 n)))) ;; ---------------------------------------------------------------------------------------------------- (let ((fun (mjr_mxp_string-or-func-to-lambda fun "x"))) (if closed (let ((wts #(#(1/2 1/2) #(1/3 4/3 1/3) #(3/8 9/8 9/8 3/8) #(14/45 64/45 8/15 64/45 14/45) #(95/288 125/96 125/144 125/144 125/96 95/288) #(41/140 54/35 27/140 68/35 27/140 54/35 41/140) #(5257/17280 25039/17280 343/640 20923/17280 20923/17280 343/640 25039/17280 5257/17280) #(3956/14175 23552/14175 -3712/14175 41984/14175 -3632/2835 41984/14175 -3712/14175 23552/14175 3956/14175) #(25713/89600 141669/89600 243/2240 10881/5600 26001/44800 26001/44800 10881/5600 243/2240 141669/89600 25713/89600) #(80335/299376 132875/74844 -80875/99792 28375/6237 -24125/5544 89035/12474 -24125/5544 28375/6237 -80875/99792 132875/74844 80335/299376) #(4777223/17418240 49450643/29030400 -35608243/87091200 6166523/1935360 -17591827/14515200 28404871/14515200 28404871/14515200 -17591827/14515200 6166523/1935360 -35608243/87091200 49450643/29030400 4777223/17418240) #(1364651/5255250 150048/79625 -1264644/875875 3572512/525525 -3432753/350350 14586048/875875 -2090408/125125 14586048/875875 -3432753/350350 3572512/525525 -1264644/875875 150048/79625 1364651/5255250) #(106364763817/402361344000 731649485593/402361344000 -22582626859/22353408000 144926245243/28740096000 -78862978129/16094453760 298542743759/44706816000 -46704658663/33530112000 -46704658663/33530112000 298542743759/44706816000 -78862978129/16094453760 144926245243/28740096000 -22582626859/22353408000 731649485593/402361344000 106364763817/402361344000) #(631693279/2501928000 311056753/156370500 -5395044599/2501928000 765940609/78185250 -46375653541/2501928000 5525678207/156370500 -39205297537/833976000 712193069/13030875 -39205297537/833976000 5525678207/156370500 -46375653541/2501928000 765940609/78185250 -5395044599/2501928000 311056753/156370500 631693279/2501928000) #(25221445/98402304 442589775/229605376 -388226775/229605376 1746295975/229605376 -372104925/32800768 4103141115/229605376 -10001664025/688816128 1698012675/229605376 1698012675/229605376 -10001664025/688816128 4103141115/229605376 -372104925/32800768 1746295975/229605376 -388226775/229605376 442589775/229605376 25221445/98402304) #(120348894184/488462349375 1021012852736/488462349375 -31952201728/10854718875 1331538968576/97692469875 -280654342912/8881133625 3713412349952/54273594375 -54452275263488/488462349375 14990200029184/97692469875 -606473420576/3618239625 14990200029184/97692469875 -54452275263488/488462349375 3713412349952/54273594375 -280654342912/8881133625 1331538968576/97692469875 -31952201728/10854718875 1021012852736/488462349375 120348894184/488462349375) #(85455477715379/342372925440000 170070083613041/83691159552000 -2303599607287621/941525544960000 294905863842803/26900729856000 -37533141482453/1743565824000 2688251482636513/67251824640000 -9269566080827161/188305108992000 4755096326451617/104613949440000 -6391636155891919/376610217984000 -6391636155891919/376610217984000 4755096326451617/104613949440000 -9269566080827161/188305108992000 2688251482636513/67251824640000 -37533141482453/1743565824000 294905863842803/26900729856000 -2303599607287621/941525544960000 170070083613041/83691159552000 85455477715379/342372925440000) #(611197056507/2534852320000 55461906657/25348523200 -1927646624637/506970464000 1455921591759/79214135000 -1588823650719/31685654000 2777157949977/22632610000 -2128376034297/9053044000 5981920958997/15842827000 -628691926671297/1267426160000 34599550825551/63371308000 -628691926671297/1267426160000 5981920958997/15842827000 -2128376034297/9053044000 2777157949977/22632610000 -1588823650719/31685654000 1455921591759/79214135000 -1927646624637/506970464000 55461906657/25348523200 611197056507/2534852320000) #(1311546499957236437/5377993912811520000 18351023301032567/8604790260498432 -10034170825244579/3064383995904000 495333631700360617/32593902501888000 -5459170901406129527/149388719800320000 1309851809386170571/16598746644480000 -76111777777911661/609749876736000 13908243317345626267/89633231880192000 -2609922351726730283/19918495973376000 143909204406256715953/2688996956405760000 143909204406256715953/2688996956405760000 -2609922351726730283/19918495973376000 13908243317345626267/89633231880192000 -76111777777911661/609749876736000 1309851809386170571/16598746644480000 -5459170901406129527/149388719800320000 495333631700360617/32593902501888000 -10034170825244579/3064383995904000 18351023301032567/8604790260498432 1311546499957236437/5377993912811520000) #(1145302367137/4842604238472 3355823042500/1470076286679 -97339548544375/20581068013506 82748714972500/3430178002251 -2069649611963125/27441424018008 101305879622128/490025428893 -1557905611303750/3430178002251 2869553648930000/3430178002251 -2511881305088125/1960101715572 17040565224805000/10290534006753 -1684005984173647/935503091523 17040565224805000/10290534006753 -2511881305088125/1960101715572 2869553648930000/3430178002251 -1557905611303750/3430178002251 101305879622128/490025428893 -2069649611963125/27441424018008 82748714972500/3430178002251 -97339548544375/20581068013506 3355823042500/1470076286679 1145302367137/4842604238472))) (h (/ (- end start) order))) (* h (loop for n from 0 upto order for w across (aref wts (1- order)) for x = start then (+ x h) sum (* w (funcall fun x))))) (let ((wts #(#(2) #(3/2 3/2) #(8/3 -4/3 8/3) #(55/24 5/24 5/24 55/24) #(33/10 -21/5 39/5 -21/5 33/10) #(4277/1440 -1057/480 1967/720 1967/720 -1057/480 4277/1440) #(736/189 -848/105 1952/105 -19672/945 1952/105 -848/105 736/189) #(16083/4480 -25227/4480 44703/4480 -15399/4480 -15399/4480 44703/4480 -25227/4480 16083/4480) #(20225/4536 -4175/324 41675/1134 -137675/2268 169555/2268 -137675/2268 41675/1134 -4175/324 20225/4536) #(4325321/1036800 -72635189/7257600 4310317/181440 -11746361/453600 48901919/3628800 48901919/3628800 -11746361/453600 4310317/181440 -72635189/7257600 4325321/1036800) #(9626/1925 -35771/1925 123058/1925 -266298/1925 427956/1925 -494042/1925 427956/1925 -266298/1925 123058/1925 -35771/1925 9626/1925) #(4527766399/958003200 -4881098833/319334400 43831288241/958003200 -1651951561/21288960 13602071249/159667200 -1158778153/31933440 -1158778153/31933440 13602071249/159667200 -1651951561/21288960 43831288241/958003200 -4881098833/319334400 4527766399/958003200) #(2303435659/416988000 -1746642583/69498000 7067957029/69498000 -11311677097/41698800 1670445427/3088800 -9326809303/11583000 32009894483/34749000 -9326809303/11583000 1670445427/3088800 -11311677097/41698800 7067957029/69498000 -1746642583/69498000 2303435659/416988000) #(905730205/172204032 -3689759795/172204032 2226544555/28700672 -2147796905/12300288 47259022115/172204032 -15560694095/57401344 1684042445/14350336 1684042445/14350336 -15560694095/57401344 47259022115/172204032 -2147796905/12300288 2226544555/28700672 -3689759795/172204032 905730205/172204032) #(11555275136/1915538625 -62273397568/1915538625 290404217984/1915538625 -919494024608/1915538625 2192567376256/1915538625 -3992965568192/1915538625 1897702417792/638512875 -2131717882256/638512875 1897702417792/638512875 -3992965568192/1915538625 2192567376256/1915538625 -919494024608/1915538625 290404217984/1915538625 -62273397568/1915538625 11555275136/1915538625) #(362555126427073/62768369664000 -1782924415592401/62768369664000 7581522218007113/62768369664000 -7030469481873619/20922789888000 6061162410557827/8966909952000 -60117144155671973/62768369664000 56621644374387437/62768369664000 -7822950144565727/20922789888000 -7822950144565727/20922789888000 56621644374387437/62768369664000 -60117144155671973/62768369664000 6061162410557827/8966909952000 -7030469481873619/20922789888000 7581522218007113/62768369664000 -1782924415592401/62768369664000 362555126427073/62768369664000) #(62209540161/9529520000 -193893576633/4764760000 186372621/866320 -750890558259/952952000 208563085359/95295200 -22447758358269/4764760000 9590151541221/1191190000 -2105424409923/190590400 11691873278343/952952000 -2105424409923/190590400 9590151541221/1191190000 -22447758358269/4764760000 208563085359/95295200 -750890558259/952952000 186372621/866320 -193893576633/4764760000 62209540161/9529520000) #(57424625956493833/9146248151040000 -51499668595907713/1422749712384000 2831206484232759701/16005934264320000 -373766725994602289/640237370572800 337230604409417/237124952064 -266340535976610163/103934638080000 986233830545113243/291016986624000 -5408953361592551857/1778437140480000 1586799402762999283/1280474741145600 1586799402762999283/1280474741145600 -5408953361592551857/1778437140480000 986233830545113243/291016986624000 -266340535976610163/103934638080000 337230604409417/237124952064 -373766725994602289/640237370572800 2831206484232759701/16005934264320000 -51499668595907713/1422749712384000 57424625956493833/9146248151040000) #(4373703751565/623668727682 -20649923434495/415779151788 61065477359735/207889575894 -253893379558405/207889575894 402631871310710/103944787947 -142866589182035/14849255421 284667128261570/14849255421 -6459981265186655/207889575894 4301173329207145/103944787947 -28374376947229255/623668727682 4301173329207145/103944787947 -6459981265186655/207889575894 284667128261570/14849255421 -142866589182035/14849255421 402631871310710/103944787947 -253893379558405/207889575894 61065477359735/207889575894 -20649923434495/415779151788 4373703751565/623668727682))) (h (/ (- end start) order))) (* h (loop for n from 1 upto (1- order) for w across (aref wts (- order 2)) for x = (+ start h) then (+ x h) sum (* w (funcall fun x)))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-gauss-legendre (fun &key start end (order 10)) "Compute the definite integral using Gaussian quadrature (1 start end) (error "mjr_intg_simple-gauss-legendre: Lower limit (START) must be numerically less than upper limit (END)!")) ((not (integerp order)) (error "mjr_intg_simple-gauss-legendre: ORDER must be an integer!")) ((< order 2) (error "mjr_intg_simple-gauss-legendre: ORDER must be greater than 2!")) ((> order 20) (error "mjr_intg_simple-gauss-legendre: ORDER must be less than 21!"))) (let ((wid (/ (- end start) 2)) (mid (/ (+ end start) 2))) (* wid (loop for a in (aref abscissa order) for w in (aref weights order) sum (* w (funcall fun (+ (* wid a) mid))))))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_simple-gauss-kronrod (fun &key start end (order 31)) "Return two estimates (a high order one and a lower order one) for the integral of f on [start, end] using Gauss-Kronrod. The order argument determines the order of the Kronrod rule and the embedded Gauss rule: * order=15 => 15-point Gauss-Kronrod + 7-point Gauss * order=21 => 21-point Gauss-Kronrod + 10-point Gauss * order=31 => 31-point Gauss-Kronrod + 15-point Gauss * order=41 => 41-point Gauss-Kronrod & 20-point Gauss * order=51 => 51-point Gauss-Kronrod & 25-point Gauss * order=61 => 61-point Gauss-Kronrod & 30-point Gauss References: The quadrature weights and abscissas (computed to 80 decimal digits) are courtesy of L. W. Fullerton, Bell Labs, Nov. 1981." (cond ((not (member order '(15 21 31 41 51 61))) (error "mjr_intg_simple-gauss-kronrod: Order must be 15, 21, 31, 41, 51, or 61!")) ((not (numberp start)) (error "mjr_intg_simple-gauss-kronrod: start must be a number!")) ((not (numberp end)) (error "mjr_intg_simple-gauss-kronrod: end must be a number!")) ((< end start) (error "mjr_intg_simple-gauss-kronrod: start must be less than end!"))) (let* ((fun (mjr_mxp_string-or-func-to-lambda fun "x")) (oidx (1- (truncate order 10))) (wg (aref #(#(0.129484966168869693270611432679082d0 0.279705391489276667901467771423780d0 0.381830050505118944950369775488975d0 0.417959183673469387755102040816327d0) #(0.066671344308688137593568809893332d0 0.149451349150580593145776339657697d0 0.219086362515982043995534934228163d0 0.269266719309996355091226921569469d0 0.295524224714752870173892994651338d0) #(0.030753241996117268354628393577204d0 0.070366047488108124709267416450667d0 0.107159220467171935011869546685869d0 0.139570677926154314447804794511028d0 0.166269205816993933553200860481209d0 0.186161000015562211026800561866423d0 0.198431485327111576456118326443839d0 0.202578241925561272880620199967519d0) #(0.017614007139152118311861962351853d0 0.040601429800386941331039952274932d0 0.062672048334109063569506535187042d0 0.083276741576704748724758143222046d0 0.101930119817240435036750135480350d0 0.118194531961518417312377377711382d0 0.131688638449176626898494499748163d0 0.142096109318382051329298325067165d0 0.149172986472603746787828737001969d0 0.152753387130725850698084331955098d0) #(0.011393798501026287947902964113235d0 0.026354986615032137261901815295299d0 0.040939156701306312655623487711646d0 0.054904695975835191925936891540473d0 0.068038333812356917207187185656708d0 0.080140700335001018013234959669111d0 0.091028261982963649811497220702892d0 0.100535949067050644202206890392686d0 0.108519624474263653116093957050117d0 0.114858259145711648339325545869556d0 0.119455763535784772228178126512901d0 0.122242442990310041688959518945852d0 0.123176053726715451203902873079050d0) #(0.007968192496166605615465883474674d0 0.018466468311090959142302131912047d0 0.028784707883323369349719179611292d0 0.038799192569627049596801936446348d0 0.048402672830594052902938140422808d0 0.057493156217619066481721689402056d0 0.065974229882180495128128515115962d0 0.073755974737705206268243850022191d0 0.080755895229420215354694938460530d0 0.086899787201082979802387530715126d0 0.092122522237786128717632707087619d0 0.096368737174644259639468626351810d0 0.099593420586795267062780282103569d0 0.101762389748405504596428952168554d0 0.102852652893558840341285636705415d0)) oidx)) (xgk (aref #(#(0.991455371120812639206854697526329d0 0.949107912342758524526189684047851d0 0.864864423359769072789712788640926d0 0.741531185599394439863864773280788d0 0.586087235467691130294144838258730d0 0.405845151377397166906606412076961d0 0.207784955007898467600689403773245d0 0.000000000000000000000000000000000d0) #(0.995657163025808080735527280689003d0 0.973906528517171720077964012084452d0 0.930157491355708226001207180059508d0 0.865063366688984510732096688423493d0 0.780817726586416897063717578345042d0 0.679409568299024406234327365114874d0 0.562757134668604683339000099272694d0 0.433395394129247190799265943165784d0 0.294392862701460198131126603103866d0 0.148874338981631210884826001129720d0 0.000000000000000000000000000000000d0) #(0.998002298693397060285172840152271d0 0.987992518020485428489565718586613d0 0.967739075679139134257347978784337d0 0.937273392400705904307758947710209d0 0.897264532344081900882509656454496d0 0.848206583410427216200648320774217d0 0.790418501442465932967649294817947d0 0.724417731360170047416186054613938d0 0.650996741297416970533735895313275d0 0.570972172608538847537226737253911d0 0.485081863640239680693655740232351d0 0.394151347077563369897207370981045d0 0.299180007153168812166780024266389d0 0.201194093997434522300628303394596d0 0.101142066918717499027074231447392d0 0.000000000000000000000000000000000d0) #(0.998859031588277663838315576545863d0 0.993128599185094924786122388471320d0 0.981507877450250259193342994720217d0 0.963971927277913791267666131197277d0 0.940822633831754753519982722212443d0 0.912234428251325905867752441203298d0 0.878276811252281976077442995113078d0 0.839116971822218823394529061701521d0 0.795041428837551198350638833272788d0 0.746331906460150792614305070355642d0 0.693237656334751384805490711845932d0 0.636053680726515025452836696226286d0 0.575140446819710315342946036586425d0 0.510867001950827098004364050955251d0 0.443593175238725103199992213492640d0 0.373706088715419560672548177024927d0 0.301627868114913004320555356858592d0 0.227785851141645078080496195368575d0 0.152605465240922675505220241022678d0 0.076526521133497333754640409398838d0 0.000000000000000000000000000000000d0) #(0.999262104992609834193457486540341d0 0.995556969790498097908784946893902d0 0.988035794534077247637331014577406d0 0.976663921459517511498315386479594d0 0.961614986425842512418130033660167d0 0.942974571228974339414011169658471d0 0.920747115281701561746346084546331d0 0.894991997878275368851042006782805d0 0.865847065293275595448996969588340d0 0.833442628760834001421021108693570d0 0.797873797998500059410410904994307d0 0.759259263037357630577282865204361d0 0.717766406813084388186654079773298d0 0.673566368473468364485120633247622d0 0.626810099010317412788122681624518d0 0.577662930241222967723689841612654d0 0.526325284334719182599623778158010d0 0.473002731445714960522182115009192d0 0.417885382193037748851814394594572d0 0.361172305809387837735821730127641d0 0.303089538931107830167478909980339d0 0.243866883720988432045190362797452d0 0.183718939421048892015969888759528d0 0.122864692610710396387359818808037d0 0.061544483005685078886546392366797d0 0.000000000000000000000000000000000d0) #(0.999484410050490637571325895705811d0 0.996893484074649540271630050918695d0 0.991630996870404594858628366109486d0 0.983668123279747209970032581605663d0 0.973116322501126268374693868423707d0 0.960021864968307512216871025581798d0 0.944374444748559979415831324037439d0 0.926200047429274325879324277080474d0 0.905573307699907798546522558925958d0 0.882560535792052681543116462530226d0 0.857205233546061098958658510658944d0 0.829565762382768397442898119732502d0 0.799727835821839083013668942322683d0 0.767777432104826194917977340974503d0 0.733790062453226804726171131369528d0 0.697850494793315796932292388026640d0 0.660061064126626961370053668149271d0 0.620526182989242861140477556431189d0 0.579345235826361691756024932172540d0 0.536624148142019899264169793311073d0 0.492480467861778574993693061207709d0 0.447033769538089176780609900322854d0 0.400401254830394392535476211542661d0 0.352704725530878113471037207089374d0 0.304073202273625077372677107199257d0 0.254636926167889846439805129817805d0 0.204525116682309891438957671002025d0 0.153869913608583546963794672743256d0 0.102806937966737030147096751318001d0 0.051471842555317695833025213166723d0 0.000000000000000000000000000000000d0)) oidx)) (wgk (aref #(#(0.022935322010529224963732008058970d0 0.063092092629978553290700663189204d0 0.104790010322250183839876322541518d0 0.140653259715525918745189590510238d0 0.169004726639267902826583426598550d0 0.190350578064785409913256402421014d0 0.204432940075298892414161999234649d0 0.209482141084727828012999174891714d0) #(0.011694638867371874278064396062192d0 0.032558162307964727478818972459390d0 0.054755896574351996031381300244580d0 0.075039674810919952767043140916190d0 0.093125454583697605535065465083366d0 0.109387158802297641899210590325805d0 0.123491976262065851077958109831074d0 0.134709217311473325928054001771707d0 0.142775938577060080797094273138717d0 0.147739104901338491374841515972068d0 0.149445554002916905664936468389821d0) #(0.005377479872923348987792051430128d0 0.015007947329316122538374763075807d0 0.025460847326715320186874001019653d0 0.035346360791375846222037948478360d0 0.044589751324764876608227299373280d0 0.053481524690928087265343147239430d0 0.062009567800670640285139230960803d0 0.069854121318728258709520077099147d0 0.076849680757720378894432777482659d0 0.083080502823133021038289247286104d0 0.088564443056211770647275443693774d0 0.093126598170825321225486872747346d0 0.096642726983623678505179907627589d0 0.099173598721791959332393173484603d0 0.100769845523875595044946662617570d0 0.101330007014791549017374792767493d0) #(0.003073583718520531501218293246031d0 0.008600269855642942198661787950102d0 0.014626169256971252983787960308868d0 0.020388373461266523598010231432755d0 0.025882133604951158834505067096153d0 0.031287306777032798958543119323801d0 0.036600169758200798030557240707211d0 0.041668873327973686263788305936895d0 0.046434821867497674720231880926108d0 0.050944573923728691932707670050345d0 0.055195105348285994744832372419777d0 0.059111400880639572374967220648594d0 0.062653237554781168025870122174255d0 0.065834597133618422111563556969398d0 0.068648672928521619345623411885368d0 0.071054423553444068305790361723210d0 0.073030690332786667495189417658913d0 0.074582875400499188986581418362488d0 0.075704497684556674659542775376617d0 0.076377867672080736705502835038061d0 0.076600711917999656445049901530102d0) #(0.001987383892330315926507851882843d0 0.005561932135356713758040236901066d0 0.009473973386174151607207710523655d0 0.013236229195571674813656405846976d0 0.016847817709128298231516667536336d0 0.020435371145882835456568292235939d0 0.024009945606953216220092489164881d0 0.027475317587851737802948455517811d0 0.030792300167387488891109020215229d0 0.034002130274329337836748795229551d0 0.037116271483415543560330625367620d0 0.040083825504032382074839284467076d0 0.042872845020170049476895792439495d0 0.045502913049921788909870584752660d0 0.047982537138836713906392255756915d0 0.050277679080715671963325259433440d0 0.052362885806407475864366712137873d0 0.054251129888545490144543370459876d0 0.055950811220412317308240686382747d0 0.057437116361567832853582693939506d0 0.058689680022394207961974175856788d0 0.059720340324174059979099291932562d0 0.060539455376045862945360267517565d0 0.061128509717053048305859030416293d0 0.061471189871425316661544131965264d0 0.061580818067832935078759824240055d0) #(0.001389013698677007624551591226760d0 0.003890461127099884051267201844516d0 0.006630703915931292173319826369750d0 0.009273279659517763428441146892024d0 0.011823015253496341742232898853251d0 0.014369729507045804812451432443580d0 0.016920889189053272627572289420322d0 0.019414141193942381173408951050128d0 0.021828035821609192297167485738339d0 0.024191162078080601365686370725232d0 0.026509954882333101610601709335075d0 0.028754048765041292843978785354334d0 0.030907257562387762472884252943092d0 0.032981447057483726031814191016854d0 0.034979338028060024137499670731468d0 0.036882364651821229223911065617136d0 0.038678945624727592950348651532281d0 0.040374538951535959111995279752468d0 0.041969810215164246147147541285970d0 0.043452539701356069316831728117073d0 0.044814800133162663192355551616723d0 0.046059238271006988116271735559374d0 0.047185546569299153945261478181099d0 0.048185861757087129140779492298305d0 0.049055434555029778887528165367238d0 0.049795683427074206357811569379942d0 0.050405921402782346840893085653585d0 0.050881795898749606492297473049805d0 0.051221547849258772170656282604944d0 0.051426128537459025933862879215781d0 0.051494729429451567558340433647099d0)) oidx)) (centr (/ (+ start end) 2)) (hlgth (/ (- end start) 2)) (fc (funcall fun centr)) (resg (if (oddp oidx) 0.0d0 (* fc (aref wg (1- (length wg)))))) (resk (* fc (aref wgk (1- (length wgk)))))) (loop for j from 1 upto (truncate order 2) for absc = (* hlgth (aref xgk (1- j))) for fsum = (+ (funcall fun (- centr absc)) (funcall fun (+ centr absc))) do (incf resk (* (aref wgk (1- j)) fsum)) when (evenp j) do (incf resg (* (aref wg (1- (/ j 2))) fsum))) (values (* resk hlgth) ;; Gauss-Kronrod integral value (* resg hlgth)))) ;; Gauss integral value ;; (loop for ord in '(15 21 31 41 51 61) ;; for (k g) = (multiple-value-list (gk (lambda (x) (sin (* 20 x))) :start 1 :end 2 :order ord)) ;; for e = (abs (- k g)) ;; do (format 't "~10d :: ~30,20f ~30,20f ~30,20f~%" ord k g e)) ;; 15 :: 0.05375100448700819000 0.03329022396630121000 0.02046078052070697500 ;; 21 :: 0.05375100604962934000 0.05376956145745758600 0.00001855540782824683 ;; 31 :: 0.05375100608944126400 0.05375100600323407000 0.00000000008620719416 ;; 41 :: 0.05375100610990241000 0.05375100617328299000 0.00000000006338057856 ;; 51 :: 0.05375100612233012000 0.05375100607136562600 0.00000000005096449651 ;; 61 :: 0.05375100613067555000 0.05375100617328321000 0.00000000004260766046 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun mjr_intg_gbl-adp-factory-order (fun &key start end order-min order-max order-step (err 0.0002) (imethod #'mjr_intg_simple-gauss-legendre) show-progress) "Compute the definite integral of FUN between A and B using the given integration method. The order is stepped from order-min to order-max by order-step until the difference in two successive iterations is less than err. Here are some possible values for method: * #'mjr_intg_simple-newton-cotes -- :order-min 3 :order-max 15 :order-step 1 * #'mjr_intg_simple-gauss-legendre -- :order-min 2 :order-max 20 :order-step 1 * #'mjr_intg_simple-monte-carlo -- :order-min 1000 :order-max 100000 :order-step 1000 * #'mjr_intg_simple-gauss-kronrod -- :order-min 21 :order-max 61 :order-step 10 The values for order-min, order-max, and order-step will be set as above if all three are NIL." (if (not (or order-min order-max order-step)) (cond ((equalp imethod #'mjr_intg_simple-gauss-legendre) (mjr_intg_gbl-adp-factory-order fun :start start :end end :err err :imethod imethod :show-progress show-progress :order-min 2 :order-max 20 :order-step 1)) ((equalp imethod #'mjr_intg_simple-newton-cotes) (mjr_intg_gbl-adp-factory-order fun :start start :end end :err err :imethod imethod :show-progress show-progress :order-min 2 :order-max 15 :order-step 1)) ((equalp imethod #'mjr_intg_simple-monte-carlo) (mjr_intg_gbl-adp-factory-order fun :start start :end end :err err :imethod imethod :show-progress show-progress :order-min 1000 :order-max 100000 :order-step 1000)) ((equalp imethod #'mjr_intg_simple-gauss-kronrod) (mjr_intg_gbl-adp-factory-order fun :start start :end end :err err :imethod imethod :show-progress show-progress :order-min 21 :order-max 61 :order-step 10))) (loop for cur-ord from order-min upto order-max by order-step for prv-val = nil then cur-val for cur-val = (funcall imethod fun :start start :end end :order cur-ord) do (if show-progress (format 't "~10d : ~25,10f~%" cur-ord cur-val)) when (and prv-val (mjr_eps_= prv-val cur-val err)) do (return cur-val))))