Quadrature — sicmutils 0.14.0

Quadrature

Quadrature is the process of computing definite integrals of functions. A sugared default procedure for quadrature is provided, and we hope that it is adequate for most purposes.

(definite-integral <integrand>
                   <lower-limit> <upper-limit>
                   {:compile? true})

The integrand must be a real-valued function of a real argument. The limits of integration are specified as additional arguments.

Optionally you can supply a map of keyword arguments. The top level definite-integral function uses the following three arguments:

:compile? can be used to suppress compilation of the integrand, thus forcing it to be interpreted. This is usually to be ignored.
:info?: If true, definite-integral will return a map of integration information returned by the underlying integrator. Else, returns an estimate of the definite integral.
:method: Specifies the integration method used. Must be:
- a keyword naming one of the available methods in available-methods
- a function with the proper integrator signature
- a dictionary of integrator options with a :method key

:method defaults to :open, which specifies an adaptive bulirsch-stoer quadrature method. The other allowed / supported methods are:

:open
:closed
:closed-open
:open-closed
:bulirsch-stoer-open
:bulirsch-stoer-closed
:adaptive-bulirsch-stoer
:left-riemann
:right-riemann
:lower-riemann
:upper-riemann
:midpoint
:trapezoid
:boole
:milne
:simpson
:simpson38
:romberg
:romberg-open

The quadrature methods are all based on extrapolation. The Romberg method is a Richardson extrapolation of the trapezoid rule. It is usually worse than the other methods, which are adaptive rational function extrapolations of Trapezoid and Midpoint rules.

Closed integrators are best if we can include the endpoints of integration. This cannot be done if the endpoint is singular: thus the open formulas. Also, open formulas are forced when we have infinite limits.

Let’s do an example, it is as easy as pi!

(defn witch [x]
  (/ 4.0 (+ 1.0 (* x x))))

(definite-integral witch 0.0 1.0
  {:method :romberg :tolerance 1e-12})
;; => 3.141592653589793

Here’s another example for fun:

(defn foo [n]
  (let [f (fn [x] (expt (log (/ 1 x)) n))]
    (definite-integral
      f 0.0 1.0
      {:tolerance 1e-12
       :method :open-closed})))

(foo 0)
;;=> 1.0

(foo 1)
;;=> 0.9999999999983304

(foo 2)
;;=> 1.999999999998337

(foo 3)
;;=> 5.999999999998272

(foo 4)
;;=> 23.99999999949962

(foo 5)
;;=> 119.99999998778476

Do you recognize this function? What is (foo 6)?