250 lines
17 KiB
HTML
250 lines
17 KiB
HTML
|
|
<!DOCTYPE html>
|
|
|
|
<html xmlns="http://www.w3.org/1999/xhtml">
|
|
<head>
|
|
<meta charset="utf-8" />
|
|
<title>numerical.integration — pylib 2019.5.19 documentation</title>
|
|
<link rel="stylesheet" href="../../_static/alabaster.css" type="text/css" />
|
|
<link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
|
|
<link rel="stylesheet" href="../../_static/custom.css" type="text/css" />
|
|
<script type="text/javascript" id="documentation_options" data-url_root="../../" src="../../_static/documentation_options.js"></script>
|
|
<script type="text/javascript" src="../../_static/jquery.js"></script>
|
|
<script type="text/javascript" src="../../_static/underscore.js"></script>
|
|
<script type="text/javascript" src="../../_static/doctools.js"></script>
|
|
<script type="text/javascript" src="../../_static/language_data.js"></script>
|
|
<script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
|
|
<link rel="index" title="Index" href="../../genindex.html" />
|
|
<link rel="search" title="Search" href="../../search.html" />
|
|
|
|
<link rel="stylesheet" href="../../_static/custom.css" type="text/css" />
|
|
|
|
|
|
<meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />
|
|
|
|
</head><body>
|
|
|
|
|
|
<div class="document">
|
|
<div class="documentwrapper">
|
|
<div class="bodywrapper">
|
|
|
|
|
|
<div class="body" role="main">
|
|
|
|
<h1>Source code for numerical.integration</h1><div class="highlight"><pre>
|
|
<span></span><span class="c1"># -*- coding: utf-8 -*-</span>
|
|
<span class="sd">"""Numerical integration, numerical quadrature.</span>
|
|
|
|
<span class="sd">de: numerische Integration, numerische Quadratur.</span>
|
|
|
|
<span class="sd">:Date: 2015-10-15</span>
|
|
|
|
<span class="sd">.. module:: integration</span>
|
|
<span class="sd"> :platform: *nix, Windows</span>
|
|
<span class="sd"> :synopsis: Numerical integration.</span>
|
|
|
|
<span class="sd">.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de></span>
|
|
<span class="sd">"""</span>
|
|
|
|
<span class="kn">from</span> <span class="nn">__future__</span> <span class="k">import</span> <span class="n">division</span>
|
|
<span class="kn">from</span> <span class="nn">numpy</span> <span class="k">import</span> <span class="n">linspace</span><span class="p">,</span> <span class="n">trapz</span><span class="p">,</span> <span class="n">zeros</span>
|
|
|
|
<div class="viewcode-block" id="trapez"><a class="viewcode-back" href="../../numerical.html#numerical.integration.trapez">[docs]</a><span class="k">def</span> <span class="nf">trapez</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">N</span><span class="o">=</span><span class="mi">10</span><span class="p">,</span> <span class="n">x</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
|
|
<span class="n">save_values</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
|
|
<span class="sa">r</span><span class="sd">"""</span>
|
|
<span class="sd"> Integration of :math:`f(x)` using the trapezoidal rule</span>
|
|
<span class="sd"> (Simpson's rule, Kepler's rule).</span>
|
|
|
|
<span class="sd"> de: Trapezregel, Simpsonregel (Thomas Simpson), Keplersche</span>
|
|
<span class="sd"> Fassregel (Johannes Kepler)</span>
|
|
|
|
<span class="sd"> :param f: function to integrate.</span>
|
|
<span class="sd"> :type f: function or list</span>
|
|
<span class="sd"> :param a: lower limit of integration (default = 0).</span>
|
|
<span class="sd"> :type a: float</span>
|
|
<span class="sd"> :param b: upper limit of integration (default = 1).</span>
|
|
<span class="sd"> :type b: float</span>
|
|
<span class="sd"> :param N: specify the number of subintervals.</span>
|
|
<span class="sd"> :type N: int</span>
|
|
<span class="sd"> :param x: variable of integration, necessary if f is a list</span>
|
|
<span class="sd"> (default = None).</span>
|
|
<span class="sd"> :type x: list</span>
|
|
<span class="sd"> :param verbose: print information (default = False)</span>
|
|
<span class="sd"> :type verbose: bool</span>
|
|
|
|
<span class="sd"> :returns: the definite integral as approximated by trapezoidal</span>
|
|
<span class="sd"> rule.</span>
|
|
<span class="sd"> :rtype: float</span>
|
|
|
|
<span class="sd"> The trapezoidal rule approximates the integral by the area of a</span>
|
|
<span class="sd"> trapezoid with base h=b-a and sides equal to the values of the</span>
|
|
<span class="sd"> integrand at the two end points.</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> f_n(x) = f(a)+\frac{f(b)-f(a)}{b-a}(x-a)</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> I &= \int\limits_a^b f(x) \,\mathrm{d}x \\</span>
|
|
<span class="sd"> I &\approx \int\limits_a^b f_n(x) \,\mathrm{d}x \\</span>
|
|
<span class="sd"> &= \int\limits_a^b</span>
|
|
<span class="sd"> \left( f(a)+\frac{f(b)-f(a)}{b-a}(x-a) \right)</span>
|
|
<span class="sd"> \mathrm{d}x \\</span>
|
|
<span class="sd"> &= \left.\left( f(a)-a\frac{f(b)-f(a)}{b-a} \right)</span>
|
|
<span class="sd"> x \right\vert_a^b +</span>
|
|
<span class="sd"> \left. \frac{f(b)-f(a)}{b-a} \frac{x^2}{2}</span>
|
|
<span class="sd"> \right\vert_a^b \\</span>
|
|
<span class="sd"> &= \frac{b-a}{2}\left[f(a)+f(b)\right]</span>
|
|
|
|
<span class="sd"> The composite trapezium rule. If the interval is divided into n</span>
|
|
<span class="sd"> segments (not necessarily equal)</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> a = x_0 \leq x_1 \leq x_2 \leq \ldots \leq x_n = b</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> I &\approx \sum\limits_{i=0}^{n-1} \frac{1}{2} (x_{i+1}-x_i)</span>
|
|
<span class="sd"> \left[f(x_{i+1})+f(x_i)\right] \\</span>
|
|
|
|
<span class="sd"> Special Case (Equaliy spaced base points)</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> x_{i+1}-x_i = h \quad \forall i</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> I \approx h \left\{ \frac{1}{2} \left[f(x_0)+f(x_n)\right] +</span>
|
|
<span class="sd"> \sum\limits_{i=1}^{n-1} f(x_i) \right\}</span>
|
|
|
|
<span class="sd"> .. rubric:: Example</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> I &= \int\limits_a^b f(x) \,\mathrm{d}x \\</span>
|
|
<span class="sd"> f(x) &= x^2 \\</span>
|
|
<span class="sd"> a &= 0 \\</span>
|
|
<span class="sd"> b &= 1</span>
|
|
<span class="sd"> </span>
|
|
<span class="sd"> analytical solution</span>
|
|
|
|
<span class="sd"> .. math::</span>
|
|
<span class="sd"> I = \int\limits_{0}^{1} x^2 \,\mathrm{d}x</span>
|
|
<span class="sd"> = \left. \frac{1}{3} x^3 \right\vert_0^1</span>
|
|
<span class="sd"> = \frac{1}{3}</span>
|
|
<span class="sd"> </span>
|
|
<span class="sd"> numerical solution</span>
|
|
|
|
<span class="sd"> >>> f = lambda(x): x**2</span>
|
|
<span class="sd"> >>> trapez(f, 0, 1, 1)</span>
|
|
<span class="sd"> 0.5</span>
|
|
<span class="sd"> >>> trapez(f, 0, 1, 10)</span>
|
|
<span class="sd"> 0.3350000000000001</span>
|
|
<span class="sd"> >>> trapez(f, 0, 1, 100)</span>
|
|
<span class="sd"> 0.33335000000000004</span>
|
|
<span class="sd"> """</span>
|
|
<span class="n">N</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">N</span><span class="p">)</span>
|
|
<span class="c1"># f is function or list</span>
|
|
<span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="s1">'__call__'</span><span class="p">):</span>
|
|
<span class="c1"># h width of each subinterval</span>
|
|
<span class="n">h</span> <span class="o">=</span> <span class="p">(</span><span class="n">b</span><span class="o">-</span><span class="n">a</span><span class="p">)</span><span class="o">/</span><span class="n">N</span>
|
|
|
|
<span class="c1"># x variable of integration</span>
|
|
<span class="n">x</span> <span class="o">=</span> <span class="n">linspace</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
|
|
<span class="k">if</span> <span class="n">save_values</span><span class="p">:</span>
|
|
<span class="c1"># ff contribution from the points</span>
|
|
<span class="n">ff</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">))</span>
|
|
<span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">N</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
|
|
<span class="n">ff</span><span class="p">[</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">])</span>
|
|
<span class="n">T</span> <span class="o">=</span> <span class="p">(</span><span class="n">ff</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">/</span><span class="mf">2.</span><span class="o">+</span><span class="nb">sum</span><span class="p">(</span><span class="n">ff</span><span class="p">[</span><span class="mi">1</span><span class="p">:</span><span class="n">N</span><span class="p">])</span><span class="o">+</span><span class="n">ff</span><span class="p">[</span><span class="n">N</span><span class="p">]</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span><span class="o">*</span><span class="n">h</span>
|
|
<span class="k">else</span><span class="p">:</span>
|
|
<span class="n">TL</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
|
|
<span class="n">TR</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">N</span><span class="p">])</span>
|
|
<span class="n">TI</span> <span class="o">=</span> <span class="mi">0</span>
|
|
<span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">N</span><span class="p">):</span>
|
|
<span class="n">TI</span> <span class="o">=</span> <span class="n">TI</span> <span class="o">+</span> <span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">])</span>
|
|
<span class="n">T</span> <span class="o">=</span> <span class="p">(</span><span class="n">TL</span><span class="o">/</span><span class="mf">2.</span><span class="o">+</span><span class="n">TI</span><span class="o">+</span><span class="n">TR</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span><span class="o">*</span><span class="n">h</span>
|
|
<span class="k">else</span><span class="p">:</span>
|
|
<span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">f</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span>
|
|
<span class="n">T</span> <span class="o">=</span> <span class="mi">0</span>
|
|
<span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
|
|
<span class="n">T</span> <span class="o">=</span> <span class="n">T</span> <span class="o">+</span> <span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">])</span><span class="o">/</span><span class="mi">2</span><span class="o">*</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">+</span><span class="n">f</span><span class="p">[</span><span class="n">n</span><span class="p">])</span>
|
|
|
|
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
|
|
<span class="nb">print</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>
|
|
|
|
<span class="k">return</span> <span class="n">T</span></div>
|
|
|
|
|
|
<span class="k">if</span> <span class="vm">__name__</span> <span class="o">==</span> <span class="s1">'__main__'</span><span class="p">:</span>
|
|
<span class="n">func</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">x</span><span class="o">**</span><span class="mi">2</span>
|
|
<span class="n">trapez</span><span class="p">(</span><span class="n">func</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="mf">1e6</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
|
|
<span class="c1">#print(trapz(func, linspace(0,1,10)))</span>
|
|
|
|
<span class="n">trapez</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="mi">4</span><span class="p">,</span><span class="mi">9</span><span class="p">],</span> <span class="n">x</span><span class="o">=</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="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
|
|
<span class="c1">#print(trapz([0,1,4,9]))</span>
|
|
|
|
<span class="n">trapez</span><span class="p">([</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">x</span><span class="o">=</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="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
|
|
|
|
<span class="n">trapez</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="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">x</span><span class="o">=</span><span class="p">[</span><span class="o">-</span><span class="mi">1</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="mi">1</span><span class="p">],</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
|
|
</pre></div>
|
|
|
|
</div>
|
|
|
|
</div>
|
|
</div>
|
|
<div class="sphinxsidebar" role="navigation" aria-label="main navigation">
|
|
<div class="sphinxsidebarwrapper">
|
|
<h1 class="logo"><a href="../../index.html">pylib</a></h1>
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
<h3>Navigation</h3>
|
|
|
|
<div class="relations">
|
|
<h3>Related Topics</h3>
|
|
<ul>
|
|
<li><a href="../../index.html">Documentation overview</a><ul>
|
|
<li><a href="../index.html">Module code</a><ul>
|
|
</ul></li>
|
|
</ul></li>
|
|
</ul>
|
|
</div>
|
|
<div id="searchbox" style="display: none" role="search">
|
|
<h3 id="searchlabel">Quick search</h3>
|
|
<div class="searchformwrapper">
|
|
<form class="search" action="../../search.html" method="get">
|
|
<input type="text" name="q" aria-labelledby="searchlabel" />
|
|
<input type="submit" value="Go" />
|
|
</form>
|
|
</div>
|
|
</div>
|
|
<script type="text/javascript">$('#searchbox').show(0);</script>
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
</div>
|
|
</div>
|
|
<div class="clearer"></div>
|
|
</div>
|
|
<div class="footer">
|
|
©2019, Daniel Weschke.
|
|
|
|
|
|
|
Powered by <a href="http://sphinx-doc.org/">Sphinx 2.1.2</a>
|
|
& <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
|
|
|
|
</div>
|
|
|
|
|
|
|
|
|
|
</body>
|
|
</html> |