iconEuler Examples

Introduction to Matlab in Euler

by R. Grothmann

This is a translation of an introduction to Matlab by Cleve Moler. I
found the file via Googe on the following address.

http://www.mathworks.com/moler/intro.pdf

I try to do the same in Euler, so you can study the differences and
similarities.

So the content of this notebook and the ideas are Moler's.

Golden Ratio

He begins by exploring the Golden Ratio, which is in fact the zero of
the polynomial

Molders Introduction to Matlab in Euler

>p=[-1,-1,1]
[ -1  -1  1 ]
The Matlab function for polynomial roots is roots(p).
>r=polysolve(p)
[ -0.61803398875+0i   1.61803398875+0i  ]
Of course, this can be solved exactly with Maxima. In Matlab with the
symbolic toolbox, this would be a Maple solution.
>&solve(1/x=x-1)
                       1 - sqrt(5)      sqrt(5) + 1
                  [x = -----------, x = -----------]
                            2                2

Maxima prints nicer than Matlab. However, it even better using Latex.

Molders Introduction to Matlab in Euler

>real(r[2])
1.61803398875
Let us fix the value for phi.
>phi=(sqrt(5)+1)/2
1.61803398875
Maxima and Matlab can evaluate this to an arbitrary number of digits.
>&bfloat((sqrt(5)+1)/2)
                 1.6180339887498948482045868343656b0

The computation becomes arbitrarily slow, if you do this.

The following one line function can be done in Matlab in an ugly way
only.
>function f(x) := 1/x-(x-1)
>plot2d("f",0,4);
>solve("f",1)
1.61803398875
>plot2d(phi,0,>points,>add); insimg;

Molders Introduction to Matlab in Euler

The following functions plots a sketch in Euler. The code is a
translation from the Matlab code. In any case, it is not a good idea
to use such programs for sketches.
>function plotgr ...
phi = (1+sqrt(5))/2;
x = [0,phi,phi,0,0];
y = [0,0,1,1,0];
u = [1,1];
v = [0,1];
plot2d(x,y,a=0,b=1.7,c=-0.7,d=1,<frame,<axis,<grid);
plot2d(u,v,style="--",>add);
ctext("phi",toscreen([phi/2,1.1]));
ctext("phi-1",toscreen([(1+phi)/2,-.05]));
ctext("1",toscreen([-.05,.5]));
ctext("1",toscreen([.5,-.05]));
endfunction
>plotgr; insimg(crop=0.7);

Molders Introduction to Matlab in Euler

This is an abuse of Euler or Matlab. A graphical program like C.a.R.
can do a nicer image with much less effort.

Molders Introduction to Matlab in Euler

A continued Fraction

The next example is a continued fraction. Euler and Matlab can
construct a string in the following way.
>p="1"; loop 1 to 6; p="1+1/("+p+")"; end; p,
1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))
The string can be evaluated.
>p()
1.61538461538
This string is called a continued fraction. The fraction converges to
the Golden Ratio phi.
>phi-p()
0.00264937336528
The following function is a direct translation of Matlab code. It
determines the numerator and denominator of the continued fraction.
>function test(n) ...
  p = 1;
  q = 1;
  for k = 1:n
    s = p;
    p = p + q;
    q = s;
  end
  return{p,q}
endfunction
>{p,q}=test(6);
>p+"/"+q
21/13
This is an approximation of phi.

Euler has the contfrac function, which produces the continued fraction
of a value. It also has the fracprint() function, which prints a value
in fractional form (using continued fractions, by the way).
>contfrac((1+sqrt(5))/2,6), contfracval(%), fracprint(%)
[ 1  1  1  1  1  1  1 ]
1.61538461538
21/13
The following is a well known approximation of pi. It has been known
since ancient times.
>contfrac(pi,3), api=contfracbest(pi,3), fracprint(api), api-pi
[ 3  7  15  1 ]
3.14159292035
355/113
2.66764189405e-007

Fibonacci Numbers

The following is a direct translation of the Matlab code.
>function fibonacci (n) ...
  f=zeros(1,n);
  f[1] = 1;
  f[2] = 2;
  for k = 3:n
    f[k] = f[k-1] + f[k-2];
  end
  return f
endfunction
>fibonacci(12)
[ 1  2  3  5  8  13  21  34  55  89  144  233 ]
Euler has the sequence command to produce this.
>sequence("x[n-2]+x[n-1]",[1,2],12)
[ 1  2  3  5  8  13  21  34  55  89  144  233 ]
The following is the famous inefficient recursion.
>function fibnum (n) ...
  if n<=1 then return 1
  else return fibnum(n-1)+fibnum(n-2)
endfunction
It works for small n. But the effort grows exponentially like 2^n.
>tic; fibnum(24), toc;
75025

Used 0.858 seconds
The quotients of successive Fibonacci numbers converges to phi.
>n=40; f=fibonacci(n);
>longestformat; (f[2:n]/f[1:n-1]-phi)', longformat;
     0.3819660112501051 
    -0.1180339887498949 
    0.04863267791677184 
   -0.01803398874989481 
   0.006966011250105098 
  -0.002649373365279484 
   0.001013630297724166 
 -0.0003869299263654646 
  0.0001478294319232631 
-5.646066000730698e-005 
  2.15668056606777e-005 
-8.237676933475768e-006 
 3.146528619657474e-006 
-1.201864648914253e-006 
 4.590717870289751e-007 
-1.753497695933248e-007 
 6.697765919660981e-008 
-2.558318845657936e-008 
 9.771908393574336e-009 
-3.732536946188247e-009 
 1.425702222945802e-009 
-5.445699446937624e-010 
 2.080071670462758e-010 
-7.945177848966978e-011 
 3.034772433352373e-011 
-1.159183860011126e-011 
 4.427569422205124e-012 
-1.691313755713964e-012 
 6.459277557269161e-013 
-2.466915560717098e-013 
 9.414691248821328e-014 
-3.597122599785507e-014 
 1.376676550535194e-014 
-5.329070518200751e-015 
 1.998401444325282e-015 
-8.881784197001252e-016 
 2.220446049250313e-016 
-2.220446049250313e-016 
                      0 
Indeed there is the famous formula for the Fibonacci numbers.

Molders Introduction to Matlab in Euler

>n=1:40; 1/(2*phi-1)*(phi^(n+1)-(1-phi)^(n+1))
[ 1  2  3  5  8  13  21  34  55  89  144  233  377  610  987  1597  2584
4181  6765  10946  17711  28657  46368  75025  121393  196418  317811
514229  832040  1346269  2178309  3524578  5702887  9227465  14930352
24157817  39088169  63245986  102334155  165580141 ]

Fractal Fern

This is a fractal generated by randomly selecting one of four
transformations of the plane, and plotting the iterated points. The
code is a tranlation of Matlab code.
>function fern (n) ...
p = [ .85, .92, .99, 1.00];
A1 = [ .85, .04; -.04, .85]; b1 = [0; 1.6];
A2 = [ .20, -.26; .23, .22]; b2 = [0; 1.6];
A3 = [-.15, .28; .26, .24]; b3 = [0; .44];
A4 = [ 0, 0 ; 0, .16];
x=[0.5;0.5];
X=zeros(n,2);
r=random(1,n);
loop 1 to n;
  if r[#] < p[1]
    x = A1.x + b1;
  elseif r[#] < p[2]
    x = A2.x + b2;
  elseif r[#] < p[3]
    x = A3.x + b3;
  else
    x = A4.x;
  endif
  X[#]=x';
end;
return X;
endfunction
We generate 100 thousand points.
>X=fern(100000)'; ...
>plot2d(X[1],X[2],a=-3,b=3,c=0,d=10, ...
>  >points,<frame,<grid,style=".",color=green,>insimg);

Molders Introduction to Matlab in Euler

It is clear, that neither Matlab nor Euler is good for type of
iterations. In fact, a multi-core algorithm could generate this
fractal much faster. Also it should really be done on the level of
machine code.
>remvalue X;

Magic Squares

Matlab had magic squares since the beginning. It is not clear why,
since the program is essentially a numerical program like Euler.
>shortformat; A=magic(3)
            8             1             6 
            3             5             7 
            4             9             2 
>sum(A)'
[ 15  15  15 ]
>sum(A')'
[ 15  15  15 ]
>sum(diag(A,0))
15
>sum(diag(flipy(A),0))
15
The sums must clearly be equal to the following.
>sum(1:9)/3
15
Like in Matlab, we can produce lots of permuted magic squares in the
following way.
>for k=0:3; rot(A,k), "", end;
            8             1             6 
            3             5             7 
            4             9             2 

            6             7             2 
            1             5             9 
            8             3             4 

            2             9             4 
            7             5             3 
            6             1             8 

            4             3             8 
            9             5             1 
            2             7             6 

>for k=0:3; rot(A',k), "", end;
            8             3             4 
            1             5             9 
            6             7             2 

            4             9             2 
            3             5             7 
            8             1             6 

            2             7             6 
            9             5             1 
            4             3             8 

            6             1             8 
            7             5             3 
            2             9             4 

This magic square is regular.
>det(A)
-360
>inv(A)
     0.147222     -0.144444     0.0638889 
   -0.0611111     0.0222222      0.105556 
   -0.0194444      0.188889     -0.102778 
>fracprint(inv(A))
   53/360    -13/90    23/360 
  -11/180      1/45    19/180 
   -7/360     17/90   -37/360 
Matlab and Euler can be set to print in fractional form.
>fracformat; inv(A), longformat;
             53/360              -13/90              23/360 
            -11/180                1/45              19/180 
             -7/360               17/90             -37/360 
The following Frobenius norm is defined in Matlab. We would have to
define it in Euler.
>sqrt(max(abs(eigenvalues(A.A'))))
15
>real(eigenvalues(A))
[ 4.89897948557  -4.89897948557  15 ]
Euler and Matlab can compute a decomposition into singular values.
>{U,v,W}=svd(A); v,
[ 15  3.46410161514  6.92820323028 ]
The ::= defines a matrix for Maxima and Euler.
>A ::= magic(3)
                  8                   1                   6 
                  3                   5                   7 
                  4                   9                   2 
So we can symbolically compute the eigenvalues.
>&eigenvalues(A)[1]
                     [- 2 sqrt(6), 2 sqrt(6), 15]

>shortestformat; A=magic(4)
       16         2         3        13 
        5        11        10         8 
        9         7         6        12 
        4        14        15         1 
By chance, commuting does not change the sum of the diagonals.
>A=A[:,[1,3,2,4]]
       16         3         2        13 
        5        10        11         8 
        9         6         7        12 
        4        15        14         1 
>sum(A)', sum(A')', sum(diag(A,0)), sum(diag(flipy(A),0))
[ 34  34  34  34 ]
[ 34  34  34  34 ]
34
34
>det(A)
0
>rank(A)
3
We compute a table of ranks.
>r=(3:24)'|0
        3         0 
        4         0 
        5         0 
        6         0 
        7         0 
        8         0 
        9         0 
       10         0 
       11         0 
       12         0 
       13         0 
       14         0 
       15         0 
       16         0 
       17         0 
       18         0 
       19         0 
       20         0 
       21         0 
       22         0 
       23         0 
       24         0 
>for i=1 to rows(r); r[i,2]=rank(magic(r[i,1])); end;
>r
        3         3 
        4         3 
        5         5 
        6         5 
        7         7 
        8         3 
        9         9 
       10         7 
       11        11 
       12         3 
       13        13 
       14         9 
       15        15 
       16         3 
       17        17 
       18        11 
       19        19 
       20         3 
       21        21 
       22        13 
       23        23 
       24         3 
>plot2d(r'[1],r'[2],>bar,>insimg);

Molders Introduction to Matlab in Euler

The following plot shows the structure of the 9x9 magic square.
>n=9; A=magic(n); plot3d(A/totalmax(A)*n/2,>bar,<frame, ...
>  height=50°,angle=20°,>insimg);

Molders Introduction to Matlab in Euler

Cryptography

Another example of Matlab used for integer computation. I think this
is an abuse of the program. But I tried to do this nevertheless.

First, the Euler function char works only for numbers, not for vectors
of numbers. So we write such a function. 

This function should be done in C code. It would be very easy to do in
Euler using the TC Compiler.
>function x2s (x:real vector) ...
s="";
for i=1 to length(x)
  s=s+char(x[i]);
end;
return s;
endfunction
For a test, translate the ASCII codes 48 to 57 to letters.
>x2s(48:57)
0123456789
>s=x2s(32:96)
 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`
There is also a predefined function chartostr, which does just what we
need.
>chartostr(32:96)
 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`
This is the converse function written in Euler.
>function s2x (s:string) ...
  n=strlen(s);
  x=zeros(1,n);
  for i=1 to n;
    x[i]=ascii(s);
    s=substring(s,2,n-i+1);
  end;
  return x;
endfunction
>s2x(s)
[ 32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49
50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67
68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
86  87  88  89  90  91  92  93  94  95  96 ]
The builtin function strtochar does the same.
>strtochar(s)
[ 32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49
50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67
68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
86  87  88  89  90  91  92  93  94  95  96 ]
Now the crypto algorithms maps pairs of letters to pairs.

Molders Introduction to Matlab in Euler

>longformat; p=97
97
>x=[52;54]
                 52 
                 54 
>A=[71,2;2,26]
                 71                   2 
                  2                  26 
>y=mod(A.x,97)
                 17 
                 53 
>mod(A.y,97)
                 52 
                 54 
This function has its own inverse. The reason is the following.
>mod(A.A,97)
                  1                   0 
                  0                   1 
Here is a shorter version of the Matlab code.
>function crypto (s:string) ...
global A,p;
x=strtochar(s)-32;
X=redim(x,2,(length(x)+1)/2);
X=mod(A.X,p);
x=redim(X,1,length(x))+32;
return chartostr(x);
endfunction
>s=crypto("Hello World")
H-?36 WZU{s
>crypto(s)
Hello World
>crypto(crypto("Hello World!"))
Hello World!

The 3n+1 Sequence

Another integer problem. It is not known, if the sequence

Molders Introduction to Matlab in Euler

reaches 1 from all start values.
>function threenplus1 (n) ...
N=n;
repeat
  if mod(n,2)==0 then n=n/2; N=N|n;
  else n=3*n+1; N=N|n;
  endif
  until n==1;
end;
return N
endfunction
>N=threenplus1(7)
[ 7  22  11  34  17  52  26  13  40  20  10  5  16  8  4  2  1 ]
>plot2d(N); insimg;

Molders Introduction to Matlab in Euler

The values become large rather quickly. So it is best to use
logarithmic plots.
>plot2d(N,>logplot); insimg;

Molders Introduction to Matlab in Euler

Euler does not have a nice user interfaces for interactions. However,
it can receive keys and mouse actions from the user. The following
waits for key up or down, or for return to stop.
>function test ...
n=7;
repeat
  N=threenplus1(n);
  plot2d(N,>logplot,title="n="+ ...
     n+". Press up, down or return");
  k=key;
  if k==1 then n=n+1;
  elseif k==2 then n=max(n-1,1);
  elseif k==13 then return;
  endif;
end;
endfunction
>test;
>insimg;

Molders Introduction to Matlab in Euler

Floating Point Arithmetic

This is a section about floating point arithmetic. Here is the epsilon
Euler uses for internal checks. It can be changed.
>epsilon
2.22044604925e-012
It is 1000 times the smallest number such that 1+x is not x.
>printhex(epsilon)
2.7100000000000*16^-10
printhex prints the hexadecimal representation of the IEEE number in
Euler.
>printhex(2^-52)
1.0000000000000*16^-13
There is also a dual print.
>printdual(2^-52)
1.0000000000000000000000000000000000000000000000000000*2^-52
Note that 0.1 is not exactly representable.
>printhex(0.1)
1.999999999999A*16^-1
Thus adding 0.1 does not meet 1 exactly.
>x=0; loop 1 to 10; x=x+0.1; end; printdual(x)
1.1111111111111111111111111111111111111111111111111111*2^-1
However, the following works.
>0:0.1:1
[ 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 ]
As an example, the following expression does not evaluate to 0
exactly.
>longestformat; 1-(4/3-1)*3
 2.220446049250313e-016 
>&1-(4/3-1)*3
                                  0

In contrast to Matlab, Euler returns an error when solving the
following linear system. The reason is the internal epsilon.
>[17,5;1.7,0.5]\[22;2.2]
Determinant zero!

Error in :
[17,5;1.7,0.5]\[22;2.2]
                       ^
>fit([17,5;1.7,0.5],[22;2.2])
      1.294117647058824 
                      0 
>svdsolve([17,5;1.7,0.5],[22;2.2])
      1.191082802547771 
     0.3503184713375797 
An interesting example is the evaluation of the following polynomial.
First an exact plot.
>x = 0.988:.0001:1.012;
>plot2d(x,(x-1)^7,color=red,thickness=2);
>p &= expand((x-1)^7)
          7      6       5       4       3       2
         x  - 7 x  + 21 x  - 35 x  + 35 x  - 21 x  + 7 x - 1

If we evaluated the expanded polynomial, we get lots of errors. Note
that the y-scale is rather small.
>y = p(x);
>plot2d(x,y,color=blue,>add); insimg;

Molders Introduction to Matlab in Euler

Euler has xpolyval, which uses an exact residuum to iterate to the
correct evaluation.
>longformat; p=polycons(ones(1,7))
[ -1  7  -21  35  -35  21  -7  1 ]
>plot2d(x,xpolyval(p,x,eps=1e-16),color=red,thickness=2); insimg;

Molders Introduction to Matlab in Euler

Biorythm

This is one of the exercises. Bio rythm is of course nonsense. But it
makes a great programming exercise.
>function bio (birth, now) ...
  t=linspace(now-8,now+8,200);
  plot2d(t-now,100*sin(2*pi/33*(t-birth)),color=red,);
  plot2d(t-now,100*sin(2*pi/28*(t-birth)),color=green,>add);
  plot2d(t-now,100*sin(2*pi/23*(t-birth)),color=blue,>add);
  title("red=intellectual, green=emotional, blue=physical");
endfunction
Today I am not at the height of my mental power. But emotionally, I am
okay. So maybe this is the reason, why I solve this unrational
example.
>bio(day(1956,5,15),day(2012,14,7)); insimg;

Molders Introduction to Matlab in Euler

Examples Homepage