The main transformation made by combine/radical is
  ... * x^(m/d)*y^(n/d) * ...  ==>  ... * (x^m*y^n)^(1/d) * ...         
where x and y are both positive, (i.e. signum(x)=signum(y)=1) and 0<m, n<d where m, n, d are integers.
For example: 2^(1/2)*3^(1/2) ==> 6^(1/2)

combine(2^(1/2)*3^(1/2),radical);

                                  1/2
                                 6

combine(2^(1/8)*3^(1/8)*35^(1/8),radical);

                                  (1/8)
                               210

combine(2^(7/8)*3^(1/8)*35^(1/8),radical);

                                   (1/8)
                              13440

combine(2^(7/8)*3^(1/4)*35^(1/8),radical);

                            (1/4)     (1/8)
                           3      4480

Doesn't the above satisfy the conditions stated in the Help? e.g. n<d

combine(x^(m/d)*y^(n/d),radical) assuming x::positive,y::positive,m::integer,n::integer,d::integer,n<d;

                             (m/d)  (n/d)
                            x      y

Shouldn't the last line be ?
                               m  n (1/d)
                             (x  y )

combine(x^(m/d)*y^(n/d),radical) assuming n::posint,m::posint,d::posint, n=m;

                             (m/d)  (n/d)
                            x      y

maybe I'm not using asssuming properly here, am I?

prompted by Robert's answer I searched mapleprimes and stumbled upon this very nice blog entry by Stephen Forrest, which discusses prime numbers and all that, you may find it useful and dig for ideas there:

http://www.mapleprimes.com/blog/stephen-forrest/finding-the-set-of-prime-numbers-le-n

prompted by Robert's answer I searched mapleprimes and stumbled upon this very nice blog entry by Stephen Forrest, which discusses prime numbers and all that, you may find it useful and dig for ideas there:

http://www.mapleprimes.com/blog/stephen-forrest/finding-the-set-of-prime-numbers-le-n

The last plot produced by the getarrows procedure:

you're welcome to hijack my blog any time

;-)

Using the latest version of the getarrows procedure, if I run the code with _Env_in_maplet set to true, it takes 8 seconds to execute, while if I comment it out and add interface(warnlevel=0): or interface('warnlevel'=0): outside the loop, the execution time jumps to 30 seconds. That's quite a difference: _Env_in_maplet must be advertized at once!

Placing a try/finally inside the loop does not change execution time in a noticeable way.

Placing the dsolve outside of the loop instead of inside (as you suggested below)  does not change execution time in a noticeable way, but it is cleaner and I certainly like it this way.

The isochrones with a few arrows, together with the stable and unstable manifold, in color:

Edit. The problem below with the tt came about because of the do loop, which turns out to be useless, as pointed out by Joe, so I have moved this comment to the bottom. In other words, ignore the paragraphs below this line!

___________________________________________________

The following fix was useful:

L2 := seq( select( s->s[1]=tt, L ), tt=0..T ):
for tt from 0 to 10 do  
  f := tt-> {op(L2[tt+1,..][..,2..3])}:
end do:

Either my use of t in f := t-> or in s->s[1]=t (or in both) was messing things up. I just need to replace t by, say, tt and it's okay now.

thanks Joe, I hadn't noticed that the loop wasn't doing anything useful! It's gone now.

Thanks Joe, it appears that I missed this message of yours. I will incorporate your suggestions into the next incarnation of the code. Could take a few days, but I'll get there. Thanks again Joe,

Patrick.

All I had to do was to add the option connect=true:

p := seq( plots:-pointplot( f(t), connect=true, thickness=3 ),t=0..T ):
plots:-display(p);

One thing I haven't had time to look into, my next trivial assignment, is to extrapolate between the points and draw smooth lines, as opposed to the current pointplot.

The isochrones at time t=1 is a bit ugly, but that's life.

I don't think I will attempt to add "t=1", "t=2", "t=3" labels. I'll probably just add it later with pstricks, if I ever need to "publish" a figure with isochrones.

Today I want to blog about how I managed to connect the points created by Robert Israel's getarrows procedure and thus create a plot of isochrone lines. The getarrows procedure was originally intended to create a list of points at which to draw the arrows of a vector field. This procedure I'm going to call a "salvo" because the desired effect is to simulate a salvo of arrows thrown from a neighborhood of a given point (which in my case of interest is chosen on the stable manifold of a saddle-point system) and to track the movement of the salvo at fixed time intervals such as t=0,1,2, 3, etc. One arrow being thrown from the stable manifold "exactly" will converge to the saddle-point, while the other arrows being thrown just off of it will diverge more or less fast. The advantage of the salvo device, as opposed to the usual vector field on a fixed grid, is that it shows graphically the speed of evolution of the system. A minor disadvantage, I suppose you can call it that, is that a point of origin must be selected and all isochrones thereof depend on that point of origin. All lines connecting the ends of arrows thrown from the same neighborhood and sampled at the same time are called isochrone lines, since iso=same and chronos=time (and no I didn't major in Greek so please correct me if I err in this neologism). If you understand this explanation, you're a most excellent reader and I do urge you to read on.

Since Robert and Joe and pagan and others had done all the hard work, adapting the procedure turned out to be easy. It took me a mere 5 minutes while playing with my son. Yeah right, you can say that again. Okay. It took me a mere 5 hours while being frequently disturbed by my baby son.

The first thing I needed was to output time as well as the point, i.e. [t, x(t), y(t)]. That was, I mean it this time, trivial.

The second thing I needed was to create a list of the iso-points, that is to separate the [t1,x(t1),y(t1)] from the [t2,x(t2),y(t2)] and to organize them in such a way as to be usable with pointplot. That took forever. I had to learn to use "select" and, most difficult, how to apply "seq" to "select". I also had to learn how to extract lists from lists with calls like this L[1,..][..,2..3]. That took a while too. Then I needed to put the mess together. Now please be warned: I didn't know how to avoid a loop so my code contains the much abhorred for/do structure. However, I'm pleased to report that it is only 4 lines and does seem to be quite fast. Eureka! (and here the other Greek word I know).

I know there's a way to write all this more succinctly and more efficiently. Feel free to make suggestions. But don't feel urged, for the wonderful help I have so far received has made my code good enough for my current purpose.

restart;
# Using the undocumented _Env_in_maplet command
getarrows0 := proc(x0, y0, xd, yd, T, x1, y1)
local DEsol, tm,j;
global x,y,t;
    _Env_in_maplet := true; # silence event warnings during integration
    DEsol:= dsolve({xd,yd,x(0)=x0,y(0)=y0}
                   , numeric
                   , 'events' = [[y(t)-y1,'halt'],[x(t)-x1,'halt'],[x(t),'halt']]
                   , 'output=procedurelist'
                  );
    tm:= floor(subs(DEsol(T),t));
    seq(subs(DEsol(j), [t, x(t), y(t)]), j = 0 .. tm);
end proc:

xdot := diff(x(t),t) = x(t)-y(t):
ydot := diff(y(t),t) = y(t)*(1-x(t)):

T := 10:
xmax := 3:
ymax := 3:
DEsol:= dsolve({xdot,ydot,x(0)=1-1e-10,y(0)=1-1e-10},numeric,output=listprocedure):
t1 := fsolve(rhs(DEsol(t)[2])=0.05, t=-100..0):
x0 := rhs(DEsol(t1)[2]):
y0 := rhs(DEsol(t1)[3]):

y0 := [ seq( seq( y0*(1+i/10*(1e-1)^j),i=-90..99 ), j=1..10) ]:

st := time[real]():
L := map2(getarrows0, x0, y0, xdot, ydot, T, xmax, ymax): # using _Env_in_maplet
time[real]() - st;

                                11.076

L2 := seq( select( s->s[1]=t, L ), t=0..T ):

for t from 0 to 10 do  
  f := t-> {op(L2[t+1,..][..,2..3])}:
end do:
 
with(plots):
p := seq( pointplot( f(t) ),t=0..T ):
display(p);

Thanks Joe, I had totally missed pagan's hint.

I would instinctively avoid an undocumented variable (it's hard enough for me to keep up with the documented ones), but I noticed a significant difference in speed of execution, with

st := time[real]():
init1 := map2(getarrows0, x0, y0, xdot, ydot, T, xmax, ymax): # using _Env_in_maplet

time[real]() - st;

                                3.246


st := time[real]():
init1 := map2(getarrows, x0, y0, xdot, ydot, T, xmax, ymax): # using try/finally
time[real]() - st;

                                27.114

for one particular grid of points (the last one I posted about). I will need to create a fine grid (I think) and I had noticed that the time of execution was getting slow as I increased the number of points in the grid and/or raised digits. So this efficient method is very good to have.

Thanks a lot for showing these great tools, Joe; after all  _Env_in_maplet will come in very handy.

Christmas gift-giving, that is:

http://graphics8.nytimes.com/images/blogs/freakonomics/pdf/WaldfogelDeadweightLossXmas.pdf

The dismal science at its best

;-)

came across this fairly recent post while searching "seq", and remembered your fga (frequently given answers) project on sum v. add, thought I'd mark it here for the record:

http://www.mapleprimes.com/forum/sumoverseqsoutput#comment-26085