2
$\begingroup$

For the closed Newton-Cotes quadrature over $[x_1, x_n]$, the coefficients $H_{n,i}$ for $$ \int_{x_1}^{x_n} f(x)\:\text{d}x = h \sum_{i=1}^n H_{n,i} \; f(x_i) $$ are given explicitly by $$ H_{n,r+1} =\frac{(-1)^{n-r}}{r!(n-r)!}\int_0^n \frac{\prod_{k=0}^n (t-k)}{t-r}\:\text{d}t; $$ see http://mathworld.wolfram.com/Newton-CotesFormulas.html.

To compute those values numerically, one could go ahead and evaluate the integral, but perhaps that's not the most efficient thing to do. SciPy does something else, but I don't quite get it.

Any hints?

$\endgroup$
  • 1
    $\begingroup$ Is your question why SciPy does something different? Are they just using a different equation that, however, evaluates to the same number? $\endgroup$ – Wolfgang Bangerth Oct 28 '16 at 15:22
  • $\begingroup$ They are the same numbers, but I don't understand what they are doing (any why). No reference given in the source code. $\endgroup$ – Nico Schlömer Oct 28 '16 at 15:25
2
$\begingroup$

I can explain what the scipy code is doing (https://github.com/scipy/scipy/blob/v0.18.0/scipy/integrate/quadrature.py#L833-L842). I haven't checked this directly, but this should be correct in principle.

The line ti = 2 * yi - 1 remaps the interval $(0,1)$ to $(-1,1)$, so now the weights will be the integrals of Lagrange elementary interpolants over $(-1,1)$ with nodes $(-1:2/N:1)$.

The line C = ti ** nvec[:, np.newaxis] constructs the (transpose of) the Vandermonde matrix, with entries $C_{ij} = t_i ^ j$. This matrix has the property that when multiplied on the left by a vector $\alpha$, it produces $$\alpha^\top C = (p_\alpha(t_0), \ldots, p_\alpha(t_N))$$ which is the vector of $p_\alpha$ (polynomial with coefficients being the entries of $\alpha$) evaluated at each $t_0,\ldots,t_N$.

The inverse of this, $C^{-1}$, will then consist of vectors $\alpha_k$ that give the coefficients of Lagrange interpolants, because $\alpha_k^\top C = e_k$, and having $p_{\alpha_k}(t_l) = \delta_{kl}$ is the definition of Lagrange interpolants.

The Vandermonde matrix is very ill-conditioned, so scipy applies twice the Newton-Raphson iteration for the (matrix) reciprocal: Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv).

Once the coefficients of Lagrange interpolants, $C^{-1}$, are known, it is straightforward to integrate them, by multiplying $C^{-1}$ on the right with $(2/1, 0, 2/3, 0, 2/5, \ldots)^\top$ (the integrals of $x^m$ for odd $m$ are zero—the integral is over $(-1,1)$ because of ti above), which gives the Newton-Cotes weights: ai = Cinv[:, ::2].dot(vec) * (N / 2.)

$\endgroup$
  • 2
    $\begingroup$ You should propose a pull request that adds your explanation to the documentation of that function. $\endgroup$ – Wolfgang Bangerth Oct 28 '16 at 20:37
0
$\begingroup$

In the pass I did something similarly. At the time I applied Gauss-Legendre quadrature to evaluate $H_{n,r+1}$ intergals; this calculation adopted (1) Gauss-Legendre weight and points calculated in advance (then, change the variable $t$ in the integrand to $\xi \in [-1,1]$), and (2) multi-precision floating point library fmlib.

I guess that you obtain $H_{n,i} $ from Lagrange Interpolation Formulae. Let me review some basic thing so that you may recognize the same derivation. From the Lagrange Interpolation, we have \begin{align} f(x) = \sum_{i=1}^{n} l_i (x) f_i + R_n (x) \end{align} where \begin{align} l_i (x) &= \frac{ \Pi_n( x) }{ (x-x_i) } . \frac{ 1 }{ \Pi _n ' (x_i) } = \frac{ ( x - x_1 ) \ldots ( x - x_{i-1} ) ( x - x_{i+1} ) \ldots ( x - x_{n} ) }{ ( x_i - x_1 ) \ldots ( x_i - x_{i-1} ) ( x_i - x_{i+1} ) \ldots ( x_i - x_{n} ) } , \\ R_n (x) &= \frac{ f^{(n)} (\xi) }{ n! } \cdot \Pi _n (x) , \end{align} for some $\xi = \xi(x) \in (x_1, x_n )$. Note that we defined $ \Pi_n (x) = \prod_{i=1}^{n} (x-x_i) $, and we also have \begin{align} \left| R_n (x) \right| &\le \frac{ ( x_n - x_1 )^{n} }{ n! } \max_{ x_1 \le x \le x_n } \left| f^{(n)} (x) \right| \end{align} which can be proved using the Role’s theorem.

Now we take \begin{align} \int_{x_1}^{x_n} f(x) d x &= \sum_{i=1}^{n} \underbrace{ \left( \int_{x_1}^{x_n} l_i (x) d x \right) }_{ A_ i } f_i + \underbrace{ \int_{x_1}^{x_n} R_n (x) d x }_{ B } \end{align} Define \begin{align} A_ i &= \int_{x_1}^{x_n} l_i (x) d x ~~~(\equiv H_{n,i} ? ) \\ B &= \int_{x_1}^{x_n} R_n (x) d x \end{align} We are going to estimate those integrals. Change variable, \begin{align} h &= \frac{ x_n - x_1 }{ n - 1 } ,~~~~~ x = x_1 + (t-1) h ,~~ \forall t \in [1,n] \end{align} Using the latter one, $x - x_j = x_1 + (t-1)h - x_j $ and $ x_j = x_1 + (j-1) h $, we have $$ x - x_j = x_1 + (t-1)h - \left[ x_1 + (j-1)h \right] = (t-j)h $$ Thus \begin{align} A_ i &= \int_{x_1}^{x_n} l_i (x) d x = \int_{x_1}^{x_n} \frac{ ( x - x_1 ) \ldots ( x - x_{i-1} ) ( x - x_{i+1} ) \ldots ( x - x_{n} ) }{ ( x_i - x_1 ) \ldots ( x_i - x_{i-1} ) ( x_i - x_{i+1} ) \ldots ( x_i - x_{n} ) } d x \\ &= \int_{1}^{n} \frac{ h ( t - 1 ) \cdot h ( t - 2 ) \cdots h ( t-i+1 ) \cdot h (t-i-1) \cdots h(t-n) }{ h ( i - 1 ) \cdot h ( i - 2 ) \cdots h ( i-i+1 ) \cdot h (i-i-1) \cdots h(i-n) } \cdot h d t \end{align} Or \begin{align} A_ i &= h \int_{1}^{n} \frac{ ( t - 1 ) \cdot ( t - 2 ) \cdots ( t-i+1 ) \cdot (t-i-1) \cdots (t-n) }{ ( i - 1 ) \cdot ( i - 2 ) \cdots ( i-i+1 ) \cdot (i-i-1) \cdots (i-n) } d t \end{align} And \begin{align} B &= \int_{x_1}^{x_n} R_n (x) d x = \int_{x_1}^{x_n} ( x - x_1 ) \ldots ( x - x_{i} ) \ldots ( x - x_{n} ) \cdot \frac{ f^{(n)} (\xi (x) ) }{ n! } d x \\ &= \int_{1}^{n} h ( t - 1 ) \cdot h ( t - 2 ) \cdots h ( t - i) \cdots h(t-n) \cdot \frac{ f^{(n)} (\xi (x(t)) ) }{ n! } \cdot h d t \end{align} We estimate \begin{align} | B | &\le M h^{n+1} \end{align} where $M = |C_n | \sup_{ \xi \in [x_1,x_n] } \left| f^{(n)} (\xi) \right| $ and \begin{align} C_n &= \frac{1}{ n! } \int_{1}^{n} ( t - 1 ) \cdots ( t - i) \cdots ( t - n ) d t \label{c} \end{align} Here integrals $A_i$ and $C_n$ can be estimated EXACTLY using Gauss-Legendre formulas with $n/2$ points.

As I could see in numerical result at the time, we have (a) $A_i $ are increased to very large values as $n$ increasing, for example $\max A_i \sim 10^{0}, 10^{2}, 10^{4}, 10^{7}, 10^{10}$ and $10^{24}$ when $n=10,20,30,40,50$, and $100$, respectively; (b) $A_i$ is symmetry with repect to $i$, i.e. $A_{i} = A_{n-i+1}$; (c) $\sum_{i=1}^{n} A_{i} = n-1$; and (c) $C_n = 0$ for $n$ odd, and $C_n$ is bounded by a small value when $n$ even.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.