Often, when studying iterated functions, one needs to find and classify the fixed points. Graphical methods go one an approximate idea of where the fixed points are but we usually need precise values. In Maple, this is very easy. If we have a function F, all we need to do to find the fixed points is solve(F(x)=x,x);. From there, you can easily plug the values you find into the deivative of F in order to classify all the fixed points.
Maple also allows you to easily find the points on a cycle as well. For example, if you want to find points on a 2-cycle for an iterated function F, you can do solve((F@@2)(x)=x,x);. (The symbol (F@@2)(x) means F iterated twice.) Remember that this also gives you the fixed points as well so to find the 2-cycle, you must take the points which are solutions of this last equation but not of F(x)=x.
with(plots):
graphAnalysis:=proc(func,range,pt,n)
local A1,A2,i,pic,p,B,C;
p:=pt;
A1:=plot(func(x),x=range,color=black, scaling=unconstrained):
A2:=plot(x,x=range,color=black,scaling=unconstrained):
pic:={A1,A2}:
for i from 0 to n do
B:=PLOT(CURVES([[p,p],[p,func(p)]],COLOR(HUE,0))):
C:=PLOT(CURVES([[p,func(p)],[func(p),func(p)]], COLOR(HUE,0.7))):
p:=func(p):
pic:=pic union {B,C}:
od:
display(pic);
end:
f:= x -> x^3-2.5*x;
graphAnalysis(f,-2..2,0.8,60);
| After pressing Enter following the call to the graphAnalysis procedure, Maple will graph the line y=x as well as the graph of the function y=f(x). Then, it will start at the point (pt,pt) and perform (in this case) 60 steps of the graphical analysis process. Maple displays the output all at once whereas some packages use animation displaying one step at a time. | |
with(plot):
kochf:=(A,B)->[A,[(2/3)*A[1]+(1/3)*B[1],(2/3)*A[2]+(1/3)*B[2]], [(1/2)*(B[1]+A[1])+(1/3)*(A[2]-B[2]),(1/2)*(B[2]+A[2])+(1/3)*(B[1]-A[1])], [(1/3)*A[1]+(2/3)*B[1],(1/3)*A[2]+(2/3)*B[2]],B];
iterate:=proc(n,f)
local i,points,C;
points:=[[0,0],[1,0]];
for i from 1 to n do
C:=zip((x,y)->[x,y], points[1..-2], points[2..-1]);
points:=evalf(map(x->op(f(x[1],x[2])[1..-2]),C));
points:=[op(points),[1,0]];
od;
PLOT(CURVES(points),AXESSTYLE(NONE));
end:
iterate(5,kochf);
sqrf:=(A,B)->[A,[(2/3)*A[1]+(1/3)*B[1],(2/3)*A[2]+(1/3)*B[2]], [A[1]+(1/3)*(B[1]-A[1])+(1/3)*(A[2]-B[2]),A[2]+(1/3)*(B[2]-A[2])+(1/3)*(B[1]-A[1])], [A[1]+(2/3)*(B[1]-A[1])+(1/3)*(A[2]-B[2]),A[2]+(2/3)*(B[2]-A[2])+(1/3)*(B[1]-A[1])],[(1/3)*A[1]+(2/3)*B[1],(1/3)*A[2]+(2/3)*B[2]],B];
iterate(5,sqrf);
The key to this Maple file is the procedure iterate. It takes as input an integer n and a function f. The function f should be a two variable function and it is assumed that each of the two variable is a point. The function must map the two points into a list of points, each defined uniquely in terms of the two input points. Then, the procedure iterate starts with the unit interval, in this case the list of points [[0,0],[1,0]].
Then iterate repeats the following n times: it takes a list of points and applies f to any two adjacent points in the list. Then iterate plots the list of points we obtain after that loop and draws a line between adajcent points.
Note that in the above example we have defined two different functions on points which satisfy the proper criteria. The koch function produces the standard Koch snowflake while the function sqrf adds a square to the interval instead of an equilateral triangle. Below, we've shown the output of this iteration with the koch function.
with(plots):
f:= (x,i) -> i*x*(1-x);
p0:=0.25;
points:={}:
for i from 2.4 by 0.02 to 4 do
p:=p0:
for j from 0 to 100 do
p:=f(p,i):
if j>40 then points:=points union {[i,p]} fi:
od:
od:
PLOT(POINTS(op(points),SYMBOL(POINT)));
Pressing Enter after the last line will cause Maple to display the orbit diagram as constructed by the given algorithm.
with(linalg):
with(plot):
iterfunc := proc(n,reps)
local i,p,j,points,die;
die:=rand(1..n):
points:=[]:
p:=vector([0,0]):
for i from 1 to reps do
j:=die();
p:=evalm(((A.j)&*(p-(p.j)))+(p.j));
if (i>50) then points:=[op(points),[p[1],p[2]]] fi;
od;
PLOT(POINTS(op(points),SYMBOL(POINT)),AXESSTYLE(NONE));
end:
A1:=array([[0.5,0],[0,0.5]]);
A2:=A1:
A3:=A1:
p1:=vector([0,0]);
p2:=vector([0,1]);
p3:=vector([1,0]);
iterfunc(3,5000);
iterfunc is a procedure taking two variables. This procedure assumes that global variables with names A1, A2, A3, ... and p1, p2, p3, ... exist. The Ai should by 2x2 arrays and the pi's should be points in the plane, declared prior to any use of the procedure iterfunc.
Here, iterfunc is told to plot the first 5000 points in this iterated function system which uses fi(X)=Ai(X)+pi as the affine maps.
Since Maple struggles when one tries to display over 5000 points, here's a trick to avoid bogging down your computer:
A1:=iterfunc(3,5000):
A2:=iterfunc(3,5000):
A3:=iterfunc(3,5000):
display(A1,A2,A3);
This will display 15000 points of the fractal. Since A1, A2, and A3 are created randomly. each one will be different. This process is slow but accomplishes the task - eventually.
(Remember that the display command only works when we include the plots package.)
restart:
with(plots):
f1:=p->[0.5,0.27*p[2]];
f2:=p->[-0.14*p[1]+0.26*p[2]+0.57,0.25*p[1]+0.22*p[2]-0.03];
f3:=p->[0.17*p[1]-0.22*p[2]+0.4,0.22*p[1]+0.176*p[2]+0.09];
f4:=p->[0.78*p[1]+0.07*p[2]+0.11,-0.07*p[1]+0.74*p[2]+0.27];
probf:=proc()
local z,n;
z:=rand()/1e10;
if z<=2 then n:=1
elif z<=17 then n:=2
elif z<=30 then n:=3
else n:=4
fi;
n;
end:
iterfunc:=proc(imax)
local pts,p,n,i;
pts:=NULL; p:=[0,0];
for i from 1 to imax do
n:=probf();
p:=(f.n)(p);
if (i>50) then pts:=pts, p
fi;
od;
plot([pts],axes=NONE, style=point, symbol=POINT,color=COLOR(RGB,0,1,0));
end:
planner:=proc(num)
local i,AA,pts,pics;
pics:=NULL;
for i from 1 to num do
pts:=[(f.i)([0,0]), (f.i)([1,0]), (f.i)([1,1]), (f.i)([0,1]), (f.i)([0,0])];
pics:=pics, PLOT(CURVES(pts,COLOR(RGB,1,0,0)), AXESSTYLE(NONE));
od;
pics:= pics, PLOT(CURVES([[0,0],[1,0],[1,1],[0,1],[0,0]]), AXESSTYLE(NONE));
display({pics});
end:
AA:=planner(4):
BB:=iterfunc(2000):
display({AA,BB});
The picture to the right displays the result of the last command. In the next few paragraphs, I will step through the parts of this Maple script.
First, the command restart: simply clears Maple's active memory. All lines need to be re-executed.
Next, I declare the four functions I am going to use in this example of IFS. The input of each function is assumed to be a sequence of length 2, i.e. a point in the plane. Note that compared to the script under the previous heading (Iterated Function Systems), we have done away with the use of matrices. It saves a little time in Maple to avoid matric multiplications.
Then the next piece of data you need is the probability function probf. The line z:=rand()/1e10; gives us a random integer between 0 and 100. As written above, probf output the integer 1 2% of the time, the integer 2 15% of the time, the integer 3 13% of the time and the number 4 70% of the time. To create new IFS, you would need to modify, add or remove the functions and redefine probf as desired.
Now iterfunc does pretty much what it did in the IFS script above this one but in a slightly different manner. The procedure planner takes as an argument the number of functiosn you are working with and will output a plot aobject. The plot object shows only a few things: 1) the unit square in black; and, for each function, 2) the parallelogram to which the corresponding functions maps the unit square. The output of planner helps you to see how each function deforms the unit square. Improvements could be made to this code because it is sometime difficult to see by how much the unit square gets rotated and whether it gets flipped (reflection through a line). One could display either the fractal put out by iterfunc or planner independently, as the case may be. However, showing them together helps one to see how the whole fractal (all that appears in an around the unit square) is affine self-similar to various parts of itself (that which appears in or near each of the transformed unit squares). This aide can also help you predict the look of a fractal based on the functions.
Note that in our example of the fractal fern, one of the functions not only shrinks and rotates but also flips the unit square. That's why the leaves on the right of the whole fern point upwards. Also, we created the stem of the fern with the function f1. You may notice that f1 collapses the whole unit square onto a vertical segment of length 0.27 with x coordinates equal to 0.5.
To produce a picture with more points, use the trick described in the previous subsection.
Behind our estimates and tricks lurks a theorem called the Collage Theorem that relates the overall structure of the fractal to the "transformed" unit squares, unit square after the fi are applied to it.