Numerical Quadrature

ToC

TOC

Back to index/Generated with Clerk from src/emmy/numerical/quadrature.cljc@f807468

Numerical Quadrature

This namespace unites all of the work inside emmy.numerical.quadrature behind a single interface, fronted by the all-powerful definite-integral function.

The interface takes f, an integrand, along with bounds a and b:

(definite-integral f a b)

Optionally, you can provide a dictionary of customizing options. These are passed down to whatever method you supply via the :method key.

(definite-integral f a b opts)

Implementation

The keys in quad-methods below define the full range of integration methods available in the package. Each entry in this dictionary is either:

An 'integrator' function that matches the interface above for definite-integral (possibly created with qc/defintegrator)
a dictionary of extra options. This must contain a :method key.

This latter style is used when the method itself is a specialization of a more general method.

(def ^:private quadrature-methods

  {:open                    {:method :adaptive-bulirsch-stoer

                             :interval qc/open}

   :closed                  {:method :adaptive-bulirsch-stoer

                             :interval qc/closed}

   :closed-open             {:method :adaptive-bulirsch-stoer

                             :interval qc/closed-open}

   :open-closed             {:method :adaptive-bulirsch-stoer

                             :interval qc/open-closed}

   :bulirsch-stoer-open     bs/open-integral

   :bulirsch-stoer-closed   bs/closed-integral

   :adaptive-bulirsch-stoer (qa/adaptive bs/open-integral bs/closed-integral)

   :left-riemann            riemann/left-integral

   :right-riemann           riemann/right-integral

   :lower-riemann           riemann/lower-integral

   :upper-riemann           riemann/upper-integral

   :midpoint                mid/integral

   :trapezoid               trap/integral

   :boole                   boole/integral

   :milne                   milne/integral

   :simpson                 simp/integral

   :simpson38               simp38/integral

   :romberg                 romberg/closed-integral

   :romberg-open            romberg/open-integral})

{:adaptive-bulirsch-stoer #object[emmy.numerical.quadrature.adaptive$adaptive$rec__53942 0x6f73f5e4 "

emmy.numerical.quadrature.adaptive$adaptive$rec__53942@6f73f5e4"

] :boole #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x79d68875 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@79d68875"

] :bulirsch-stoer-closed #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x63996e8e "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@63996e8e"

] :bulirsch-stoer-open #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x35c5d807 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@35c5d807"

] :closed {:interval [:emmy.numerical.quadrature.common/closed :emmy.numerical.quadrature.common/closed] :method :adaptive-bulirsch-stoer} :closed-open {:interval [:emmy.numerical.quadrature.common/closed :emmy.numerical.quadrature.common/open] :method :adaptive-bulirsch-stoer} :left-riemann #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x45e9d78e "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@45e9d78e"

] :lower-riemann #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x8b17a66 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@8b17a66"

] :midpoint #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x207488e0 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@207488e0"

] :milne #object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x3b6a6909 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@3b6a6909"

] 9 more elided}

(def available-methods

  (into #{} (keys quadrature-methods)))

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

The user can specify a method by providing the :method key in their options with:

a key in the above dictionary
another dict
a custom integration function

The latter two are the allowed value types in quadrature-methods.

(defn- extract-method

  "Attempts to turn the supplied argument into an integration method; returns nil

  if method doesn't exist."

  [method]

  (cond (fn? method)

        [method {}]

        (keyword? method)

        (extract-method

         (quadrature-methods method))

        (map? method)

        (let [[f m] (extract-method

                     (:method method))]

          [f (merge (dissoc method :method) m)])))

#object[emmy.numerical.quadrature$extract_method 0x30b49c4d "

emmy.numerical.quadrature$extract_method@30b49c4d"

]

(defn get-integrator

  "Takes:

  - An integration method, specified as either:

    - 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

  - `a` and `b` integration endpoints

  - an optional dictionary of options `m`

  And returns a pair of an integrator function and a possibly-enhanced options

  dictionary.

  (Some integration functions require extra options, so the returned dictionary

  may have more entries than the `m` you pass in.)

  If either endpoint is infinite, the returned integrator is wrapped in

  `qi/improper` and able to handle infinite endpoints (as well as non-infinite

  endpoints by passing through directly to the underlying integrator)."

  ([method a b] (get-integrator method a b {}))

  ([method a b m]

   (when-let [[integrate opts] (extract-method method)]

     (let [integrate (if (or (g/infinite? a)

                             (g/infinite? b))

                       (qi/improper integrate)

                       integrate)]

       [integrate (dissoc (merge opts m) :method)]))))

#object[emmy.numerical.quadrature$get_integrator 0x3355a277 "

emmy.numerical.quadrature$get_integrator@3355a277"

]

Final API

Here we are! The one function you need care about if you're interested in definite integrals. Learn to use this, and then dig in to the details of individual methods if you run into trouble or want to learn more. Enjoy!

(defn definite-integral

  "Evaluates the definite integral of integrand `f` across the interval $a, b$.

  Optionally accepts a dictionary `opts` of customizing options; All `opts` will

  be passed through to the supplied `integrate` functions.

  If you'd like more control, or to retrieve the integration function directly

  without looking it up via `:method` each time, see `get-integrator`.

  All supplied options are passed through to the underlying integrator; see the

  specific integrator for information on what options are available.

  ## Keyword arguments:

  `: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

  Defaults to `:open`, which specifies an adaptive bulirsch-stoer quadrature method.

  `:compile?` If true, the generic function will be simplified and compiled

  before execution.

  `: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."

  ([f a b] (definite-integral f a b {}))

  ([f a b {:keys [method compile? info?]

           :or {method :open

                compile? false

                info? false}

           :as opts}]

   (if-let [[integrate m] (get-integrator method a b opts)]

     (let [f      (if compile?

                    (c/compile-fn f 1)

                    #?(:cljs (comp u/double f)

                       :clj f))

           result (integrate f a b m)]

       (if info? result (:result result)))

     (u/illegal (str "Unknown method: " method

                     ". Try one of: "

                     available-methods)))))

#object[emmy.numerical.quadrature$definite_integral 0x37e21f75 "

emmy.numerical.quadrature$definite_integral@37e21f75"

]