This is related to post in Math Less Travelled by Brent Yorgey (see also here by John D. Cook). They call the curve „mystery curve“. Of course, I had to try this in EMT (Euler Math Toolbox).
The first attempt was the following.
>function colsum (A) := sum(A')';
>function plotone (a:vector,f:vector,r=none,grid=0,thickness=1) ...
$ t=linspace(0,2*pi,5000);
$ z=exp(I*t)+colsum(a'*exp(I*f'*t));
$ if r==none then r=totalmax(abs(z))*1.05; endif;
$ plot2d(re(z),im(z),=r,=grid,=thickness);
$ endfunction
>fullwindow; plotone([1/2,I/3],[6,-14],thickness=2):
The function colsum() should be pre-defined. I will add it to the next version of EMT. Matlab sums columns by default, but EMT sums rows of matrices, which is more convenient in most cases. The function plotone() plots
\(\gamma(t) = e^{it} + \sum_k a_k e^{i f_k t}.\)
in the complex plane. In the code, t is a row vector and f and are used as column vectors, so that we have to add the rows of the matrix to get the sum. The function gamma(t) is in fact a sequence of circular movements. One popular toy that does the same for one element in the sum is the Spirograph. But the three terms of the mystery curve look a lot more interesting.
The curve is evaluated in 5000 points for simplicity. This is not a good way to evaluate such functions. Since the plot can get arbitrarily complex, it should be evaluated adaptively. This avoids unnecessarily many points and allows more complex curves by increasing the number of points. the function plot2d() used adpativeeval() for this. So you can simply do the plot like this:
>function map g (t;a,f) := exp(I*t)+sum(a*exp(I*f*t));
>function map gre (t;a,f) := re(g(t,a,f))
>function map gim (t;a,f) := im(g(t,a,f))
>plot2d({{"gre",[1/2,I/3],[6,-14]}},{{"gim",[1/2,I/3],[6,-14]}}, ...
> xmin=0,xmax=2pi,r=2):
This is indeed faster, because the plot of that many points takes most of the time. So let us update our plotone() function with this trick.
>function plotone (a:vector,f:vector,grid=0,thickness=1) ...
$plot2d({{"gre",a,f}},{{"gim",a,f}},xmin=0,xmax=2pi,
$ >square,=grid,=thickness);
$endfunction
We see, that we no longer can compute the value r in the code, since it is lost inside plot2d(). But we can use the >square parameter of plot2d() instead.
For most flexibility, we have to call adaptiveeval() ourselves.
>function plotone (a:vector,f:vector,r=none,grid=0,thickness=1,eps=0.01) ...
$ {z1,z2}=adaptiveeval({{"gre",a,f}},{{"gim",a,f}},0,2pi,=eps);
$ if r==none then r=totalmax(sqrt(z1^2+z2^2))*1.05; endif;
$ plot2d(z1,z2,=r,=grid,=thickness);
$endfunction
>plotone([1/2,I/3],[6,-14],thickness=2):
This allows to make 100 good look random samples. We use the window() command to split the screen into 10×10 plots. Note that eps=0.1 only because more is not necessary for such small plots.
>function makeone (n=3,m=10,eps=0.1) ...
$ a=intrandom(n,m)/m;
$ f=intrandom(n,m);
$ plotone(a,f,grid=0,=eps);
$ endfunction
>function makeall (n=10) ...
$ clg;
$ for i=0 to n-1;
$ for j=0 to n-1;
$ window([i,j,i+1,j+1]*1024/n);
$ hold on; makeone(); hold off;
$ end;
$ end;
$ shrinkwindow();
$ endfunction
>makeall():
