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 Google 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

Moler - Introduction to Matlab in EMT

>p=[-1,-1,1]
[-1,  -1,  1]

The function of EMT for polynomial roots is polysolve(p).

>r=polysolve(p)
[ -0.618034+0i ,  1.61803+0i  ]

Of course, this can be solved exactly with Maxima. In Matlab with the symbolic toolbox, this would be a solution by Maple.

>&solve(1/x=x-1)
                       1 - sqrt(5)      sqrt(5) + 1
                  [x = -----------, x = -----------]
                            2                2

Maxima prints nicer than Matlab. However, it a lot better to use Latex. EMT can display Latex formulas in comments or in the command line. In comments, we enter

  maxima: solve(1/x=x-1)

and get

Moler - Introduction to Matlab in EMT

The polysolve() function returns complex roots. So let us make the second root real.

>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.

Functions like the following one line function can be done in Matlab too (a bit more ugly, however).

>function f(x) := 1/x-(x-1)
>plot2d("f",0,4);
>solve("f",1)
1.61803398875
>plot2d(phi,0,>points,>add):

Moler - Introduction to Matlab in EMT

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,<grid);
 plot2d(u,v,style="--",>add);
 text(latex("\phi"),toscreen([phi/2,1.1]));
 text(latex("\phi-1"),toscreen([(1+phi)/2,-.05]));
 ctext("1",toscreen([-.05,.5]));
 ctext("1",toscreen([.5,-.05]));
 endfunction
>plotgr; insimg(crop=0.7);

Moler - Introduction to Matlab in EMT

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

Moler - Introduction to Matlab in EMT

A continued Fraction

The next example in Moler's book 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 (or the "fraction" operator), 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), fraction 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.188 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;
Real 39 x 1 matrix

     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 
    ...

Indeed there is the famous formula for the Fibonacci numbers.

Moler - Introduction to Matlab in EMT

>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 translation 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);

Moler - Introduction to Matlab in EMT

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 matrix norm (belonging to the Euclidean norm) is defined in Matlab. We would have to define it in Euler.

>sqrt(max(abs(eigenvalues(A.A'))))
15
>real(eigenvalues(A))
[15,  4.89897948557,  -4.89897948557]

Euler and Matlab can compute a decomposition into singular values.

>{U,v,W}=svd(A); v,
[15,  3.46410161514,  6.92820323028]

Since orthogonal multiplications do not change the Euclidean norm of a matrix, we can define the norm in the following way.

>function euclnorm (A) ...
 {U,v,W}=svd(A);
 return max(v);
 endfunction

Check.

>euclnorm(A)
15

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;
>for i=1 to rows(r); r[i,2]=rank(magic(r[i,1])); end;
>showlarge r
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);

Moler - Introduction to Matlab in EMT

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

Moler - Introduction to Matlab in EMT

Cryptography

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

The function char(x) works only for a single character. But chartostr(v) converts a vector of ASCII numbers.

>chartostr(48:57)
0123456789
>s=chartostr(32:96)
 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`

The built-in function strtochar() is the converse.

>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.

Moler - Introduction to Matlab in EMT

>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

Moler - Introduction to Matlab in EMT

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):

Moler - Introduction to Matlab in EMT

The values become large rather quickly. So it is best to use logarithmic plots.

>plot2d(N,>logplot):

Moler - Introduction to Matlab in EMT

Euler has an interfaces for interactions with the user.

>function test (n) ...
 N=threenplus1(n);
 plot2d(N,>logplot,title="3n+1 sequence");
 endfunction
>dragvalues("test","n",27,[7,107],100,0.5,0.8):

Moler - Introduction to Matlab in EMT

Floating Point Arithmetic

This is a section about floating point arithmetic. Here is the epsilon Euler uses for internal checks. It can be changed.

>longest epsilon
 2.220446049250313e-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

As a symbolic expression, it does.

>&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] ...
                       ^

The determinant if exactly 0. But we can find a best fit for solution.

>fit([17,5;1.7,0.5],[22;2.2])
      1.294117647058824 
                      0 

The following finds the fit among the fits that has minimal norm.

>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):

Moler - Introduction to Matlab in EMT

>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):

Moler - Introduction to Matlab in EMT

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):

Moler - Introduction to Matlab in EMT

Biorhythm

This is one of the exercises. Bio rhythm 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);
   labelbox(["intellectual","emotional","physical"],
     colors=[red,green,blue],x=0.5,w=0.4);
 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 tackle this irrational example.

>bio(day(1956,5,15),day(2012,14,7)):

Moler - Introduction to Matlab in EMT

The Ulam Spiral

The Ulam spiral places the prime numbers in a spiral around the center of a rectangular grid. This is a demo from Molder's introduction.

We need to find the columns and rows of the n-th element in the spiral. That is essentially a programming exercise.

>function map col (n) ...
  if n==1 then return 0; endif;
  k=ceil(sqrt(n)); // on bounary of k*k square
  if mod(k,2)==0 then k=k+1; endif; // k must be odd
  d=n-(k-2)^2; // d-th element there
  u=(k^2-(k-2)^2)/4; // 4*u elements on this boundary
  k=(k-1)/2; // distance from center of the boundary
  if d<=u then return k;
  elseif d<=2*u then return k-(d-u);
  elseif d<=3*u then return -k;
  else return -k+(d-3*u)
  endif;
  endfunction
>function map row (n) ...
  if n==1 then return 0; endif;
  k=ceil(sqrt(n)); // on bounary of k*k square
  if mod(k,2)==0 then k=k+1; endif; // k must be odd
  d=n-(k-2)^2; // d-th element there
  u=(k^2-(k-2)^2)/4; // 4*u elements on this boundary
  k=(k-1)/2; // distance from center of the boundary
  if d<=u then return k-d;
  elseif d<=2*u then return -k;
  elseif d<=3*u then return -k+(d-2*u);
  else return k
  endif;
  endfunction

Once we have these functions, we need to test.

>function testspiral (n) ...
  k=2*n+1;
  A=zeros(k,k);
  c=col(1:k^2); r=row(1:k^2);
  for i=1 to k^2; 
     A[r[i]+n+1,c[i]+n+1]=i;
  end;
  return A
  endfunction
>shortest testspiral(3)
       37        36        35        34        33        32        31 
       38        17        16        15        14        13        30 
       39        18         5         4         3        12        29 
       40        19         6         1         2        11        28 
       41        20         7         8         9        10        27 
       42        21        22        23        24        25        26 
       43        44        45        46        47        48        49 

Seems to work.

>function makespiral (n) ...
  k=2*n+1;
  A=zeros(k,k);
  c=col(1:k^2); r=row(1:k^2);
  p=primes(k^2);
  for i=p;
     A[r[i]+n+1,c[i]+n+1]=1;
  end;
  return A
  endfunction

We test this. Note that we find the primes up to 251001 and generate a matrix of size 501x501. The matrix is then used like an image.

>A=makespiral(250);
>insrgb(rgb(!A,!A,!A));

Moler - Introduction to Matlab in EMT

Examples