Advanced Introduction to C++, Scientific Computing and Machine Learning




Claudius Gros, WS 2019/20

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

Numerical Evaluation of Integrals

numerical integration




$$ I\ =\ \int_a^b\, f(x)\, dx $$
$$ \int_{x_0}^{x_1} f(x)dx \ \approx\ \frac{f(x_0)+f(x_1)}{2} \left(x_1-x_0\right) $$


$$ T_n\ \equiv\ \sum_{k=1}^n\, \frac{f(x_{k-1})+f(x_{k})}{2}\,\Delta x $$

trapezoidal rule at work

$$ y_m\ =\ \int_0^1\, \frac{x^m}{x+a}\, dx$$
       n     p                   T(p)        (4T(p)-T(p-1)/3
       1     0     0.0454545454545455
       2     1     0.0227273635533981     0.0151516362530156
       4     2     0.0114620139299330     0.0077068973887780
       8     3     0.0066417103644638     0.0050349425093074
      16     4     0.0051140844181988     0.0046048757694438
      32     5     0.0047045044781142     0.0045679778314194
      64     6     0.0046002267297100     0.0045654674802419
     128     7     0.0045740370560200     0.0045653071647900
     256     8     0.0045674820820493     0.0045652970907257
     512     9     0.0045658428656951     0.0045652964602438
    1024    10     0.0045654330320428     0.0045652964208253
    2048    11     0.0045653305717818     0.0045652964183615
    4096    12     0.0045653049566010     0.0045652964182075
    8192    13     0.0045652985527986     0.0045652964181978
   16384    14     0.0045652969518476     0.0045652964181972
                             ^                             ^
                             9                            16

convergence scaling



$$ I \ =\ \int_a^b\, f(x)\, dx\ \approx\ T_n \,+\, O(n^{-2}) $$
$$ T_n\ =\ \frac{\Delta x}{2}\,\sum_{k=1}^n\, \Big(f(x_{k-1})+f(x_{k})\Big), \qquad \Delta x\,=\, \frac{b-a}{n} $$
$$ \begin{array}{rcl} f(x_0+\delta x) & =& f(x_0) + f'(x_0)\delta x + O\big(\delta x^2\big), \qquad\quad (\mbox{general}) \\ f(x_1) & =& f(x_0) + f'(x_0)\Delta x + O\big(\Delta x^2\big), \qquad\quad \Delta x = x_1-x_0 \end{array} $$ $$ \begin{array}{rcl} T_n^{(1)} & =& \frac{\Delta x}{2}\Big(f(x_0)+f(x_1)\Big) \ =\ \frac{\Delta x}{2}\Big(f(x_0)+f(x_0)+f'(x_0)\Delta x +O\big(\Delta x^2\big)\Big) \\ & =& f(x_0)\Delta x+f'(x_0)\frac{\Delta x^2}{2} +O\big(\Delta x^3\big) \end{array} $$ $$ \begin{array}{rcl} I^{(1)} & =& \int_{x_0}^{x_1} f(x)dx \ =\ \int_{0}^{\Delta x}dx' \Big(f(x_0)+f'(x_0)x'+O\big((x')^2\big)\Big) \\ & =& f(x_0)\Delta x+f'(x_0)\frac{\Delta x^2}{2} +O\big(\Delta x^3\big) \end{array} $$ $$ I^{(1)}-T_n^{(1)} \ \propto\ \big(\Delta x\big)^3, \qquad\quad I-T_n \ \propto\ n\big(I^{(1)}-T_n^{(1)}\big) \ \propto\ \big(\Delta x\big)^2, \qquad\quad \Delta x \propto 1/n $$

convergence acceleration









$$ T_n\ =\ I \,+\,\frac{c}{n^2} \,+\, \,O(n^{-3}),\quad\qquad T_{2n}\ =\ I \,+\,\frac{c}{4n^2} \,+\, \,O(n^{-3}) $$ $$ \frac{4}{3}T_{2n}-\frac{1}{3}T_{n} \ =\ \left(\frac{4}{3}-\frac{1}{3}\right)\, I \,+\, \left(\frac{4}{3}\frac{1}{4}-\frac{1}{3}\right)\, \frac{c}{n^2} \,+\, \,O(n^{-3}) $$

a single term

$$ T_n\ \to\ \Delta x\frac{f_0+f_1}{2},\qquad\quad T_{2n}\ \to\ \frac{\Delta x}{2}\Big( \frac{f_0+f_{1/2}}{2}+\frac{f_{1/2}+f_1}{2}\Big) ,\qquad\quad $$ $$ S_n\ \to\ \frac{\Delta x}{3}\Big(f_0+2f_{1/2}+f_1 -\frac{f_0+f_1}{2}\Big) \ = \frac{\Delta x}{6}\Big(f_0+4f_{1/2}+f_1\Big) $$

Simpson rule - error scaling

Compare

$$ S_n^{(1)}\ = \frac{\Delta x}{6}\Big(f_0+4f_{1/2}+f_1\Big), \qquad\quad I^{(1)} \ =\ \int_{x_0}^{x_1} f(x) dx $$

Romberg method

$$ A_n\ =\ A + \frac{c}{n^z},\qquad\quad A_{2n} \ =\ A + \frac{1}{2^z}\frac{c}{n^z},\qquad\quad aA_{2n}-bA_n\ =\ (a-b)A + \Big(\frac{a}{2^z}-b\Big)\frac{c}{n^z} $$

Romberg iteration

$$ S_n \ =\ \frac{4}{3}T_{2n}-\frac{1}{3}T_{n} \ =\ I \,+\,O(n^{-4}) $$

accelerated integration

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>      // for printf
#include <math.h>       // for pow
 
using namespace std;
 
/** Integrating
  * y_m = \int_0^1 x^m/(x+a)dx
  * for various m and n of the trapezoidal integration method
  */
 
// ---
// --- Straight summing up trapezoids
// ---
double sumTrapezoids(int n, int m, double a) 
{
 double result = 0.0;
 double x0, x1;
 double DeltaX = 1.0/n;
 for (int i=0; i<n; i++)
   {
   x0 = i*DeltaX;
   x1 = (i+1)*DeltaX;
   result += 0.5*( pow(x0,1.0*m)/(x0+a) + pow(x1,1.0*m)/(x1+a) );
   }
 return result*DeltaX;
}    //  end of sumTrapezoids
 
// ---
// --- main
// ---
int main() 
{
 int    m = 19;     
 double a = 10.0;
// *** ****************************** ***
// *** loop over numbers of trapzoids ***
// *** ****************************** ***
 double     Rn = 0.0;
 double     Sn = 0.0;
 double lastSn = 0.0;
 double     Tn = 0.0;
 double lastTn = 0.0;
 int         n = 1;
 printf("%8s %5s %22s %22s %22s\n", "n","p","T(p)","(4T(p)-T(p-1)/3","(16S(p)-S(p-1)/15");
 for (int jj=0; jj<15; jj++)
   {
   Tn = sumTrapezoids(n, m, a);
   Sn = ( 4.0*Tn-lastTn)/3.0;
   Rn = (16.0*Sn-lastSn)/15.0;
   if (jj==0)
     printf("%8d %5d %22.16f\n", n, jj, Tn);
   if (jj==1)
     printf("%8d %5d %22.16f %22.16f\n", n, jj, Tn, Sn);
   if (jj>1)
     printf("%8d %5d %22.16f %22.16f %22.16f\n", n, jj, Tn, Sn, Rn);
//
   n = n*2;
   lastTn = Tn;
   lastSn = Sn;
   }
  printf("\n");
//
 return 1;
}