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

Claudius Gros, SS 2023

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

# problem setting

• evaluation of

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

• by decomposition and summation
• by recursion and inverse recursion
• for various $\ n\$ and $\quad a\gg1$

# decomposition into a series

$$y_n\ =\ \int_0^1\, \frac{x^n}{x+a}\, dx$$
• term decomposition
$$\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}$$
• polynomial devision
$$\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}$$
• integration of partial terms
$$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

• summing up and down yields different results for   $n>2$
• substraction of very large terms leads to large relative errors
• $a=10$
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

• numbers/measurements $\ a_\alpha,\ b_\beta\$
with small relative errors $\ \alpha,\ \beta$
• $$a_\alpha = a\,(1\pm\alpha) \quad\quad b_\beta = b\,(1\pm\beta) \quad\quad a,b,\alpha,\beta>0 \quad\quad \alpha,\beta\ll0 \qquad$$
• substracting two large numbers $\ a,b\gg 0$
• $$a_\alpha-b_\beta = (a-b)\,\pm\,(a\,\alpha+b\,\beta) = (a-b)\left(1\,\pm\,\frac{a\,\alpha+b\,\beta}{|a-b|}\right)$$
• may lead to very large relative errors
• $$\frac{a\,\alpha+b\,\beta}{|a-b|}$$
• for $\ a\approx b$

# 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$$
• recursion relation
$$y_n \ =\ \frac{1}{n}\,-\,a\,y_{n-1}, \qquad\quad y_0 = \log\left(\frac{a+1}{a}\right)$$
• recursion is unstable
  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}$$
• unstable (exponentially) for $\quad |-a|>1$
• Lyapunov exponent $\ \lambda>0$
• $$e^{\lambda} \,\equiv\, |-a|$$
• error growth exponentially with number $n$ of iterations
(butterfly effect)
• $$\epsilon_n = e^{n\lambda}\epsilon_0$$

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)$$
• error $\ \epsilon_n\ll 1$
• $$y_n\ \to\ y_n\pm\epsilon_n, \qquad\quad y_{n-1}\pm\epsilon_{n-1} = \frac{1}{a}\,\left(\frac{1}{n}-y_n\pm\epsilon_n\right)$$
• stable for $\quad a>1$
• $$\epsilon_{n-1} = \frac{\epsilon_n}{a}$$
• starting with $\quad y_{50}=0$
  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

• all four methods for evaluating $\ y_n = \int_0^1 \frac{x^n}{x+a}dx$
• examine the dependence of the results on the value of $\ a$
• what happens for $\ a\in\,]0,1[\$ ?
#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

• recursions are equivalent to fixpoints of maps
• example logistic map
• $y_n\in[0,1]$ : population density in year $n$
• $g\in[0,4]$ : reproduction rate per year
• $(1-y_n)$ : resource limitation
$$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))$$
• fixpoints: $y_n\equiv y^*$
$$y^*=0, \quad\qquad 1=g(1-y^*)\qquad \Longleftrightarrow\qquad y^*=1-1/g$$

### stable/unstable fixpoints: recursion does/does not converge

• $y^*=y^{(0)} = 0$
$$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.$$
• $y^*=y^{(1)} = 1-1/g\quad$ for $\quad g>1\quad$ (since $y>0$) $$\epsilon_{n+1} \ =\ y_{n+1}-y^{(1)} \ \approx \ g\left(1-2y^{(1)}\right)\left(y_n-y^{(1)}\right) \ = \ g\left(1-2y^*\right)\epsilon_n$$ which leads to
$$\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.$$
• fixpoints for $g\in[1,3]$
• $y^{(0)}=0$: unstable/stable forward/backward iteration
• $y^{(1)}=1-1/g$: stable/unstable forward/backward iteration

# finding roots via iterative bisection

• find root of $\qquad\qquad f(x)\,=\,0$
• start on both sides $$\mathrm{sign}(f(x_{lower}))\,=\, -\mathrm{sign}(f(x_{upper})) \qquad\qquad \qquad\qquad \qquad\qquad \qquad\qquad$$
• consider midpoint $$x_{mid}\,=\, \frac{x_{lower}+x_{upper}}{2} \qquad\qquad \mathrm{sign}(f(x_{mid})) \,=\, ? \qquad\qquad \qquad\qquad \qquad\qquad \qquad\qquad \qquad\qquad \qquad\qquad$$
• exponential convergence
test with gnuplot
here: $\qquad\qquad\qquad f(x)=\prod_i (x-x_i)$

/* 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 initialization
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 %20.16f \n", nIter, xMean, absError, absError/xMean);
}
}