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




Claudius Gros, WS 2019/20

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

Summation, Recursion and Stability

problem setting




$$ y_n\ =\ \int_0^1\, \frac{x^n}{x+a}\, dx$$

decomposition into a series

$$ y_n\ =\ \int_0^1\, \frac{x^n}{x+a}\, dx$$
$$ \frac{x^n}{x+a} \,=\, \frac{x^n-(-a)^n+(-a)^n}{x+a} \,=\, \frac{x^n-(-a)^n}{x-(-a)} \,+\, \frac{(-a)^n}{x+a} $$ $$ \frac{x^n-(-a)^n}{x-(-a)} \,=\, x^{n-1} \,+\, (-a)x^{n-2} \,+\, (-a)^2x^{n-3} \,+\, \ \dots\ \,+\, (-a)^{n-2}x \,+\, (-a)^{n-1} $$ $$ y_n \ =\ \frac{1}{n} \,+\,\frac{(-a)}{n-1} \,+\,\frac{(-a)^2}{n-2} \,+\, \ \dots\ \,+\,\frac{(-a)^{n-2}}{2} \,+\,\frac{(-a)^{n-1}}{1} \ +\ (-a)^n\,\log\left(\frac{1+a}{a}\right) $$

careful when subtracting large terms

n, y_n (up/down):     0     0.0953101798043249     0.0953101798043249
n, y_n (up/down):     1     0.0468982019567507     0.0468982019567507
n, y_n (up/down):     2     0.0310179804324928     0.0310179804324928
n, y_n (up/down):     3     0.0231535290083968     0.0231535290084016
n, y_n (up/down):     4     0.0184647099159747     0.0184647099160125
n, y_n (up/down):     5     0.0153529008402984     0.0153529008389668
n, y_n (up/down):     6     0.0131376582721714     0.0131376582624891
n, y_n (up/down):     7     0.0114805602934212     0.0114805602325932
n, y_n (up/down):     8     0.0101943966001272     0.0101943983716508
n, y_n (up/down):     9     0.0091671049594879     0.0091671236256787
n, y_n (up/down):    10     0.0083286762237549     0.0083287335918189
n, y_n (up/down):    11     0.0076236724853516     0.0076222290336939
n, y_n (up/down):    12     0.0071258544921875     0.0071148651228605
n, y_n (up/down):    13     0.0057373046875000     0.0057442655768129
n, y_n (up/down):    14     0.0156250000000000     0.0142329159691716
n, y_n (up/down):    15    -0.0937500000000000    -0.0717638696291832
n, y_n (up/down):    16     0.8750000000000000     0.7797726530809349
n, y_n (up/down):    17    -6.0000000000000000    -7.9918313757219790
n, y_n (up/down):    18    80.0000000000000000    67.9816784898641800
n, y_n (up/down):    19  -512.0000000000000000  -807.3893324901815000
n, y_n (up/down):    20  4096.0000000000000000  6540.9419220573470000

relative errors

straightforward recursion

$$ y_n\,=\,\int_0^1\,\frac{x^n}{x+a}\, dx \,=\,\int_0^1\,\frac{x^{n-1}(x+a-a)}{x+a}\,dx \,=\,\int_0^1\,x^{n-1}\,dx \,-\,a\,\int_0^1\,\frac{x^{n-1}}{x+a}\,dx $$ $$ y_n \ =\ \frac{1}{n}\,-\,a\,y_{n-1}, \qquad\quad y_0 = \log\left(\frac{a+1}{a}\right) $$
  n      forward recursion  

  0     0.0953101798043249  
  1     0.0468982019567507  
  2     0.0310179804324935 
  3     0.0231535290083985 
  4     0.0184647099160155 
  5     0.0153529008398455 
  6     0.0131376582682119
  7     0.0114805601750236 
  8     0.0101943982497636 
  9     0.0091671286134746 
 10     0.0083287138652537
 11     0.0076219522565537 
 12     0.0071138107677960 
 13     0.0057849692451174 
 14     0.0135788789773972 
 15    -0.0691221231073049
 16     0.7537212310730486 
 17    -7.4783887813187210  
 18    74.8394433687427600  
 19  -748.3418021084802000 
 20  7483.4680210848010000

stability analyis of recursion formula

$$ y_n\ =\ -a\,y_{n-1} \,+\,\frac{1}{n} $$

small perturbation
consider time evolution of two closed-by trajetories: $\ y_n,\ \tilde y_n \qquad y_n-\tilde y_n=\epsilon_n\ll 1$

$$ \epsilon_n \ =\ -a\,\epsilon_{n-1} $$

dynamical system theory
dynamical systems with positive (negative) Lyapunov exponents
are denoted chaotic (regular)

forward vs backward recursion

$$ y_n\ =\ -a\,y_{n-1} \,+\,\frac{1}{n}, \qquad \qquad y_{n-1} \ =\ \frac{1}{a}\,\left(\frac{1}{n}-y_n\right) $$
  n      forward recursion     backward recursion

  0     0.0953101798043249     0.0953101798043249
  1     0.0468982019567507     0.0468982019567514
  2     0.0310179804324935     0.0310179804324860
  3     0.0231535290083985     0.0231535290084733
  4     0.0184647099160155     0.0184647099152671
  5     0.0153529008398455     0.0153529008473289
  6     0.0131376582682119     0.0131376581933773
  7     0.0114805601750236     0.0114805609233700
  8     0.0101943982497636     0.0101943907661000
  9     0.0091671286134746     0.0091672034481114
 10     0.0083287138652537     0.0083279655188863
 11     0.0076219522565537     0.0076294357202278
 12     0.0071138107677960     0.0070389761310555
 13     0.0057849692451174     0.0065333156125223
 14     0.0135788789773972     0.0060954153033482
 15    -0.0691221231073049     0.0057125136331850
 16     0.7537212310730486     0.0053748636681501
 17    -7.4783887813187210     0.0050748927302638
 18    74.8394433687427600     0.0048066282529171
 19  -748.3418021084802000     0.0045652964181972

summation/recursion code

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>       // for printf
#include <math.h>        // for pow
 
using namespace std;
 
//   y_n = \int_0^1 dx x^n/(x+a)
 
// ---
// --- Straight summing up of integrated terms
// ---
double performSumUp(int N, double a) 
{
  double result = pow(-a,1.0*N)*log((1+a)/a);
  for (int i=0; i<N; i++)
    result +=  pow(-a,1.0*i)/(N-i);
  return result;
 }    //  end of summationUP
 
// ---
// --- Straight summing down of integrated terms
// ---
double performSumDown(int N, double a) 
{
  double result = pow(-a,1.0*N)*log((1+a)/a);
  for (int i=N-1; i>=0; i--)
    result += pow(-a,1.0*i)/(N-i);
  return result;
}    //  end of summationDown
 
// ---
// --- Solving the forward recursion
// ---
double forwardRecursion(int N, double a) 
{
  double result = log((1+a)/a);
  for (int i=1; i<=N; i++)
    result =  -a*result + 1.0/i;
  return result;
}    //  end of forwardRecursion
 
// ---
// --- Solving the backward recursion
// ---
double backwardRecursion(int N, double a) 
{
  int startN = 50;
  double result = 0.0;
  for (int i=startN; i>N; i--)
    result = (1.0/i-result)/a;
  return result;
}    //  end of forwardRecursion
 
// ---
// --- main
// ---
int main() 
{
  int maxN = 21;
  double a = 10.0;
  printf("\na = %6.1f\n\n",a);
// *** ****************** ***
// *** summation of terms ***
// *** ****************** ***
  for (int n=0; n<maxN; n+=2)
    {
    double resultSumUp   = performSumUp(n,a);
    double resultSumDown = performSumDown(n,a);
    printf("n, y_n (up/down):   %3d %22.16f %22.16f\n",
                             n,resultSumUp,resultSumDown);
    }
  printf("\n");
// *** ****************** ***
// *** recursion relation ***
// *** ****************** ***
  printf("%3s %22s %22s\n","n","forward recursion","backward recursion");
  for (int n=0; n<maxN; n+=2)
    {
    double resultForward  = forwardRecursion(n,a);
    double resultBackward = backwardRecursion(n,a);
    printf("%3d %22.16f %22.16f\n",
                      n,resultForward,resultBackward);
    }
//
 return 1;
}  

fixpoints of maps

$$ y_{n+1} \ =\ gy_n(1-y_n), \qquad\quad y_{n+2} \ =\ g^2y_n(1-y_n)(1-gy_n(1-y_n)) $$ $$ y^*=0, \quad\qquad 1=g(1-y^*)\qquad \Longleftrightarrow\qquad y^*=1-1/g $$

stable/unstable fixpoints: recursion does/does not converge


$$ y_{n+1} \ \approx\ gy_n + O(y_n^2), \qquad\quad \left|\frac{y_{n+1}}{y_n}\right| = g \ \ \to\ \ \left\{ \begin{array}{rl} <1 & (g<1) \\ >1 & (g>1) \end{array} \right. $$ $$ \left|\frac{y_{n+1}-y^{(1)}}{y_n-y^{(1)}}\right| \ \approx \ g\left|1-2y^{(1)}\right| \ = \ \left|2-g\right| \ \ \to\ \ \left\{ \begin{array}{rl} <1 & (g\in[1,3]) \\ >1 & (g>3) \end{array} \right. $$
  

finding roots via iterative bisection




Copy Copy to clipboad
Downlaod Download
/* compile with g++  -Wall -O3 -o roots.e roots.c
                           -O3             : optimization
                     -Wall                 : all Warnings 
 */
#include <iostream>  // input/output handling
#include <stdio.h>   // for formatted input/output
#include <stdlib.h>  // srand, rand, exit
#include <cmath>     // for log
using namespace std;
 
const int nRoots  = 30000000;        // number of roots of polynomial
const int maxIter = 10000;           // maximal number of bisections
const double accuracy  = 1.0E-12;    // desired accuracy
 
// ---
// --- a polynomial function function
// ---
double polynomial(double x) 
{
 double result = 1.0;
 for (int iRoots=0; iRoots<nRoots; iRoots++)
   result = result*(x-iRoots);
 return result;
}   
 
// ---
// --- start of main
// ---
int main()
{
 srand(time(NULL));   // random number initiatlization
 double xLower = nRoots*(rand()/(double)RAND_MAX);
 double xUpper = xLower + 1.0;
 double xMean  = 0.5*(xUpper+xLower);
 double absError = xUpper-xLower;
//
 printf("# starting xLower/xUpper: %10.6f %10.6f\n", xLower, xUpper);
//
 double fLower, fUpper, fMean;
 int nIter = 0;
 while ( (absError>accuracy*nRoots)&&(nIter<maxIter) )
   {
   nIter++;
   fLower = polynomial(xLower);
   fUpper = polynomial(xUpper);
   fMean  = polynomial(xMean);
//
   if ( fLower*fUpper > 1.0 )
     {
     printf(" ! error, no bisection possible %10.6f %10.6f \n",xLower, xUpper);
     exit(EXIT_FAILURE);                  // abort program
     }
//
   if ( fMean*fUpper > 0.0 )
     xUpper = xMean;
   else
     xLower = xMean;
//
   xMean  = 0.5*(xUpper+xLower);
   absError = xUpper-xLower;
//
   printf("%7d %24.12f %16.12f %16.12f \n", nIter, xMean, absError,log(absError));
   }
}