iconEuler Home

Iteration and Sequences

This notebook demonstrates iteration in Euler.

Our first example is the growth of an amount of money at interest
rate of 5%.
>q=1.05; iterate("x*q",1000,n=10)'
               1050 
             1102.5 
           1157.625 
         1215.50625 
       1276.2815625 
      1340.09564063 
      1407.10042266 
      1477.45544379 
      1551.32821598 
      1628.89462678 
Let us subtract the amount of 30 in each year, while adding interests,
and plot the result for 30 years.
>plot2d(iterate("x*q-30",1000,n=30)); insimg;

09 - Introduction to Iteration

Iterating like this is a simple form of recursion,

09 - Introduction to Iteration

For more complicated sequences, we use the function "sequence". This
function computes the value x[n] from all previous values
x[1],...,x[n-1] given in a vector. Addtionally the variable n can be
used.

Here is an example computing the Fibonacci numbers.

09 - Introduction to Iteration

>sequence("x[n-1]+x[n-2]",[1,1],10)
[ 1  1  2  3  5  8  13  21  34  55 ]
The n-th root of the Fibonacci numbers converges to the golden
ratio.
>plot2d(sequence("x[n-1]+x[n-2]",[1,1],500)^(1/(1:500))); insimg;

09 - Introduction to Iteration

We can compute more complicate recursions. In the following example,
the n-th element is the sum of all previous elements. This time, we
start with a scalar value 1. Of course, the result is 2^n.
>sequence("sum(x)",1,10)
[ 1  1  2  4  8  16  32  64  128  256 ]
Instead of an expression in x and n, we can also use a function.

In the following example, we iterate

09 - Introduction to Iteration

where A is a Matrix, and each x[n] is a 2x1 vector.
>A=[1,1;1,2]; function step(x,n) := A.x[,n-1]
>sequence("step",[1;1],5)
Column 1 to 3:
                  1                   2                   5 
                  1                   3                   8 
Column 4 to 5:
                 13                  34 
                 21                  55 
There is the function "matrixpower" which yields the same result. This
function is faster, since it uses only log_2(n) matrix
multiplications.
>matrixpower(A,4).[1;1]
                 34 
                 55 
Of course, we could also use a loop in a function, or in the command
line.
>x=[1;1]; loop 1 to 4; x=A.x; end; x,
                 34 
                 55 

Convergence

Often, we want to iterate until convergence. If the function does not
converge, the user has to press the escape key to stop the iteration!
>iterate("cos(x)",1)
0.739085133216
If such an iteration converges, the limit is a solution of

09 - Introduction to Iteration

This examples does also work with interval arithmetic.
>res := iterate("cos(x)",~1,2~)
~0.7390851332164,0.739085133218~
If we make res a little larger, we are able to prove that h contains
a fixed point of cos.
>h=expand(res,100), cos(h) << h
~0.73908513314,0.7390851333~
1
Or we could set up a function for the iteration.
>function f(x) := (x+2/x)/2
The iteration x(n+1)=f(x(n)) converges quickly to
sqrt(2).
>iterate("f",2)
1.41421356237
It might be easier to use an expression in x instead
of a function.
>iterate("(x+2/x)/2",2)
1.41421356237
The function "iterate" with the parameter "n=..." returns the
history of 5 iterations.
>iterate("f",2,5,n=5)
[ 1.5  1.41666666667  1.41421568627  1.41421356237  1.41421356237 ]
Starting with an interval does not work.
>niterate("f",~1,2~,5)
[ ~1,2~  ~1,2~  ~1,2~  ~1,2~  ~1,2~ ]
The reason is that the interval arithmetic is too approximate. We can
improve the evaluation of the expression using a sub-division of the
interval.
>function s(x) := ieval("(x+2/x)/2",x,10)
Then we can iterate until the result is optimal, and the interval does
not get smaller. The result is an interval, which must contain a
solution of

09 - Introduction to Iteration

And the only solution of this is 

09 - Introduction to Iteration

>iterate("s",~1,2~)
~1.414213562371,1.414213562375~
The function "iterate" works for vectors too. We try the
arithmetic-geometric mean

09 - Introduction to Iteration

The n-th value is stored in the column vector x[n].
>function g(x) := [(x[1]+x[2])/2;sqrt(x[1]*x[2])]
Using the above function, the iteration converges to the arithmetic
geometric mean of the start values.
>iterate("g",[1;2])
      1.45679103105 
      1.45679103105 
The convergence is rather quick.
>iterate("g",[1;2],4,n=6)
Column 1 to 3:
                1.5       1.45710678119       1.45679104815 
      1.41421356237       1.45647531512       1.45679101394 
Column 4:
      1.45679103105 
      1.45679103105 
Again, we can use an interval iteration. It shows that the iteration
is stable. It does not prove, that the limit is in the computed
bounds, however.
>iterate("g",[~1~;~2~],4,n=5)
Column 1:
  ~1.4999999999999991,1.5000000000000013~ 
  ~1.4142135623730943,1.4142135623730958~ 
Column 2:
   ~1.4571067811865461,1.457106781186549~ 
  ~1.4564753151219691,1.4564753151219716~ 
Column 3:
  ~1.4567910481542572,1.4567910481542607~ 
  ~1.4567910139395532,1.4567910139395566~ 
Column 4:
   ~1.4567910310469048,1.456791031046909~ 
   ~1.4567910310469048,1.456791031046909~ 
The following iteration converges very slowly.

09 - Introduction to Iteration

>iterate("sqrt(x)",2,10)
[ 1.41421356237  1.189207115  1.09050773267  1.04427378243
1.02189714865  1.01088928605  1.00542990111  1.00271127505
1.00135471989  1.00067713069 ]
Applying the Steffenson accelarator, speeds up the convergence.
>steffenson("sqrt(x)",2,10)
[ 1.04887809808  1.00027970075  1.00000000978  1 ]

Tables of Functions

There is an interesting way to generate a sequence using Maxima
expressions. The mxmtable command will by default print and plot
the sequence, and it returns the sequence as a column vector.

For an example we generate the n-th derivatives of x^x at 1.
>mxmtable("diffat(x^x,x=1,n)","n",1,8,frac=1);
        1 
        2 
        3 
        8 
       10 
       54 
      -42 
      944 

09 - Introduction to Iteration

Euler Home