<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Numerical integration (quadrature) — mpmath v0.17 documentation</title>
<link rel="stylesheet" href="../_static/default.css" type="text/css" />
<link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: '../',
VERSION: '0.17',
COLLAPSE_MODINDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="../_static/jquery.js"></script>
<script type="text/javascript" src="../_static/doctools.js"></script>
<link rel="top" title="mpmath v0.17 documentation" href="../index.html" />
<link rel="up" title="Numerical calculus" href="index.html" />
<link rel="next" title="Ordinary differential equations" href="odes.html" />
<link rel="prev" title="Differentiation" href="differentiation.html" />
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="../modindex.html" title="Global Module Index"
accesskey="M">modules</a> |</li>
<li class="right" >
<a href="odes.html" title="Ordinary differential equations"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="differentiation.html" title="Differentiation"
accesskey="P">previous</a> |</li>
<li><a href="../index.html">mpmath v0.17 documentation</a> »</li>
<li><a href="index.html" accesskey="U">Numerical calculus</a> »</li>
</ul>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="numerical-integration-quadrature">
<h1>Numerical integration (quadrature)<a class="headerlink" href="#numerical-integration-quadrature" title="Permalink to this headline">¶</a></h1>
<div class="section" id="standard-quadrature-quad">
<h2>Standard quadrature (<tt class="docutils literal"><span class="pre">quad</span></tt>)<a class="headerlink" href="#standard-quadrature-quad" title="Permalink to this headline">¶</a></h2>
<dl class="function">
<dt id="mpmath.quad">
<tt class="descclassname">mpmath.</tt><tt class="descname">quad</tt><big>(</big><em>ctx</em>, <em>f</em>, <em>*points</em>, <em>**kwargs</em><big>)</big><a class="headerlink" href="#mpmath.quad" title="Permalink to this definition">¶</a></dt>
<dd><p>Computes a single, double or triple integral over a given
1D interval, 2D rectangle, or 3D cuboid. A basic example:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">mpmath</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span><span class="p">;</span> <span class="n">mp</span><span class="o">.</span><span class="n">pretty</span> <span class="o">=</span> <span class="bp">True</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">sin</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">pi</span><span class="p">])</span>
<span class="go">2.0</span>
</pre></div>
</div>
<p>A basic 2D integral:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="o">+</span><span class="n">y</span><span class="o">/</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="n">pi</span><span class="o">/</span><span class="mi">2</span><span class="p">,</span> <span class="n">pi</span><span class="o">/</span><span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">pi</span><span class="p">])</span>
<span class="go">4.0</span>
</pre></div>
</div>
<p><strong>Interval format</strong></p>
<p>The integration range for each dimension may be specified
using a list or tuple. Arguments are interpreted as follows:</p>
<p><tt class="docutils literal"><span class="pre">quad(f,</span> <span class="pre">[x1,</span> <span class="pre">x2])</span></tt> – calculates
<img class="math" src="../_images/math/095e636757db8673ec95b39676d4474011cbf411.png" alt="\int_{x_1}^{x_2} f(x) \, dx"/></p>
<p><tt class="docutils literal"><span class="pre">quad(f,</span> <span class="pre">[x1,</span> <span class="pre">x2],</span> <span class="pre">[y1,</span> <span class="pre">y2])</span></tt> – calculates
<img class="math" src="../_images/math/fbf69f02d553b361aea2535bdfe760c00c47430d.png" alt="\int_{x_1}^{x_2} \int_{y_1}^{y_2} f(x,y) \, dy \, dx"/></p>
<p><tt class="docutils literal"><span class="pre">quad(f,</span> <span class="pre">[x1,</span> <span class="pre">x2],</span> <span class="pre">[y1,</span> <span class="pre">y2],</span> <span class="pre">[z1,</span> <span class="pre">z2])</span></tt> – calculates
<img class="math" src="../_images/math/ad6e7b009c965ade8dcc10e366197f421eede7cc.png" alt="\int_{x_1}^{x_2} \int_{y_1}^{y_2} \int_{z_1}^{z_2} f(x,y,z)
\, dz \, dy \, dx"/></p>
<p>Endpoints may be finite or infinite. An interval descriptor
may also contain more than two points. In this
case, the integration is split into subintervals, between
each pair of consecutive points. This is useful for
dealing with mid-interval discontinuities, or integrating
over large intervals where the function is irregular or
oscillates.</p>
<p><strong>Options</strong></p>
<p><a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> recognizes the following keyword arguments:</p>
<dl class="docutils">
<dt><em>method</em></dt>
<dd>Chooses integration algorithm (described below).</dd>
<dt><em>error</em></dt>
<dd>If set to true, <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> returns <img class="math" src="../_images/math/cc8877b37ecb041cbfa875ec771ce90175874bc4.png" alt="(v, e)"/> where <img class="math" src="../_images/math/a9f23bf124b6b2b2a993eb313c72e678664ac74a.png" alt="v"/> is the
integral and <img class="math" src="../_images/math/a3a59bb1293ee3f6dec19de4019a7178874219ae.png" alt="e"/> is the estimated error.</dd>
<dt><em>maxdegree</em></dt>
<dd>Maximum degree of the quadrature rule to try before
quitting.</dd>
<dt><em>verbose</em></dt>
<dd>Print details about progress.</dd>
</dl>
<p><strong>Algorithms</strong></p>
<p>Mpmath presently implements two integration algorithms: tanh-sinh
quadrature and Gauss-Legendre quadrature. These can be selected
using <em>method=’tanh-sinh’</em> or <em>method=’gauss-legendre’</em> or by
passing the classes <em>method=TanhSinh</em>, <em>method=GaussLegendre</em>.
The functions <tt class="xref docutils literal"><span class="pre">quadts()</span></tt> and <tt class="xref docutils literal"><span class="pre">quadgl()</span></tt> are also available
as shortcuts.</p>
<p>Both algorithms have the property that doubling the number of
evaluation points roughly doubles the accuracy, so both are ideal
for high precision quadrature (hundreds or thousands of digits).</p>
<p>At high precision, computing the nodes and weights for the
integration can be expensive (more expensive than computing the
function values). To make repeated integrations fast, nodes
are automatically cached.</p>
<p>The advantages of the tanh-sinh algorithm are that it tends to
handle endpoint singularities well, and that the nodes are cheap
to compute on the first run. For these reasons, it is used by
<a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> as the default algorithm.</p>
<p>Gauss-Legendre quadrature often requires fewer function
evaluations, and is therefore often faster for repeated use, but
the algorithm does not handle endpoint singularities as well and
the nodes are more expensive to compute. Gauss-Legendre quadrature
can be a better choice if the integrand is smooth and repeated
integrations are required (e.g. for multiple integrals).</p>
<p>See the documentation for <tt class="xref docutils literal"><span class="pre">TanhSinh</span></tt> and
<tt class="xref docutils literal"><span class="pre">GaussLegendre</span></tt> for additional details.</p>
<p><strong>Examples of 1D integrals</strong></p>
<p>Intervals may be infinite or half-infinite. The following two
examples evaluate the limits of the inverse tangent function
(<img class="math" src="../_images/math/61c9dca77f580d0d6ff196dbf084dbd46eb55c4a.png" alt="\int 1/(1+x^2) = \tan^{-1} x"/>), and the Gaussian integral
<img class="math" src="../_images/math/15a441c328cdcb543e4f1195ce6424f269cc9607.png" alt="\int_{\infty}^{\infty} \exp(-x^2)\,dx = \sqrt{\pi}"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">2</span><span class="o">/</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="mi">1</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">])</span>
<span class="go">3.14159265358979</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="n">inf</span><span class="p">,</span> <span class="n">inf</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span>
<span class="go">3.14159265358979</span>
</pre></div>
</div>
<p>Integrals can typically be resolved to high precision.
The following computes 50 digits of <img class="math" src="../_images/math/f2ca003a7da0de4994b4733e203b74ff52d42553.png" alt="\pi"/> by integrating the
area of the half-circle defined by <img class="math" src="../_images/math/87adb30179d12676da06d9c6a7a88648265fd94f.png" alt="x^2 + y^2 \le 1"/>,
<img class="math" src="../_images/math/f260c78fb4c935f83297e9c2cb2cf7c3c3e3adfc.png" alt="-1 \le x \le 1"/>, <img class="math" src="../_images/math/4a4e6e39a7375117fa73ed88a17f8ede37db56b0.png" alt="y \ge 0"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">50</span>
<span class="gp">>>> </span><span class="mi">2</span><span class="o">*</span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">3.1415926535897932384626433832795028841971693993751</span>
</pre></div>
</div>
<p>One can just as well compute 1000 digits (output truncated):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">1000</span>
<span class="gp">>>> </span><span class="mi">2</span><span class="o">*</span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span> <span class="c">#doctest:+ELLIPSIS</span>
<span class="go">3.141592653589793238462643383279502884...216420198</span>
</pre></div>
</div>
<p>Complex integrals are supported. The following computes
a residue at <img class="math" src="../_images/math/aa9f0d97d2f39f78f05e05da40bf04f5a7c0726c.png" alt="z = 0"/> by integrating counterclockwise along the
diamond-shaped path from <img class="math" src="../_images/math/dce34f4dfb2406144304ad0d6106c5382ddd1446.png" alt="1"/> to <img class="math" src="../_images/math/f07a8bc43a3f4fe9db689ab597112a5cf846e3a8.png" alt="+i"/> to <img class="math" src="../_images/math/bae5aba07d37ff6ff813107e76260fb31ad5794e.png" alt="-1"/> to <img class="math" src="../_images/math/0fda51aabf3fa68a72ba3b46518abc175c76dab2.png" alt="-i"/> to <img class="math" src="../_images/math/dce34f4dfb2406144304ad0d6106c5382ddd1446.png" alt="1"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">chop</span><span class="p">(</span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">z</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="n">z</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="n">j</span><span class="p">,</span><span class="mi">1</span><span class="p">]))</span>
<span class="go">(0.0 + 6.28318530717959j)</span>
</pre></div>
</div>
<p><strong>Examples of 2D and 3D integrals</strong></p>
<p>Here are several nice examples of analytically solvable
2D integrals (taken from MathWorld [1]) that can be evaluated
to high precision fairly rapidly by <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">30</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">:</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="p">((</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">*</span><span class="n">y</span><span class="p">)</span><span class="o">*</span><span class="n">log</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">y</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">0.577215664901532860606512090082</span>
<span class="gp">>>> </span><span class="o">+</span><span class="n">euler</span>
<span class="go">0.577215664901532860606512090082</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">y</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">3.17343648530607134219175646705</span>
<span class="gp">>>> </span><span class="mi">4</span><span class="o">*</span><span class="n">log</span><span class="p">(</span><span class="mi">2</span><span class="o">+</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">3</span><span class="p">))</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="mi">3</span>
<span class="go">3.17343648530607134219175646705</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">y</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">1.23370055013616982735431137498</span>
<span class="gp">>>> </span><span class="n">pi</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mi">8</span>
<span class="go">1.23370055013616982735431137498</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">*</span><span class="n">y</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">1.64493406684822643647241516665</span>
<span class="gp">>>> </span><span class="n">pi</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mi">6</span>
<span class="go">1.64493406684822643647241516665</span>
</pre></div>
</div>
<p>Multiple integrals may be done over infinite ranges:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="k">print</span><span class="p">(</span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">:</span> <span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">-</span><span class="n">y</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="n">inf</span><span class="p">]))</span>
<span class="go">0.367879441171442</span>
<span class="gp">>>> </span><span class="k">print</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="n">e</span><span class="p">)</span>
<span class="go">0.367879441171442</span>
</pre></div>
</div>
<p>For nonrectangular areas, one can call <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> recursively.
For example, we can replicate the earlier example of calculating
<img class="math" src="../_images/math/f2ca003a7da0de4994b4733e203b74ff52d42553.png" alt="\pi"/> by integrating over the unit-circle, and actually use double
quadrature to actually measure the area circle:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">y</span><span class="p">:</span> <span class="mi">1</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)])</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>
<span class="go">3.14159265358979</span>
</pre></div>
</div>
<p>Here is a simple triple integral:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">z</span><span class="p">:</span> <span class="n">x</span><span class="o">*</span><span class="n">y</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">z</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">method</span><span class="o">=</span><span class="s">'gauss-legendre'</span><span class="p">)</span>
<span class="go">0.101366277027041</span>
<span class="gp">>>> </span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">-</span><span class="n">log</span><span class="p">(</span><span class="mi">2</span><span class="p">))</span><span class="o">/</span><span class="mi">4</span>
<span class="go">0.101366277027041</span>
</pre></div>
</div>
<p><strong>Singularities</strong></p>
<p>Both tanh-sinh and Gauss-Legendre quadrature are designed to
integrate smooth (infinitely differentiable) functions. Neither
algorithm copes well with mid-interval singularities (such as
mid-interval discontinuities in <img class="math" src="../_images/math/c96dd6ec1dc4ad7520fbdc78fcdbec9edd068d0c.png" alt="f(x)"/> or <img class="math" src="../_images/math/bbad379658bd32e9b6479bbd666b93a96e6b48ba.png" alt="f'(x)"/>).
The best solution is to split the integral into parts:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="nb">abs</span><span class="p">(</span><span class="n">sin</span><span class="p">(</span><span class="n">x</span><span class="p">)),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">])</span> <span class="c"># Bad</span>
<span class="go">3.99900894176779</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="nb">abs</span><span class="p">(</span><span class="n">sin</span><span class="p">(</span><span class="n">x</span><span class="p">)),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">pi</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">])</span> <span class="c"># Good</span>
<span class="go">4.0</span>
</pre></div>
</div>
<p>The tanh-sinh rule often works well for integrands having a
singularity at one or both endpoints:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">log</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">method</span><span class="o">=</span><span class="s">'tanh-sinh'</span><span class="p">)</span> <span class="c"># Good</span>
<span class="go">-1.0</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">log</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">method</span><span class="o">=</span><span class="s">'gauss-legendre'</span><span class="p">)</span> <span class="c"># Bad</span>
<span class="go">-0.999932197413801</span>
</pre></div>
</div>
<p>However, the result may still be inaccurate for some functions:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">method</span><span class="o">=</span><span class="s">'tanh-sinh'</span><span class="p">)</span>
<span class="go">1.99999999946942</span>
</pre></div>
</div>
<p>This problem is not due to the quadrature rule per se, but to
numerical amplification of errors in the nodes. The problem can be
circumvented by temporarily increasing the precision:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">30</span>
<span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">method</span><span class="o">=</span><span class="s">'tanh-sinh'</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="o">+</span><span class="n">a</span>
<span class="go">2.0</span>
</pre></div>
</div>
<p><strong>Highly variable functions</strong></p>
<p>For functions that are smooth (in the sense of being infinitely
differentiable) but contain sharp mid-interval peaks or many
“bumps”, <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> may fail to provide full accuracy. For
example, with default settings, <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> is able to integrate
<img class="math" src="../_images/math/081bf0417da2140657aad909a14ca6c3ff8b77d5.png" alt="\sin(x)"/> accurately over an interval of length 100 but not over
length 1000:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">sin</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">100</span><span class="p">]);</span> <span class="mi">1</span><span class="o">-</span><span class="n">cos</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span> <span class="c"># Good</span>
<span class="go">0.137681127712316</span>
<span class="go">0.137681127712316</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">sin</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1000</span><span class="p">]);</span> <span class="mi">1</span><span class="o">-</span><span class="n">cos</span><span class="p">(</span><span class="mi">1000</span><span class="p">)</span> <span class="c"># Bad</span>
<span class="go">-37.8587612408485</span>
<span class="go">0.437620923709297</span>
</pre></div>
</div>
<p>One solution is to break the integration into 10 intervals of
length 100:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">sin</span><span class="p">,</span> <span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1000</span><span class="p">,</span> <span class="mi">10</span><span class="p">))</span> <span class="c"># Good</span>
<span class="go">0.437620923709297</span>
</pre></div>
</div>
<p>Another is to increase the degree of the quadrature:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">sin</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1000</span><span class="p">],</span> <span class="n">maxdegree</span><span class="o">=</span><span class="mi">10</span><span class="p">)</span> <span class="c"># Also good</span>
<span class="go">0.437620923709297</span>
</pre></div>
</div>
<p>Whether splitting the interval or increasing the degree is
more efficient differs from case to case. Another example is the
function <img class="math" src="../_images/math/9e304d281d5bb66802cb88a03ebfb328b1b276e8.png" alt="1/(1+x^2)"/>, which has a sharp peak centered around
<img class="math" src="../_images/math/2d348bde3e15456e71734dc2c56fc7425c95927f.png" alt="x = 0"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="mi">100</span><span class="p">,</span> <span class="mi">100</span><span class="p">])</span> <span class="c"># Bad</span>
<span class="go">3.64804647105268</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="mi">100</span><span class="p">,</span> <span class="mi">100</span><span class="p">],</span> <span class="n">maxdegree</span><span class="o">=</span><span class="mi">10</span><span class="p">)</span> <span class="c"># Good</span>
<span class="go">3.12159332021646</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="mi">100</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">100</span><span class="p">])</span> <span class="c"># Also good</span>
<span class="go">3.12159332021646</span>
</pre></div>
</div>
<p><strong>References</strong></p>
<ol class="arabic simple">
<li><a class="reference external" href="http://mathworld.wolfram.com/DoubleIntegral.html">http://mathworld.wolfram.com/DoubleIntegral.html</a></li>
</ol>
</dd></dl>
</div>
<div class="section" id="oscillatory-quadrature-quadosc">
<h2>Oscillatory quadrature (<tt class="docutils literal"><span class="pre">quadosc</span></tt>)<a class="headerlink" href="#oscillatory-quadrature-quadosc" title="Permalink to this headline">¶</a></h2>
<dl class="function">
<dt id="mpmath.quadosc">
<tt class="descclassname">mpmath.</tt><tt class="descname">quadosc</tt><big>(</big><em>ctx</em>, <em>f</em>, <em>interval</em>, <em>omega=None</em>, <em>period=None</em>, <em>zeros=None</em><big>)</big><a class="headerlink" href="#mpmath.quadosc" title="Permalink to this definition">¶</a></dt>
<dd><p>Calculates</p>
<div class="math">
<p><img src="../_images/math/faa759cbb4bff58ffe7cb487a95808076380f36d.png" alt="I = \int_a^b f(x) dx" /></p>
</div><p>where at least one of <img class="math" src="../_images/math/c7d457e388298246adb06c587bccd419ea67f7e8.png" alt="a"/> and <img class="math" src="../_images/math/8136a7ef6a03334a7246df9097e5bcc31ba33fd2.png" alt="b"/> is infinite and where
<img class="math" src="../_images/math/7ec57f0e091bf11c1a0afb6a3b36eb314f44af3f.png" alt="f(x) = g(x) \cos(\omega x + \phi)"/> for some slowly
decreasing function <img class="math" src="../_images/math/659a6a2610357101d77dbae8eb9bf0b3100c08fd.png" alt="g(x)"/>. With proper input, <a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a>
can also handle oscillatory integrals where the oscillation
rate is different from a pure sine or cosine wave.</p>
<p>In the standard case when <img class="math" src="../_images/math/df4eccc09cea597a3f814f8d9d58ae052961d6e6.png" alt="|a| < \infty, b = \infty"/>,
<a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a> works by evaluating the infinite series</p>
<div class="math">
<p><img src="../_images/math/cb1f62664824a1320f8446ac20a12e7d117094e7.png" alt="I = \int_a^{x_1} f(x) dx +
\sum_{k=1}^{\infty} \int_{x_k}^{x_{k+1}} f(x) dx" /></p>
</div><p>where <img class="math" src="../_images/math/f5d3dead9cb8d5f6327c737e014850be7a29d0f0.png" alt="x_k"/> are consecutive zeros (alternatively
some other periodic reference point) of <img class="math" src="../_images/math/c96dd6ec1dc4ad7520fbdc78fcdbec9edd068d0c.png" alt="f(x)"/>.
Accordingly, <a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a> requires information about the
zeros of <img class="math" src="../_images/math/c96dd6ec1dc4ad7520fbdc78fcdbec9edd068d0c.png" alt="f(x)"/>. For a periodic function, you can specify
the zeros by either providing the angular frequency <img class="math" src="../_images/math/54d7d48553f4d9e7ab418118607ea324cbfddfda.png" alt="\omega"/>
(<em>omega</em>) or the <em>period</em> <img class="math" src="../_images/math/07476a4927fd27f658fee67885db2c05fca2c298.png" alt="2 \pi/\omega"/>. In general, you can
specify the <img class="math" src="../_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/>-th zero by providing the <em>zeros</em> arguments.
Below is an example of each:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">mpmath</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span><span class="p">;</span> <span class="n">mp</span><span class="o">.</span><span class="n">pretty</span> <span class="o">=</span> <span class="bp">True</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sin</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
<span class="go">0.37833007080198</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">period</span><span class="o">=</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="mi">3</span><span class="p">)</span>
<span class="go">0.37833007080198</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">zeros</span><span class="o">=</span><span class="k">lambda</span> <span class="n">n</span><span class="p">:</span> <span class="n">pi</span><span class="o">*</span><span class="n">n</span><span class="o">/</span><span class="mi">3</span><span class="p">)</span>
<span class="go">0.37833007080198</span>
<span class="gp">>>> </span><span class="p">(</span><span class="n">ei</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">)</span><span class="o">-</span><span class="n">exp</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">*</span><span class="n">ei</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">))</span><span class="o">/</span><span class="mi">2</span> <span class="c"># Computed by Mathematica</span>
<span class="go">0.37833007080198</span>
</pre></div>
</div>
<p>Note that <em>zeros</em> was specified to multiply <img class="math" src="../_images/math/174fadd07fd54c9afe288e96558c92e0c1da733a.png" alt="n"/> by the
<em>half-period</em>, not the full period. In theory, it does not matter
whether each partial integral is done over a half period or a full
period. However, if done over half-periods, the infinite series
passed to <a title="mpmath.nsum" class="reference external" href="sums_limits.html#mpmath.nsum"><tt class="xref docutils literal"><span class="pre">nsum()</span></tt></a> becomes an <em>alternating series</em> and this
typically makes the extrapolation much more efficient.</p>
<p>Here is an example of an integration over the entire real line,
and a half-infinite integration starting at <img class="math" src="../_images/math/afac538740054d4786c170ef0123850adea5cdf9.png" alt="-\infty"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="n">inf</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="go">1.15572734979092</span>
<span class="gp">>>> </span><span class="n">pi</span><span class="o">/</span><span class="n">e</span>
<span class="go">1.15572734979092</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">,</span> <span class="p">[</span><span class="o">-</span><span class="n">inf</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">period</span><span class="o">=</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">)</span>
<span class="go">-0.0844109505595739</span>
<span class="gp">>>> </span><span class="n">cos</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span><span class="o">+</span><span class="n">si</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span><span class="o">-</span><span class="n">pi</span><span class="o">/</span><span class="mi">2</span>
<span class="go">-0.0844109505595738</span>
</pre></div>
</div>
<p>Of course, the integrand may contain a complex exponential just as
well as a real sine or cosine:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">exp</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">j</span><span class="o">*</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="n">inf</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
<span class="go">(0.156410688228254 + 0.0j)</span>
<span class="gp">>>> </span><span class="n">pi</span><span class="o">/</span><span class="n">e</span><span class="o">**</span><span class="mi">3</span>
<span class="go">0.156410688228254</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">exp</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">j</span><span class="o">*</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">+</span><span class="n">x</span><span class="o">+</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="p">[</span><span class="o">-</span><span class="n">inf</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
<span class="go">(0.00317486988463794 - 0.0447701735209082j)</span>
<span class="gp">>>> </span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">7</span><span class="p">)</span><span class="o">/</span><span class="n">exp</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="p">(</span><span class="n">j</span><span class="o">+</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">7</span><span class="p">))</span><span class="o">/</span><span class="mi">2</span><span class="p">)</span>
<span class="go">(0.00317486988463794 - 0.0447701735209082j)</span>
</pre></div>
</div>
<p><strong>Non-periodic functions</strong></p>
<p>If <img class="math" src="../_images/math/02aa2b5e1a4ae3d07fa845b04aa6ebdbc38d5ae3.png" alt="f(x) = g(x) h(x)"/> for some function <img class="math" src="../_images/math/1f745cf4286f88b236750e1e9a8abe3dacd9dc39.png" alt="h(x)"/> that is not
strictly periodic, <em>omega</em> or <em>period</em> might not work, and it might
be necessary to use <em>zeros</em>.</p>
<p>A notable exception can be made for Bessel functions which, though not
periodic, are “asymptotically periodic” in a sufficiently strong sense
that the sum extrapolation will work out:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">j0</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="n">period</span><span class="o">=</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">)</span>
<span class="go">1.0</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">j1</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="n">period</span><span class="o">=</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="p">)</span>
<span class="go">1.0</span>
</pre></div>
</div>
<p>More properly, one should provide the exact Bessel function zeros:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">j0zero</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">n</span><span class="p">:</span> <span class="n">findroot</span><span class="p">(</span><span class="n">j0</span><span class="p">,</span> <span class="n">pi</span><span class="o">*</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mf">0.25</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">j0</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="n">zeros</span><span class="o">=</span><span class="n">j0zero</span><span class="p">)</span>
<span class="go">1.0</span>
</pre></div>
</div>
<p>For an example where <em>zeros</em> becomes necessary, consider the
complete Fresnel integrals</p>
<div class="math">
<p><img src="../_images/math/227a94cc9f35ca20ce7e3efe31116de1dc3ac74d.png" alt="\int_0^{\infty} \cos x^2\,dx = \int_0^{\infty} \sin x^2\,dx
= \sqrt{\frac{\pi}{8}}." /></p>
</div><p>Although the integrands do not decrease in magnitude as
<img class="math" src="../_images/math/7d0eb486837b6e2d5c5e6f6a075ea764200db0c8.png" alt="x \to \infty"/>, the integrals are convergent since the oscillation
rate increases (causing consecutive periods to asymptotically
cancel out). These integrals are virtually impossible to calculate
to any kind of accuracy using standard quadrature rules. However,
if one provides the correct asymptotic distribution of zeros
(<img class="math" src="../_images/math/234066a9d63a3b077620a87c2e47f7ab5d421941.png" alt="x_n \sim \sqrt{n}"/>), <a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a> works:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">30</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">zeros</span><span class="o">=</span><span class="k">lambda</span> <span class="n">n</span><span class="p">:</span><span class="n">sqrt</span><span class="p">(</span><span class="n">pi</span><span class="o">*</span><span class="n">n</span><span class="p">))</span>
<span class="go">0.626657068657750125603941321203</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sin</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">zeros</span><span class="o">=</span><span class="k">lambda</span> <span class="n">n</span><span class="p">:</span><span class="n">sqrt</span><span class="p">(</span><span class="n">pi</span><span class="o">*</span><span class="n">n</span><span class="p">))</span>
<span class="go">0.626657068657750125603941321203</span>
<span class="gp">>>> </span><span class="n">sqrt</span><span class="p">(</span><span class="n">pi</span><span class="o">/</span><span class="mi">8</span><span class="p">)</span>
<span class="go">0.626657068657750125603941321203</span>
</pre></div>
</div>
<p>(Interestingly, these integrals can still be evaluated if one
places some other constant than <img class="math" src="../_images/math/f2ca003a7da0de4994b4733e203b74ff52d42553.png" alt="\pi"/> in the square root sign.)</p>
<p>In general, if <img class="math" src="../_images/math/ed3e589ab43eb026034d122e5cb99b8486b35fe4.png" alt="f(x) \sim g(x) \cos(h(x))"/>, the zeros follow
the inverse-function distribution <img class="math" src="../_images/math/45bcd9e1bf383713528adb77d14117fc8fd1942f.png" alt="h^{-1}(x)"/>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">mp</span><span class="o">.</span><span class="n">dps</span> <span class="o">=</span> <span class="mi">15</span>
<span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sin</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">zeros</span><span class="o">=</span><span class="k">lambda</span> <span class="n">n</span><span class="p">:</span> <span class="n">log</span><span class="p">(</span><span class="n">n</span><span class="p">))</span>
<span class="go">-0.25024394235267</span>
<span class="gp">>>> </span><span class="n">pi</span><span class="o">/</span><span class="mi">2</span><span class="o">-</span><span class="n">si</span><span class="p">(</span><span class="n">e</span><span class="p">)</span>
<span class="go">-0.250243942352671</span>
</pre></div>
</div>
<p><strong>Non-alternating functions</strong></p>
<p>If the integrand oscillates around a positive value, without
alternating signs, the extrapolation might fail. A simple trick
that sometimes works is to multiply or divide the frequency by 2:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">sin</span><span class="p">(</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="n">x</span><span class="o">**</span><span class="mi">4</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span> <span class="c"># Bad</span>
<span class="go">1.28642190869861</span>
<span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mf">0.5</span><span class="p">)</span> <span class="c"># Perfect</span>
<span class="go">1.28652953559617</span>
<span class="gp">>>> </span><span class="mi">1</span><span class="o">+</span><span class="p">(</span><span class="n">cos</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span><span class="o">+</span><span class="n">ci</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span><span class="o">+</span><span class="n">sin</span><span class="p">(</span><span class="mi">1</span><span class="p">))</span><span class="o">/</span><span class="mi">6</span>
<span class="go">1.28652953559617</span>
</pre></div>
</div>
<p><strong>Fast decay</strong></p>
<p><a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a> is primarily useful for slowly decaying
integrands. If the integrand decreases exponentially or faster,
<a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> will likely handle it without trouble (and generally be
much faster than <a title="mpmath.quadosc" class="reference internal" href="#mpmath.quadosc"><tt class="xref docutils literal"><span class="pre">quadosc()</span></tt></a>):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">quadosc</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="n">exp</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">],</span> <span class="n">omega</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="go">0.5</span>
<span class="gp">>>> </span><span class="n">quad</span><span class="p">(</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="p">)</span><span class="o">/</span><span class="n">exp</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="n">inf</span><span class="p">])</span>
<span class="go">0.5</span>
</pre></div>
</div>
</dd></dl>
</div>
<div class="section" id="quadrature-rules">
<h2>Quadrature rules<a class="headerlink" href="#quadrature-rules" title="Permalink to this headline">¶</a></h2>
<dl class="class">
<dt id="mpmath.calculus.quadrature.QuadratureRule">
<em class="property">class </em><tt class="descclassname">mpmath.calculus.quadrature.</tt><tt class="descname">QuadratureRule</tt><big>(</big><em>ctx</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule" title="Permalink to this definition">¶</a></dt>
<dd><p>Quadrature rules are implemented using this class, in order to
simplify the code and provide a common infrastructure
for tasks such as error estimation and node caching.</p>
<p>You can implement a custom quadrature rule by subclassing
<tt class="xref docutils literal"><span class="pre">QuadratureRule</span></tt> and implementing the appropriate
methods. The subclass can then be used by <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> by
passing it as the <em>method</em> argument.</p>
<p><tt class="xref docutils literal"><span class="pre">QuadratureRule</span></tt> instances are supposed to be singletons.
<tt class="xref docutils literal"><span class="pre">QuadratureRule</span></tt> therefore implements instance caching
in <tt class="xref docutils literal"><span class="pre">__new__()</span></tt>.</p>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.calc_nodes">
<tt class="descname">calc_nodes</tt><big>(</big><em>degree</em>, <em>prec</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.calc_nodes" title="Permalink to this definition">¶</a></dt>
<dd>Compute nodes for the standard interval <img class="math" src="../_images/math/5c3818b9565a33fd3aadba10026d32c5e3eea90f.png" alt="[-1, 1]"/>. Subclasses
should probably implement only this method, and use
<tt class="xref docutils literal"><span class="pre">get_nodes()</span></tt> method to retrieve the nodes.</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.clear">
<tt class="descname">clear</tt><big>(</big><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.clear" title="Permalink to this definition">¶</a></dt>
<dd>Delete cached node data.</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.estimate_error">
<tt class="descname">estimate_error</tt><big>(</big><em>results</em>, <em>prec</em>, <em>epsilon</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.estimate_error" title="Permalink to this definition">¶</a></dt>
<dd><p>Given results from integrations <img class="math" src="../_images/math/dbc47e6361fa8fb5c5677d7f82dfa0ee87591207.png" alt="[I_1, I_2, \ldots, I_k]"/> done
with a quadrature of rule of degree <img class="math" src="../_images/math/2f5f1a6ec7851c120425e8e3940b6414fc9b1247.png" alt="1, 2, \ldots, k"/>, estimate
the error of <img class="math" src="../_images/math/8ea972a3547d6e5b8c88bc77e0110bf571d75898.png" alt="I_k"/>.</p>
<p>For <img class="math" src="../_images/math/8e44ca2cae6be42148fe1b1fd8b5d55bd6b599bd.png" alt="k = 2"/>, we estimate <img class="math" src="../_images/math/e52229571ad29fa45fb776e4899957498f1b604b.png" alt="|I_{\infty}-I_2|"/> as <img class="math" src="../_images/math/e48888adf6a02ef922dd20ff328ffcb9fd3d266f.png" alt="|I_2-I_1|"/>.</p>
<p>For <img class="math" src="../_images/math/c1a21f8705d1c694796192dfe71e41cfad8fc07d.png" alt="k > 2"/>, we extrapolate <img class="math" src="../_images/math/1fd4fc966f94ca06f7b5a61a62a2a625d6bf5a32.png" alt="|I_{\infty}-I_k| \approx |I_{k+1}-I_k|"/>
from <img class="math" src="../_images/math/ff4e06f2942db25a00a9925670d43f2e0060fd42.png" alt="|I_k-I_{k-1}|"/> and <img class="math" src="../_images/math/53cb339d9c527c4d0809da8b5a3da266db1cd45c.png" alt="|I_k-I_{k-2}|"/> under the assumption
that each degree increment roughly doubles the accuracy of
the quadrature rule (this is true for both <tt class="xref docutils literal"><span class="pre">TanhSinh</span></tt>
and <tt class="xref docutils literal"><span class="pre">GaussLegendre</span></tt>). The extrapolation formula is given
by Borwein, Bailey & Girgensohn. Although not very conservative,
this method seems to be very robust in practice.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.get_nodes">
<tt class="descname">get_nodes</tt><big>(</big><em>a</em>, <em>b</em>, <em>degree</em>, <em>prec</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.get_nodes" title="Permalink to this definition">¶</a></dt>
<dd><p>Return nodes for given interval, degree and precision. The
nodes are retrieved from a cache if already computed;
otherwise they are computed by calling <tt class="xref docutils literal"><span class="pre">calc_nodes()</span></tt>
and are then cached.</p>
<p>Subclasses should probably not implement this method,
but just implement <tt class="xref docutils literal"><span class="pre">calc_nodes()</span></tt> for the actual
node computation.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.guess_degree">
<tt class="descname">guess_degree</tt><big>(</big><em>prec</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.guess_degree" title="Permalink to this definition">¶</a></dt>
<dd><p>Given a desired precision <img class="math" src="../_images/math/36f73fc1312ee0349b3f3a0f3bd9eb5504339011.png" alt="p"/> in bits, estimate the degree <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/>
of the quadrature required to accomplish full accuracy for
typical integrals. By default, <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> will perform up
to <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/> iterations. The value of <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/> should be a slight
overestimate, so that “slightly bad” integrals can be dealt
with automatically using a few extra iterations. On the
other hand, it should not be too big, so <a title="mpmath.quad" class="reference internal" href="#mpmath.quad"><tt class="xref docutils literal"><span class="pre">quad()</span></tt></a> can
quit within a reasonable amount of time when it is given
an “unsolvable” integral.</p>
<p>The default formula used by <tt class="xref docutils literal"><span class="pre">guess_degree()</span></tt> is tuned
for both <tt class="xref docutils literal"><span class="pre">TanhSinh</span></tt> and <tt class="xref docutils literal"><span class="pre">GaussLegendre</span></tt>.
The output is roughly as follows:</p>
<blockquote>
<table border="1" class="docutils">
<colgroup>
<col width="50%" />
<col width="50%" />
</colgroup>
<thead valign="bottom">
<tr><th class="head"><img class="math" src="../_images/math/36f73fc1312ee0349b3f3a0f3bd9eb5504339011.png" alt="p"/></th>
<th class="head"><img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/></th>
</tr>
</thead>
<tbody valign="top">
<tr><td>50</td>
<td>6</td>
</tr>
<tr><td>100</td>
<td>7</td>
</tr>
<tr><td>500</td>
<td>10</td>
</tr>
<tr><td>3000</td>
<td>12</td>
</tr>
</tbody>
</table>
</blockquote>
<p>This formula is based purely on a limited amount of
experimentation and will sometimes be wrong.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.sum_next">
<tt class="descname">sum_next</tt><big>(</big><em>f</em>, <em>nodes</em>, <em>degree</em>, <em>prec</em>, <em>previous</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.sum_next" title="Permalink to this definition">¶</a></dt>
<dd><p>Evaluates the step sum <img class="math" src="../_images/math/ca32b5510fe15c6ecc32b3a9e4f6d534e3030b29.png" alt="\sum w_k f(x_k)"/> where the <em>nodes</em> list
contains the <img class="math" src="../_images/math/94d71c036eb96eaaf04435264a5abb9b737ddb58.png" alt="(w_k, x_k)"/> pairs.</p>
<p><tt class="xref docutils literal"><span class="pre">summation()</span></tt> will supply the list <em>results</em> of
values computed by <tt class="xref docutils literal"><span class="pre">sum_next()</span></tt> at previous degrees, in
case the quadrature rule is able to reuse them.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.summation">
<tt class="descname">summation</tt><big>(</big><em>f</em>, <em>points</em>, <em>prec</em>, <em>epsilon</em>, <em>max_degree</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.summation" title="Permalink to this definition">¶</a></dt>
<dd><p>Main integration function. Computes the 1D integral over
the interval specified by <em>points</em>. For each subinterval,
performs quadrature of degree from 1 up to <em>max_degree</em>
until <tt class="xref docutils literal"><span class="pre">estimate_error()</span></tt> signals convergence.</p>
<p><tt class="xref docutils literal"><span class="pre">summation()</span></tt> transforms each subintegration to
the standard interval and then calls <tt class="xref docutils literal"><span class="pre">sum_next()</span></tt>.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.QuadratureRule.transform_nodes">
<tt class="descname">transform_nodes</tt><big>(</big><em>nodes</em>, <em>a</em>, <em>b</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.QuadratureRule.transform_nodes" title="Permalink to this definition">¶</a></dt>
<dd><p>Rescale standardized nodes (for <img class="math" src="../_images/math/5c3818b9565a33fd3aadba10026d32c5e3eea90f.png" alt="[-1, 1]"/>) to a general
interval <img class="math" src="../_images/math/da2e551d2ca2155b8d8f4935d2e9757722c9bab6.png" alt="[a, b]"/>. For a finite interval, a simple linear
change of variables is used. Otherwise, the following
transformations are used:</p>
<div class="math">
<p><img src="../_images/math/733647b1ab0cd443827387b1e630fa38a826acf0.png" alt="[a, \infty] : t = \frac{1}{x} + (a-1)
[-\infty, b] : t = (b+1) - \frac{1}{x}
[-\infty, \infty] : t = \frac{x}{\sqrt{1-x^2}}" /></p>
</div></dd></dl>
</dd></dl>
<div class="section" id="tanh-sinh-rule">
<h3>Tanh-sinh rule<a class="headerlink" href="#tanh-sinh-rule" title="Permalink to this headline">¶</a></h3>
<dl class="class">
<dt id="mpmath.calculus.quadrature.TanhSinh">
<em class="property">class </em><tt class="descclassname">mpmath.calculus.quadrature.</tt><tt class="descname">TanhSinh</tt><big>(</big><em>ctx</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.TanhSinh" title="Permalink to this definition">¶</a></dt>
<dd><p>This class implements “tanh-sinh” or “doubly exponential”
quadrature. This quadrature rule is based on the Euler-Maclaurin
integral formula. By performing a change of variables involving
nested exponentials / hyperbolic functions (hence the name), the
derivatives at the endpoints vanish rapidly. Since the error term
in the Euler-Maclaurin formula depends on the derivatives at the
endpoints, a simple step sum becomes extremely accurate. In
practice, this means that doubling the number of evaluation
points roughly doubles the number of accurate digits.</p>
<dl class="docutils">
<dt>Comparison to Gauss-Legendre:</dt>
<dd><ul class="first last simple">
<li>Initial computation of nodes is usually faster</li>
<li>Handles endpoint singularities better</li>
<li>Handles infinite integration intervals better</li>
<li>Is slower for smooth integrands once nodes have been computed</li>
</ul>
</dd>
</dl>
<p>The implementation of the tanh-sinh algorithm is based on the
description given in Borwein, Bailey & Girgensohn, “Experimentation
in Mathematics - Computational Paths to Discovery”, A K Peters,
2003, pages 312-313. In the present implementation, a few
improvements have been made:</p>
<blockquote>
<ul class="simple">
<li>A more efficient scheme is used to compute nodes (exploiting
recurrence for the exponential function)</li>
<li>The nodes are computed successively instead of all at once</li>
</ul>
</blockquote>
<p>Various documents describing the algorithm are available online, e.g.:</p>
<blockquote>
<ul class="simple">
<li><a class="reference external" href="http://crd.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf">http://crd.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf</a></li>
<li><a class="reference external" href="http://users.cs.dal.ca/~jborwein/tanh-sinh.pdf">http://users.cs.dal.ca/~jborwein/tanh-sinh.pdf</a></li>
</ul>
</blockquote>
<dl class="method">
<dt id="mpmath.calculus.quadrature.TanhSinh.calc_nodes">
<tt class="descname">calc_nodes</tt><big>(</big><em>degree</em>, <em>prec</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.TanhSinh.calc_nodes" title="Permalink to this definition">¶</a></dt>
<dd><p>The abscissas and weights for tanh-sinh quadrature of degree
<img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/> are given by</p>
<div class="math">
<p><img src="../_images/math/aeae3571408b66ee4c92fbfbf2675bd8fc2dc6c1.png" alt="x_k = \tanh(\pi/2 \sinh(t_k))
w_k = \pi/2 \cosh(t_k) / \cosh(\pi/2 \sinh(t_k))^2" /></p>
</div><p>where <img class="math" src="../_images/math/df81febf50c1ed419c6f7273a3b5f4b066b47d4c.png" alt="t_k = t_0 + hk"/> for a step length <img class="math" src="../_images/math/18b44cbd4fe661b938ee214a993203fb20e604bd.png" alt="h \sim 2^{-m}"/>. The
list of nodes is actually infinite, but the weights die off so
rapidly that only a few are needed.</p>
</dd></dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.TanhSinh.sum_next">
<tt class="descname">sum_next</tt><big>(</big><em>f</em>, <em>nodes</em>, <em>degree</em>, <em>prec</em>, <em>previous</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.TanhSinh.sum_next" title="Permalink to this definition">¶</a></dt>
<dd>Step sum for tanh-sinh quadrature of degree <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/>. We exploit the
fact that half of the abscissas at degree <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/> are precisely the
abscissas from degree <img class="math" src="../_images/math/96eec32de93b032b6d531a03aa5b346156b62a75.png" alt="m-1"/>. Thus reusing the result from
the previous level allows a 2x speedup.</dd></dl>
</dd></dl>
</div>
<div class="section" id="gauss-legendre-rule">
<h3>Gauss-Legendre rule<a class="headerlink" href="#gauss-legendre-rule" title="Permalink to this headline">¶</a></h3>
<dl class="class">
<dt id="mpmath.calculus.quadrature.GaussLegendre">
<em class="property">class </em><tt class="descclassname">mpmath.calculus.quadrature.</tt><tt class="descname">GaussLegendre</tt><big>(</big><em>ctx</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.GaussLegendre" title="Permalink to this definition">¶</a></dt>
<dd><p>This class implements Gauss-Legendre quadrature, which is
exceptionally efficient for polynomials and polynomial-like (i.e.
very smooth) integrands.</p>
<p>The abscissas and weights are given by roots and values of
Legendre polynomials, which are the orthogonal polynomials
on <img class="math" src="../_images/math/5c3818b9565a33fd3aadba10026d32c5e3eea90f.png" alt="[-1, 1]"/> with respect to the unit weight
(see <a title="mpmath.legendre" class="reference external" href="../functions/orthogonal.html#mpmath.legendre"><tt class="xref docutils literal"><span class="pre">legendre()</span></tt></a>).</p>
<p>In this implementation, we take the “degree” <img class="math" src="../_images/math/f5047d1e0cbb50ec208923a22cd517c55100fa7b.png" alt="m"/> of the quadrature
to denote a Gauss-Legendre rule of degree <img class="math" src="../_images/math/79ed1d781e90030a3e1bc71b11f5b05262bccde5.png" alt="3 \cdot 2^m"/> (following
Borwein, Bailey & Girgensohn). This way we get quadratic, rather
than linear, convergence as the degree is incremented.</p>
<dl class="docutils">
<dt>Comparison to tanh-sinh quadrature:</dt>
<dd><ul class="first last simple">
<li>Is faster for smooth integrands once nodes have been computed</li>
<li>Initial computation of nodes is usually slower</li>
<li>Handles endpoint singularities worse</li>
<li>Handles infinite integration intervals worse</li>
</ul>
</dd>
</dl>
<dl class="method">
<dt id="mpmath.calculus.quadrature.GaussLegendre.calc_nodes">
<tt class="descname">calc_nodes</tt><big>(</big><em>degree</em>, <em>prec</em>, <em>verbose=False</em><big>)</big><a class="headerlink" href="#mpmath.calculus.quadrature.GaussLegendre.calc_nodes" title="Permalink to this definition">¶</a></dt>
<dd>Calculates the abscissas and weights for Gauss-Legendre
quadrature of degree of given degree (actually <img class="math" src="../_images/math/79ed1d781e90030a3e1bc71b11f5b05262bccde5.png" alt="3 \cdot 2^m"/>).</dd></dl>
</dd></dl>
</div>
</div>
</div>
</div>
</div>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<h3><a href="../index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference external" href="#">Numerical integration (quadrature)</a><ul>
<li><a class="reference external" href="#standard-quadrature-quad">Standard quadrature (<tt class="docutils literal"><span class="pre">quad</span></tt>)</a></li>
<li><a class="reference external" href="#oscillatory-quadrature-quadosc">Oscillatory quadrature (<tt class="docutils literal"><span class="pre">quadosc</span></tt>)</a></li>
<li><a class="reference external" href="#quadrature-rules">Quadrature rules</a><ul>
<li><a class="reference external" href="#tanh-sinh-rule">Tanh-sinh rule</a></li>
<li><a class="reference external" href="#gauss-legendre-rule">Gauss-Legendre rule</a></li>
</ul>
</li>
</ul>
</li>
</ul>
<h4>Previous topic</h4>
<p class="topless"><a href="differentiation.html"
title="previous chapter">Differentiation</a></p>
<h4>Next topic</h4>
<p class="topless"><a href="odes.html"
title="next chapter">Ordinary differential equations</a></p>
<h3>This Page</h3>
<ul class="this-page-menu">
<li><a href="../_sources/calculus/integration.txt"
rel="nofollow">Show Source</a></li>
</ul>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="../search.html" method="get">
<input type="text" name="q" size="18" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
>index</a></li>
<li class="right" >
<a href="../modindex.html" title="Global Module Index"
>modules</a> |</li>
<li class="right" >
<a href="odes.html" title="Ordinary differential equations"
>next</a> |</li>
<li class="right" >
<a href="differentiation.html" title="Differentiation"
>previous</a> |</li>
<li><a href="../index.html">mpmath v0.17 documentation</a> »</li>
<li><a href="index.html" >Numerical calculus</a> »</li>
</ul>
</div>
<div class="footer">
© Copyright 2010, Fredrik Johansson.
Last updated on Feb 06, 2011.
Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 0.6.6.
</div>
</body>
</html>
distrib
>
Fedora
>
13
>
i386
>
media
>
updates
>
by-pkgid
>
f3eb4c16ba6256fe5a10e54bf649f01f
>
files
>
1221
