Carl Love

Carl Love

27599 Reputation

25 Badges

12 years, 64 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

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

@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.

@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.

@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.

@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
):

{(diff(x(t), t))/x(t) = r-b*x(t)-y(t)*(c+beta/(a+x(t))), (diff(y(t), t))/y(t) = beta/(a/x(t)+1)-mu}

sol(parameters= ParamV);

#Integrating to the end of the range at the start is sufficient to suppress
#warnings:
sol(20000):

[r = 1., b = 1., c = 0.1e-1, a = .36, mu = .4, beta = .75]

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(%%);

`Point-data dimensions` = [[1 .. 312, 1 .. 2], [1 .. 312, 1 .. 2]]

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(%);

`Point-data dimensions` = [[1 .. 8617, 1 .. 3]]

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

{x(t) = mu*a/(beta-mu), y(t) = -a*(a*b*mu-beta*r+mu*r)/(a*beta*c-a*c*mu+beta^2-2*beta*mu+mu^2)}

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

[{x(t) = .4114285714, y(t) = .5992243051}]

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(%);

`Point-data dimensions` = [[1 .. 8617, 1 .. 2]]

Task 2

ParamV:= [1, 1, 0.01, 0.36, 0.4, 1]:

sol(parameters= ParamV);
sol(20000):

[r = 1., b = 1., c = 0.1e-1, a = .36, mu = .4, beta = 1.]

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

`Point-data dimensions` = [[1 .. 1129, 1 .. 2], [1 .. 1129, 1 .. 2]]

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

`Point-data dimensions` = [[1 .. 55673, 1 .. 3]]

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

[{x(t) = .2400000000, y(t) = .4532803181}]

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);

.453, 1.115

 

Download OdeplotRefine.mw

@Kitonum I understand all that you did except for Step 4. How do you justify subtracting the angle from Pi?

@minhthien2016 Yes, I'm considering half planes. If the apex S is on the z-axis, and the base vertices are the vertices of the unit square (as we've all set up the problem), then the angle between any two faces that meet at an edge is at most Pi/2.

@minhthien2016 According to that Wikipedia page, a dihedral angle of a polyhedron is always an interior angle. I stand by the correctness of my Answer: Pi/3 and arccos(1/sqrt(3)).

@Paras31 Did you execute the code by downloading the worksheet from my Answer, or by copy-and-pasting the code from my Answer? I'm wondering if the problem is a 2D Input / Document Mode anomaly, things that I have very little experience with (or patience for).

@Andiguys The \n is the newline control character, I guess it doesn't work in 2D Input, which I never use. It's not a year-version issue, because I've been using it in Maple for over 20 years. It's purpose was to give extra space between the legend box and the w. You can just remove it: Change `w\n` to 'w'. Note that I also changed the quotes from backward to forward.

@emendes I'm not sure if you're asking a followup question, or you were about to and then figured it out. Does my selector (the one with And) work for you?

@Andiguys You could replace `πm` and `Πm` with pi__m and PI__m.

Here's how to do the plot (of course, there are a vast number of stylistic possibilities for presenting the same data):

W:= <seq(30..95)>:
TRC__min:= <seq(s[w][1], w= W)>/10^10:
Tau1__min:= <seq(eval(tau1, s[w][2]), w= W)>:
Tau2__min:= <seq(eval(tau2, s[w][2]), w= W)>:
plot(
    [<W | Tau1__min>, <W | Tau2__min>, <W | TRC__min>], 
    legend= [tau1__min, tau2__min, 'TRC__min'*10^``(-10)],
    labels= [`w\n`, ``], axes= boxed, thickness= 4
);

This plot uses the data from your original worksheet in this thread. If you've subsequently made changes that would change the numeric value of wtau1tau2, or TRC, then those new numbers aren't shown in this plot.

@emendes To select procedures whose names begin with d_, change procedure to 

And(procedure, suffixed('d_'))

in the select command.

@JAMET Are you incapable of making even the slightest generalization from my "advice"? The code in your Question contains 3 occurences of \\nAll three need to be changed to \n.

@nm The help page ?assuming says that (I'm paraphrasing) assuming a property in general without attaching it to specific subexpressions applies the property to all names in the expression. The problem is that y(x) is not considered a name, (but it is considered a function). This works:

simplify(residual) assuming positive, y(x)::positive

First 6 7 8 9 10 11 12 Last Page 8 of 703