Carl Love

## 27187 Reputation

11 years, 341 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Basins of attraction...

@C_R The Wikipedia article "Attractor" contains a definition of basin of attraction. It's the mathematical generalization of the geological / geographical / hydrological concept of drainage basin or watershed.

## Extra frames...

@C_R I haven't run this yet, but in order to get such a slow and smooth animation to post, the OP very likely included a huge number of extra frames (10 times, or more) in the animation. These extra frames are not needed when viewing in a worksheet because you can slow down the animation with the toolbar controls. Not computing the extra frames should speed up the computation significantly.

## New worksheet...

@dharr As far as I can tell, the worksheet poly.mw attached to your Answer is the same as when it was originally attached. If this is not true, please include a comment in the worksheet to indicate the new part.

## diagnostic...

@Andiguys This is to diagnose the problem, not to fix it: Please insert this command after the double loop:

lprint(indices(s, pairs, indexorder));

Then execute the entire worksheet. Then post the entire executed worksheet.

## s[1] isn't an equation...

@Andiguys s[1] isn't an equation, so it doesn't have an rhs. It's just a number. So, remove the first rhs.

## -2 <= x <= 2...

@dharr p = sinnx (identically, not approximately) only for in -2..2, because of the real domain of arcsin(x/2). I don't know why the OP is trying to extend the domain to -2.1..2.1.

simplify ignores the domain restriction, just as it does with sin(2*arcsin(x)). Indeed, sin(arcsin(x)) automatically simplifies to x.

## The 3D plot...

@Andiguys Here is one way to get that 3D plot. Do this after the double loop. I needed to remove the higher values of Cv because Maximize rejected them.

CV:= [50, 55, 60, 65, 70]:
T0:= [0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7]:
S:= Matrix((nops(CV),nops(T0)), (i,j)-> s[CV[i],T0[j]][1]):
A:= (Cv,t0)-> CurveFitting:-ArrayInterpolation(
[CV, T0], S, [[Cv], [t0]], method= spline
)[1,1]:
plot3d(A, (min..max)(CV), (min..max)(T0), labels= ['Cv', '`&tau;0`', 'TRC']);

This isn't the most efficient way to use ArrayInterpolation, but the syntax is already daunting, and it's sufficiently efficient for this purpose.

## do, od, end, end do...

@Andiguys I didn't mean to put do at the end of the whole loop; I meant to put it and the end of the the line with the for:

for i from 1 to 10 do
i^2
od

od can be replaced by end or end do (they mean exactly the same thing), but not by do. The do on the first line cannot be replaced by anything.

## Inner seq converts vectors to list...

@Scot Gould Your code is sufficient to produce a list of vectors (which are actually aliased arrays (see ?ArrayTools:-Alias)), and that's fine for working within a worksheet. I used the inner [seq](p[]) to convert those vectors to lists so that I could copy-and-paste the output line to MaplePrimes as plaintext. It's only for that reason that I used the extra seq.

Copy-and-pasting vectors as plaintext (using good-old Ctrl-C Ctrl-V) produces this:

{Vector[row](4, [1, 2, 3, 4]), Vector[row](4, [1, 2, 4, 3]), Vector[row](4, [1, 4, 2, 3]), Vector[row](4, [2, 1, 3, 4]), Vector[row](4, [2, 1, 4, 3]), Vector[row](4, [2, 3, 1, 4]), Vector[row](4, [2, 3, 4, 1]), Vector[row](4, [2, 4, 1, 3]), Vector[row](4, [2, 4, 3, 1]), Vector[row](4, [4, 1, 2, 3]), Vector[row](4, [4, 2, 1, 3]), Vector[row](4, [4, 2, 3, 1])}

(for which I edited the color in MaplePrimes).

Copy-and-pasting vectors as images (using Copy as Image from the context menu) produces this:

MaplePrimes always doubles copy-and-pasted images of the output portions of execution groups. The extra can be removed by a single backspace. But no other editing can be done, not even changing the color.

Copy-and-pasting an execution group containing both a command and its output with Ctrl-C Ctrl-V and then editing the colors and font weight and removing extra blank lines (by replacing Return with Shift-Return) produces this:

{seq}(p[], p= Iterator:-TopologicalSorts(4, {2<3}));
{[1, 2, 3, 4], [1, 2, 4, 3], [1, 4, 2, 3], [4, 1, 2, 3],
[2, 1, 3, 4], [2, 1, 4, 3], [2, 4, 1, 3], [4, 2, 1, 3],
[2, 3, 1, 4], [2, 3, 4, 1], [2, 4, 3, 1], [4, 2, 3, 1]}

Doing the same as above with the only change being using Copy as Image produces a result almost identical to above. The only difference is that there's no Return / Shift-Return distinction, so the extra blank lines can be removed with simple backspace.

My browser is Google Chrome, which doesn't allow its own context menu for the Paste, telling you to use Ctrl-V instead. Perhaps another browser would handle it differently.

## Confusing name...

@Kitonum I think that the command is not well known because its name makes little sense. No-one can tell from the name that it has anything to do with permutations. I think that Donald Knuth came up with the name.

## Misplaced parenthesis...

@charlie_fcl You need to replace (ifelse with ifelse(.

## But be wary of variable re-use...

@dharr I'm writing this for the OP's benefit; I know that you understand already.

@steveshaver5 When a for loop is finished, its control variable (the i or j) will still be assigned a value. In this case, that value will be 6 (one more than the final value used). Since you've used j as a symbolic matrix entry, this will cause problems. It won't cause Maple to throw an error; but your computation will be incorrect. So, use different loop control variables, such as ii and jj.

This problem doesn't happen with the second method shown by @dharr ; the (i, j) being on the left side of the arrow -> makes them distinct variables that have no effect on previously used i and j.

## Sets that are members of sets...

@vv However, I think that it's fairly clear (although not totally obvious) that an operation called flattening is intended to unbundle only objects (say, sets or lists) that are themselves top-level members of objects of the same type, e.g., sets that are elements of sets.

## for-do...

@Andiguys Your syntax error is not specific to nested loops. Each line that begins with for should end with do.

I'm working on that graph.

## True worksheet, not document...

@Paras31 My previous worksheet may have actually been a document in disguise (because although I removed all your 2D Input and replaced it with 1D input, it began its life as your uploaded document), and that may be why (totally guessing here) the plots didn't display. But the attachment below is a pure typed-from-scratch1D input worksheet.

Please download-and-open this worksheet directly in your Maple. Then execute it in its entirety by clicking on !!! on the toolbar. Then tell me if you get the same results as displayed below.

As always, plots will be crisper in your worksheet than the way they get rendered on MaplePrimes (which is why I usually copy-and-paste plots to MaplePrimes rather than post worksheets). In particular, the 3D plot in your "Task 2" will show fine detail rather than being a "blob".

Note that this worksheet computes and plots about 130,000 points in 6 plots all in under half a second.

 > restart:
 > #Record starting times so that overall time can be measured at the end of worksheet: (st, str):= (time, time[real])():
 > v:= (x,y): V:= ((X,Y):= v(t)): #abbreviations for dependent variables sys:= {     diff(ln~([V]), t)[] =~     (r - b*X - Y*(c + beta/(a+X)), beta/(a/X+1) - mu) }; ics:= {v(0)=~ (0.2, 0.05)}: Param:=  [r, b, c,    a,    mu,  beta]: ParamV:= [1, 1, 0.01, 0.36, 0.4, 0.75]: sol:= dsolve(     sys union ics, {V}, parameters= Param, numeric, maxfun= 0,     range= 0..20000, (abserr, relerr)=~ 1e-4 ):

 > sol(parameters= ParamV); #Integrating to the end of the range at the start is sufficient to suppress #warnings: sol(20000):

 > Opts2D:=     axes= boxed, gridlines, labelfont= [times, 16, bold], size= 100*[8,5],     titlefont= [helvetica, 18],     legendstyle= [location= right, font= [times, 16]] : P:= plots:
 > TrajPlot:= sol->     P:-odeplot(         sol, `[]`~(t, [V]), t= 0..400, color= ["Blue", "Red"], refine= 1,         legend= [V], labels= [`t\n`, ``],         title= "Trajectories         \n", Opts2D, _rest     ) : TrajPlot(sol);

 > #Number of computed points in plots: Npts:= P-> `Point-data dimensions` = [rtable_dims]~([indets(P, rtable)[]]): Npts(%%);

 > Opts3D:= labelfont= [times, 18, bold]: Plot3d:= sol->     P:-odeplot(         sol, [t,V], t= 0..20000, refine= 1, labels= [t,v], Opts3D, _rest     ) : Plot3d(sol); Npts(%);

 > eq_points_sym:= solve(eval(sys, diff= 0), {V});

 > eq_points:= eval([eq_points_sym], Param=~ ParamV);

 > PhasePlot:= sol->     P:-odeplot(         sol, [V], t= 0..20000, refine= 1, color= red, thickness= 2,         title= "Phase-space Diagram\n", labels= [V], Opts2D, _rest     ) : PhasePlot(sol); print(); Npts(%);

 > ParamV:= [1, 1, 0.01, 0.36, 0.4, 1]:
 > sol(parameters= ParamV); sol(20000):

 > TrajPlot(sol, thickness= 0.8); Npts(%);

 > Plot3d(sol, thickness= 0.2); Npts(%);

 > eq_points:= eval([eq_points_sym], Param=~ ParamV);

 > P:-display(     PhasePlot(sol, legend= ``(V), thickness= 0.7),     plot(         eval([[V]], eq_points[-1]), style= point, color= purple,         legend= `Eq. pt.`, symbol= solidcircle, symbolsize= 16     ),     scaling= constrained,     title= "Phase-space diagram         \nand Equilibrium point         \n" );

 > (time, time[real])() -~ (st,str);

 >