Advanced Search
All Posts : Education
I wrote a procedure for drawing a 1d phase portrait of an equation 
                        `x′`=f(x)
in Classic Maple, 
phase:=proc(f,r)
local a,b,bn,bs,U;
a:=op([1,1],plot(args));
bn:=-`-`(op([2,1..2],r));
bs:=[a[1,1],seq(`if`(
(a[i-1,2]-0.0001*bn*signum(a[i-1,2]))*(a[i,2]-0.0001*bn*signum(a[i,2]))<0,
`if`(a[i-1,2]*a[i,2]<0, (a[i-1,1]*a[i,2]-a[i-1,2]*a[i,1])/(a[i,2]-a[i-1,2]),
fsolve(f,subsop([2,1]=a[i,2]-0.02*bn,[2,2]=a[i,2]+0.02*bn,r))),NULL),
i=2..nops(a)),a[-1,1]];
b:=bs/bn;
U:=remove(x->abs(x[2,1])<0.04,
[seq([[b[i]+0.01,0],[b[i+1]-b[i]-0.02,0]],i=`if`(a[1,2]<0,2,1)..nops(b)-1,2),
seq([[b[i+1]-0.01,0],[b[i]-b[i+1]+0.02,0]],i=`if`(a[1,2]<0,1,2)..nops(b)-1,2)]);
plots[display](plots[arrow](
U,width=0.02,head_length=0.04,view=[b[1]..b[-1],-1..1],color=blue),
plots[pointplot](map(x->[x,0],b[2..-2]),symbol=circle,symbolsize=19),
plots[textplot](map(x->[x/bn,-0.07,sprintf("%.2f",x)],bs[2..-2])),
axes=none)
end:
It works as in the following example for f(x)=sin(x), 
phase(sin(x),x=-8..8);
For discontinuous functions coordinates of some singularities may be printed not very precisely. That can be fixed by increasing numpoints. For example, 
phase(1/sin(x),x=-8..8,numpoints=200);
produces the same picture as above. For very large intervals some of singularities may be missed, that can be also fixed by increasing numpoints, as in the following example, 
phase(-x^2,x=-1000..1000,numpoints=3000);
For Standard Maple some numbers have to be adjusted, in particular symbolsize can be set to 10 instead of 19, and something like head_width=0.05 can be added to the arrow command. Still, it doesn't produce the same picture in Standard Maple, because arrows in Standard Maple are not centered vertically - their body is located higher than it should. Finally, I would like to thank Gerald A. Edgar for inspiring me to extend the original version of the procedure that I wrote for continuous functions, to piecewise continuous functions. __________ Alec Mihailovs http://mihailovs.com/Alec/
phaseportrait plotting

Please Wait...