problem setting

• evaluation of

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

decomposition into a series

• term decomposition

• polynomial devision

• integration of partial terms

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

• recursion relation

• 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

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$

• 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

• 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

• fixpoints: $y_n\equiv y^*$

stable/unstable fixpoints: recursion does/does not converge

• $y^*=y^{(0)} = 0$

• $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

• 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);
   }
}