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;
Iterating like this is a simple form of recursion,
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.
>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;
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
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
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
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
And the only solution of this is
>iterate("s",~1,2~)
~1.414213562371,1.414213562375~
The function "iterate" works for vectors too. We try the arithmetic-geometric mean
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.
>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 ]
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