The code below creates animated phaseportraits for arbitrary, autonomous system of ODEs. Using the Lorenz equations as an example, I take the user through the creation of the animation so that he or she knows what they are really looking at when the code compiles and runs. There are two significant problems with the code that reduce the power of the program, but in the main it works.
Please see my statements of these problems, inserted in the code below. Note that code lines
begin with '>' and comment lines do not.
> restart: > with(plots): > with(DEtools): > 1+1: > sys:= [diff(x(t),t) = -s*x(t)+s*y(t), > diff(y(t),t) = r*x(t)-y(t)-x(t)*z(t), > diff(z(t),t) = -b*z(t)+x(t)*y(t)]; > s:=10:r:=28:b:=8/3: > inits:=[[x(0)=2,y(0)=3,z(0)=5]];
Programming Problem: Phaseportrait, the function in Maple that this worksheet will make the most use of, creates a two-dimensional phaseportrait for a given system of ODEs for a given set of initial conditions. You can give the program a list of initial conditions to create phaseportraits for and the program will draw the trajectories all on the same graph, but it is difficult and perhaps impossible to create a small procedure that generates a set of initial conditions for this function (maybe you can get it to work: see initialconditions or maybe DEtools. You might also try DEtools[phaseportrait]).
I would love to see a small procedure that could specify arbitrarily many initial conditions in an arbitrarily small neighborhood (or ball/cell)–it would greatly aid in the analysis of sensitive dependence on initial conditions for the system.
> Lxz:= proc(n); > phaseportrait(sys,[x(t),y(t),z(t)], > t=0..n, > inits, > stepsize=.005, > scene=[x(t),z(t)], > linecolor=t); > end: > Lxy:= proc(n); > phaseportrait(sys,[x(t),y(t),z(t)], > t=0..n, > inits, > stepsize=.005, > scene=[x(t),y(t)], > linecolor=t); > end: > lzlistxz:= [Lxz(1), Lxz(2), Lxz(3), Lxz(4), Lxz(5), Lxz(6), Lxz(7), Lxz(8), Lxz(9), Lxz(10), Lxz(11), Lxz(12), Lxz(13), Lxz(14), Lxz(15), Lxz(16), Lxz(17), Lxz(18), Lxz(19), Lxz(20), Lxz(21), Lxz(22), Lxz(23), Lxz(24), Lxz(25), Lxz(26), Lxz(27), Lxz(28), Lxz(29), Lxz(30)]: > lzlistxy:= [Lxy(1), Lxy(2), Lxy(3), Lxy(4), Lxy(5), Lxy(6), Lxy(7), Lxy(8), Lxy(9), Lxy(10), Lxy(11), Lxy(12), Lxy(13), Lxy(14), Lxy(15), Lxy(16), Lxy(17), Lxy(18), Lxy(19), Lxy(20), Lxy(21), Lxy(22), Lxy(23), Lxy(24), Lxy(25), Lxy(26), Lxy(27), Lxy(28), Lxy(29), Lxy(30)]: > display(lzlistxz, insequence=true); > display(lzlistxy, insequence=true);
Programming Problem: As hard as I have tried, I can’t get the display command to take a generated list (as opposed to a written list, such as lzlistxy above) of phaseportraits for input.
I’d like to be able to generate a list of phaseportraits to display based on a starting time t1 and an ending time t2 and a stepsize (e.g. if t1 = 1 and t2 =2 and the stepsize is 0.25, I’d like Maple to generate the phaseportraits Lxz(1), Lxz(1.25), Lxz(1.5), Lxz(1.75), and Lxz(2)).
The problem I have encountered is not getting the computer to generate such a list–the problem mainly seems to lie in the format of the elements of the list.
Instead of holding the elements in the list as actual plots, Maple holds them as instructions for plots. Unfortunately, the display command does not accept a list of instructions for plots as an input, so I’m at an impasse.
A small procedure that gets around this problem would greatly increase the usability of this program since it would give the user the ability to create arbitrarily smooth animation for an arbitrary length of time.
One solution, applicable elsewhere, is to have a procedure to generate an arbitrary linear sequence. Thus
`for` := proc(f,t,b) local x,sq; description `Compute a sequence of numbers from F to T by [B]|1.`; sq := NULL; for x from f to t by `if`(2 < nargs, b, 1) do sq := sq,x od; sq end: # `for`
You can then apply it to your problem in the following manner
lzlistxz := [seq(Lxz(xx), xx=`for`(1,2,0.2))]: display(lzlistxz, insequence=true);
I really can’t see what you find difficult here. Both problems seem quite straightforward to me.
First, here’s a procedure to randomly generate n initial conditions within epsilon (in each
coordinate) of [x0,y0,z0]:
makeinits:= proc(x0,y0,z0,n,epsilon) local i, rnd; rnd:= random(-100 .. 100); [ seq([x(0)=x0+rnd()*epsilon/100, y(0)=y0+rnd()*epsilon/100, z(0)=z0+rnd()*epsilon/100], i=1..n)]; end;
So to produce such a set of initial conditions, you can say e.g.
inits:= makeinits(3,4,5, 10, 0.1);
And to generate a list of phaseportraits, all you need is something like this:
lzlistxz:= [seq(Lxz(i), i=1..30)]: display(lzlistxz, insequence=true);
You can avoid having to type that long list by using a construction like this
lzlistxz:=eval( [ seq( 'Lxz'(k), k=1..30) ] );
It seems to me that a lot of unnecessary computation is going on, though. Here is a quick attempt at handling the problem differently. The idea is to compute only Lxz(30), and then make an animation out of that plot structure.
The procedure animatePLOT below is made to handle a 2-dimensional PLOT structure containing only one curve. It calls one of two subprocedures.
One of these (animatePLOT1) handles a PLOT structure containing one list of points (the typical output from plotting one curve), the other (animatePLOT2) handles a PLOT structure containing a sequence of lists of points, which is the output from phaseportrait, when linecolor has been specified.
The PLOT structure that can be input to animatePLOT2 must be of the form
PLOT(CURVES(...),..). The syntax is
animatePLOT( p, n );
where p is the PLOT structure and the optional n is the number of frames (default=number of points in the list (or members of the sequence)).
> restart; > animatePLOT:=proc(p::PLOT) > if nops(select(has,[op(p)],CURVES))>1 then ERROR("More than one curve") fi; > if nops(select(type,{op(op(1,p))},list))>1 then animatePLOT2(args) else animatePLOT1(args) fi > end: > > animatePLOT2:=proc(p) local n,p1,k,L,N,m,S,S1,S2,LC; > L:=select(type,[op(op(1,p))],listlist); > S:=remove(type,[op(op(1,p))],listlist); > S1:=remove(has,S,{COLOR,COLOUR}); > S2:=op(select(has,S,{COLOR,COLOUR})); > LC:=[op(S2)]; > N:=nops(L); > if nargs>1 then n:=args[-1] else n:=N fi; > m:=floor(N/n); > for k from 1 to n do > p1[k]:=subsop(1=CURVES(op(L[1..k*m]),COLOR(op(LC[1..k*m+1])),op(S1)),p); > od; > plots[display]([seq(p1[k],k=1..n)],insequence=true) > end: > > animatePLOT1:=proc(p) local n,p1,k,L,N,m; > L:=op([1,1],p); > N:=nops(L); > if nargs>1 then n:=args[-1] else n:=N fi; > m:=floor(N/n); > for k from 1 to n do > p1[k]:=subsop([1,1]=L[1..k*m],p); > od; > plots[display]([seq(p1[k],k=1..n)],insequence=true) > end: > p:=plot(sin,0..2*Pi):animatePLOT(p); > p:=plot([sin,cos,sin*cos,sin+cos],0..2*Pi): > animatePLOT(p); > p:=plot(sin,0..2*Pi,style=point): > animatePLOT(p); > with(DEtools): > sys:=[diff(x(t),t)=-s*x(t)+s*y(t),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t), diff(z(t),t)=-b*z(t)+x(t)*y(t)]; > s:=10:r:=28:b:=8/3: > inits:=[[x(0)=2,y(0)=3,z(0)=5]]; > p:=phaseportrait(sys,[x(t),y(t),z(t)],t=0..30,inits,stepsize=.005, scene=[x(t),z(t)],linecolor=t): > animatePLOT(p,30);
Needless to say, this procedure is rather raw and needs improvement. I’m aware of animatecurve, but that seems to use the same number of points for each frame, which doesn’t produce good pictures at the end unless numpoints is set to a high value.