dharr

Dr. David Harrington

5799 Reputation

21 Badges

19 years, 241 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

Here's a variation of @mmcdara's idea for a tetrahedron, since he challenged me :). Would require significantly more work for shapes where the folds aren't all on a flat base.

restart;

with(geom3d): with(GraphTheory):

Example - tetrahedron

solid := tetrahedron(t, point(o, [0,0,0]), 3);
f:=faces(solid);
pts:=vertices(solid);
seq(point(cat(p,i),pts[i]),i=1..nops(pts));

t

[[[3^(1/2), 3^(1/2), 3^(1/2)], [3^(1/2), -3^(1/2), -3^(1/2)], [-3^(1/2), 3^(1/2), -3^(1/2)]], [[3^(1/2), 3^(1/2), 3^(1/2)], [-3^(1/2), -3^(1/2), 3^(1/2)], [3^(1/2), -3^(1/2), -3^(1/2)]], [[3^(1/2), 3^(1/2), 3^(1/2)], [-3^(1/2), 3^(1/2), -3^(1/2)], [-3^(1/2), -3^(1/2), 3^(1/2)]], [[3^(1/2), -3^(1/2), -3^(1/2)], [-3^(1/2), -3^(1/2), 3^(1/2)], [-3^(1/2), 3^(1/2), -3^(1/2)]]]

[[3^(1/2), 3^(1/2), 3^(1/2)], [3^(1/2), -3^(1/2), -3^(1/2)], [-3^(1/2), 3^(1/2), -3^(1/2)], [-3^(1/2), -3^(1/2), 3^(1/2)]]

p1, p2, p3, p4

Code the faces with the indices of the points.

indextable:=table(pts=~[$1..nops(pts)]):
ff:=subsindets(f,[algebraic$3],x->indextable[x]);

[[1, 2, 3], [1, 4, 2], [1, 3, 4], [2, 4, 3]]

Make corresponding graph

triangles:=map(CycleGraph,ff):
tetgraph:=Graph([$1..nops(pts)],`union`(map(Edges,triangles)[])):
SetVertexPositions(tetgraph,evalf(pts)):

Try some cuts (red edges on drawing) - spanning tree works for the convex platonic solids

tree:=SpanningTree(tetgraph):
HighlightSubgraph(tetgraph,tree);

DrawGraph(tetgraph, dimension = 3, size = [300, 300],orientation=[50,70]); #red edges are cuts

Uncut edges are the folds, which in this case are all in the same plane.

folds:=Edges(tetgraph) minus Edges(tree);
basevertices:=convert(`union`(%[]),list);
plane(base,map2(cat,p,basevertices)):
triangle(basetriangle,map2(cat,p,basevertices)):

{{2, 3}, {2, 4}, {3, 4}}

[2, 3, 4]

Remove base from the other faces

tofold:=remove(x->{x[]}={basevertices[]},ff);

[[1, 2, 3], [1, 4, 2], [1, 3, 4]]

Unfold

for m to nops(folds) do
 fold:=folds[m];
 foldface:=select(x-> fold subset {x[]}, tofold)[]; # find face to fold
 facepts:=map2(cat,p,foldface); # points in face
 angle:=evalf(FindAngle(base,plane(pl,facepts))); # angle between base and face
 ii := select(has,foldface,fold); # fold sorted as in tofold
 line(cat(raxis,m),map2(cat,p,ii)); # fold is rotation axis
 rotation(cat(unfolded,m),triangle(cat(tr,m),facepts),Pi-angle,cat(raxis,m)); # and rotate
end do:
 

Unfolded

draw([base,unfolded1,unfolded2,unfolded3]);

Angles all the same so generate animation

n:=20:
angles:=[seq(0..Pi-angle,numelems=n)]:

for i to nops(angles) do
  ang:=angles[i];
  rotation(rot1,tr1,ang,raxis1);
  rotation(rot2,tr2,ang,raxis2);
  rotation(rot3,tr3,ang,raxis3);
  dr[i]:=draw([solid,rot1,rot2,rot3]);
end do:
 

plots:-display([seq(dr[i],i=1..n)],insequence=true)

NULL

Download tetrahedron.mw

Here's a more matrixy way.

MM:=Matrix(2,[0,1,1,0]);
eq14 := lhs(eq13) = MM% . simplify(eval(MM^(-1) . rhs(eq13), `%.` = `.`));

ras5_v1.mw

Edit: simpler is 

eq14 := lhs(eq13) = MM %. simplify(value(MM^(-1).rhs(eq13)));

 

Try

p := plot(x^2):
im:=convert(p,Image):
ImageTools:-Embed(im);

 

Use dchange for this. I did the equations separately, but you can put them in a list and do it all in one step.

restart

_local(gamma)

de1 := diff(N(T), T) = R*N(T)*(1-N(T)/K)-alpha*N(T)*P(T)/(A+N(T)); de2 := diff(P(T), T) = gamma*N(T)*P(T)/(A+N(T))+C*P(T)/(1+Q*P(T))-M*P(T); ic1 := N(0) = N_0; ic2 := P(0) = P_0

diff(N(T), T) = R*N(T)*(1-N(T)/K)-alpha*N(T)*P(T)/(A+N(T))

diff(P(T), T) = gamma*N(T)*P(T)/(A+N(T))+C*P(T)/(1+Q*P(T))-M*P(T)

N(0) = N_0

P(0) = P_0

Desired new variables and parameters in terms of the old variables

itr := {a = A/K, b = gamma/R, c = C/R, m = M/R, q = Q*R*K/alpha, t = R*T, x(t) = N(T)/K, y(t) = alpha*P(T)/(R*K)}

{a = A/K, b = gamma/R, c = C/R, m = M/R, q = Q*R*K/alpha, t = R*T, x(t) = N(T)/K, y(t) = alpha*P(T)/(R*K)}

Maple needs the old variables in terms of the new variables.

tr := solve(itr, {A, C, M, Q, T, gamma, N(T), P(T)})

{A = a*K, C = c*R, M = m*R, Q = q*alpha/(R*K), T = t/R, gamma = b*R, N(T) = x(t)*K, P(T) = y(t)*R*K/alpha}

Try it directly

PDEtools:-dchange(tr, de1, [t, x, y, a, b, c, q, m])

(diff(x(t), t))*K*R = R*x(t)*K*(1-x(t))-x(t)*K^2*y(t)*R/(a*K+x(t)*K)

We see we have to divide each side by K*R. Adding simplify simplifies the rhs, but could do this as a separate step.

ode1 := PDEtools:-dchange(tr, de1/(K*R), [t, x, y, a, b, c, q, m], simplify)

diff(x(t), t) = x(t)*(1-x(t)-y(t)/(a+x(t)))

ode2 := PDEtools:-dchange(tr, alpha*de2/(R^2*K), [t, x, y, a, b, c, q, m], simplify)

diff(y(t), t) = y(t)*(b*x(t)/(a+x(t))+c/(1+q*y(t))-m)

init1 := PDEtools:-dchange(tr, ic1/K, [t, x, y, a, b, c, q, m], simplify)

x(0) = N_0/K

init2 := PDEtools:-dchange(tr, alpha*ic2/(K*R), [t, x, y, a, b, c, q, m], simplify)

y(0) = alpha*P_0/(R*K)

NULL

Download dchange.mw

y := 2*hypergeom([1/2, -k/2], [1 - k/2], csc(omega*T)^2)

Depending on what you know about k, you may be able to convert it to a nice form. For example, for k=3

convert(eval(y, k = 3), StandardFunctions)

gives

Download hypergeom.mw

implicitplot3d makes use of the ISOSURFACE structure, see ?plot,structure. I believe that this cannot be coloured with a color function, only GRID and MESH structures:

"For specification of multiple colors, the COLOR(t, A) structure, where t is RGB, HSV or HUE, is similar to that for the 2-D case, except that m by n GRID and MESH structures  are also supported. In this situation, A must be an Array with datatype=float[8] and dimensions 1..m, 1..n when t is HUE, or 1..m, 1..n, 1..3 otherwise."

Here's a workaround and how to use a colorfunction of 3 variables (simple in Maple 2016 and after)

restart

Desired function in implicit form

C := 2*x*y*z-x^2-y^2-z^2+1;

2*x*y*z-x^2-y^2-z^2+1

Desired colour function in 3-variable form

colfunc := (x, y, z) -> (x^2 + y^2 + z^2)/3:

Cannot use a colour function for implicit plot
One (ugly) workaround is to cast C in explicit form, handling the two pieces separately.
(Numeric option not required for operator plot calls as here)

z1,z2:=solve(C,z);
f1:=unapply(z1,x,y,numeric):
f2:=unapply(z2,x,y,numeric):

x*y+(x^2*y^2-x^2-y^2+1)^(1/2), x*y-(x^2*y^2-x^2-y^2+1)^(1/2)

Find the two corresponding colorfunctions

colfunc1:=(x,y)->colfunc(x,y,f1(x,y));
colfunc2:=(x,y)->colfunc(x,y,f2(x,y));
NULL

proc (x, y) options operator, arrow; colfunc(x, y, f1(x, y)) end proc

proc (x, y) options operator, arrow; colfunc(x, y, f2(x, y)) end proc

plot3d([f1,f2],-1..1,-1..1,color=[colfunc1,colfunc2]);

For Maple 2016 and later, colorscheme can be used

plot3d([f1,f2],-1..1,-1..1,colorscheme=["xyzcoloring",colfunc]);

NULL

Download 3d_coloring.mw

They are not to any order, but they are what I would expect. If I simplify by choosing specific examples for the initial functions, then play around with the expansion order, I find the total degrees in the residuals range from order-2 to 2*order-1. Since you have terms with a function times a first derivative of a function, you can expect up to (order plus (order-1)). For the low end, differentiation reduces the order and you have some third-order derivatives. So I would say the answer is verified up to order-3.

restart:

with(DETools):

PDE1 := diff(eta(t,x),t) + 1/2*diff(u(t,x),x) + 1/2*eta(t,x)*diff(u(t,x),x) - 1/48*diff(u(t,x),x$3) + diff(eta(t,x),x)*u(t,x);

diff(eta(t, x), t)+(1/2)*(diff(u(t, x), x))+(1/2)*eta(t, x)*(diff(u(t, x), x))-(1/48)*(diff(diff(diff(u(t, x), x), x), x))+(diff(eta(t, x), x))*u(t, x)

PDE2 := diff(u(t,x),t) + u(t,x)*diff(u(t,x),x) + diff(eta(t,x),x,t,t) + diff(eta(t,x),x) - 1/6*diff(u(t,x),x,x,t);

diff(u(t, x), t)+u(t, x)*(diff(u(t, x), x))+diff(diff(diff(eta(t, x), t), t), x)+diff(eta(t, x), x)-(1/6)*(diff(diff(diff(u(t, x), t), x), x))

sys := rifsimp([PDE1, PDE2]);

table( [( Solved ) = [diff(diff(diff(eta(t, x), t), t), x) = -u(t, x)*(diff(u(t, x), x))-(diff(u(t, x), t))+(1/6)*(diff(diff(diff(u(t, x), t), x), x))-(diff(eta(t, x), x)), diff(diff(diff(u(t, x), x), x), x) = 24*eta(t, x)*(diff(u(t, x), x))+48*(diff(eta(t, x), x))*u(t, x)+48*(diff(eta(t, x), t))+24*(diff(u(t, x), x))] ] )

id := initialdata(sys[Solved]);

table( [( Finite ) = [], ( Infinite ) = [eta(t, x[0]) = _F1(t), (D[2](eta))(t[0], x) = _F2(x), (D[1, 2](eta))(t[0], x) = _F3(x), u(t, x[0]) = _F4(t), (D[2](u))(t, x[0]) = _F5(t), (D[2, 2](u))(t, x[0]) = _F6(t)] ] )

Take a specific example

id2:=eval(eval(id),{_F1(t)=t,_F2(x)=0,_F3(x)=1,_F4(t)=1,_F5(t)=0,_F6(t)=0});

table( [( Finite ) = [], ( Infinite ) = [eta(t, x[0]) = t, (D[2](eta))(t[0], x) = 0, (D[1, 2](eta))(t[0], x) = 1, u(t, x[0]) = 1, (D[2](u))(t, x[0]) = 0, (D[2, 2](u))(t, x[0]) = 0] ] )

sols := rtaylor(sys[Solved], id2, point=[t = 0, x = 0], order = 6);

[eta(t, x) = t+t*x-(1/6)*t^3*x+2*t^2*x^2-4*t^2*x^3-(1/3)*t^4*x^2+(1/120)*t^5*x+(2/9)*t^3*x^3+(22/3)*t^2*x^4, u(t, x) = 1+(112/5)*t*x^5+7*t^2*x^4-(4/3)*t^3*x^3+2*x^4+(48/5)*x^5+8*x^3+8*t*x^3+(8/5)*x^6]

Total degrees of terms in residual, and min and max of these

pdetest(sols,PDE1):
map(degree,[op(%)]);
min(%),max(%);

[5, 5, 4, 5, 5, 11, 11, 8, 10, 10, 10, 9, 9, 11, 11, 10, 9, 11, 9, 11, 11, 9, 7, 6, 7, 8, 6, 8, 7, 8, 7, 6, 6, 4, 5, 5]

4, 11

pdetest(sols,PDE2):
map(degree,[op(%)]);
min(%),max(%);

[5, 5, 5, 4, 5, 11, 10, 10, 9, 9, 10, 9, 11, 11, 9, 11, 11, 11, 7, 8, 8, 8, 7, 6, 4, 5, 5, 11, 10, 9, 8, 7, 6]

4, 11

``

Download MinWorkingExa.mw

You could use the userinfo mechanism to do this. For example, 

pkg := module() export squareme; option package; local ModuleLoad;
  ModuleLoad := proc() userinfo(2, pkg, "Hi I'm loading"); end proc;
  squareme := proc(x) x^2 end proc;
end module; 

which leads, with infolevel[pkg] >= 2:

restart;

infolevel[pkg]:=3;

3

libname:=libname,interface(worksheetdir):

 

with(pkg);

ModuleLoad: Hi I'm loading

[squareme]

``

Download pkgload.mw

You probably know which signs of c__1, c__2 etc to choose the right solution (nicest form of solutions?).

restart;

with(DEtools):

ode1 := diff(c(T),T)=-2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2);
ode2 := diff(p__c(T),T)=2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T);
sys := [ode1,ode2]:

diff(c(T), T) = -2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2)

diff(p__c(T), T) = 2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T)

First solve the system with beta__c = 0.
We want the solution for non-zero beta__c to reduce to this in the limit of beta__c -> 0.

sys0:=eval(sys,beta__c=0);
dsolve(sys0);

[diff(c(T), T) = -2*c(T), diff(p__c(T), T) = 2*p__c(T)]

{c(T) = c__2*exp(-2*T), p__c(T) = c__1*exp(2*T)}

Now solve the system, with solving order [c(T), p__c(T)] (see ?dsolve,system), where the last one is solved first. Because we didn't use the explicit option, a partial solution is given.

To get the full solution we take any of the first solutions for p__c(T) and substitute into any of the second de's for c(T) - this is what option explicit does.

The de's for c(T) have beta__c in the denominator, which seems to be a problem, but let's proceed.

sol1 := dsolve(sys,[c(T),p__c(T)]);

[{p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2)^(1/4)}, {c(T) = -(1/2)*2^(1/2)*(beta__c*p__c(T)*(diff(p__c(T), T)-2*p__c(T)))^(1/2)/beta__c, c(T) = (1/2)*2^(1/2)*(beta__c*p__c(T)*(diff(p__c(T), T)-2*p__c(T)))^(1/2)/beta__c}]

We see that for c__2 = 0, the p__c(T) solutions reduces to the beta__c = 0 one, so we can just replace c__2 by c__2*beta (or c__2*beta^2, etc) to get solutions with the right limiting behavior.
(Perhaps c__2 = beta__c is OK? Looks simpler).

sol2:=[eval(dsolve(sys,[c(T),p__c(T)],explicit),c__2=c__2*beta__c)];

[{c(T) = -(1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -(1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*(-32*beta__c^2*c__2/(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -2*(-2*beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*(-64*beta__c^2*c__2/(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}]

map(odetest,sol2,sys);

[[0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]]

The denominator issue disappears for beta__c>0.

sol3:=simplify(sol2) assuming beta__c::positive;

[{c(T) = -2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}]

Check the limiting behavior - simplify symbolic is reckless here; should consider each case for different signs of c__1, c__2 etc.

simplify(eval(sol3,beta__c=0),symbolic);

[{c(T) = -2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = 2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = ((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = 2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = ((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -(2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -(1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = (2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -(1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -(2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = (1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = (2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = (1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}]

 

Download dsolve.mw

I think I understand better now. You never need to remove anything, just add the latest model to the list if it doesn't match any of the others, or add it as a twin if it does. If that is correct, then this works. Since lists are inefficient, this uses a table containing Arrays for the twins and a Vector (could be an Array) for the list.

KeyModels := Vector(): # unique models (was list l)
Models := table():     # equivalent models (was list of twins)
nkeys := 0:
for n to nmodels do
    # see which key model the new model n is equivalent to
    # by comparing with the existing key models
    # here do it by a random number generator
    r := rand(1 .. nkeys + 1)();
    if r = nkeys + 1 then   # didn't match any
        nkeys := nkeys + 1;
        # new model - record it as key
        # and start a new Array
        KeyModels ,= n; # append to key models
        Models[n] := Array([n]); # or perhaps Array([])
    else                    # matched existing model KeyModelNum
        KeyModelNum := KeyModels[r];
        # add to existing models
        Models[KeyModelNum] ,= n;
    end if;
end do:

restart;

nmodels:=100000;

100000

st:=time():

KeyModels := Vector(): # unique models (was list l)

time()-st;

4.813

Results

KeyModels;

_rtable[36893489746279458924]

Models[4];
Models[28];

Models[4]

Array(%id = 36893489746279460364)

NULL

Download twins.mw

Notes: In principle, you could skip the KeyModels Vector and just use for i,_ in eval(Models) do to iterate through the keys (but not in numerical order). I've anticipated you might want the key model in with the twins, but if not replace Array([n]) with Array([]).

Maple does not properly export .eps for 3D plots (works fine with vector graphics for 2D plots). You can export .pdf from the plot context menu and it will work directly in overleaf with the graphicx package. If you want to do it programmatically you can use the jpeg or png driver, but of course these are bitmap formats (as is the .pdf)..

Most things in linalg have their equivalent in LinearAlgebra, so I wouldn't say it is split in two. Note that the ArrayTools package is restricted to rectangular storage array objects (at least in 2023). For these, ArrayTools:-IsZero is fast.

For sparse (or rectangular) matrices, M-> andmap(`=`,M,0) is a simple alternative.

[Edit: see @sursumCorda's response - andmap is only good for moderately sized rectangular matrices]

The non-dim Lambda is a function of 2 non-dim variables, so you will need 3D plots this time.

[Edit in blue - the choice of variables to eliminate was not working; was always the first 2]

restart;

local gamma:with(LinearAlgebra):

Want to "non-dimensionalize" this expression.

p3 := RootOf((64*gamma^2*sigma__d1^2*sigma__d2^6*sigma__v^2 + 64*gamma^2*sigma__d2^8*sigma__v^2 + 512*sigma__d2^6)*_Z^10 - 2*gamma^6*sigma__d1^2*sigma__v^6 - 8*gamma^5*sigma__d1^2*sigma__v^5*_Z + (-3*gamma^6*sigma__d1^4*sigma__v^6 - 3*gamma^6*sigma__d1^2*sigma__d2^2*sigma__v^6 + 24*gamma^4*sigma__d2^2*sigma__v^4)*_Z^2 + (-12*gamma^5*sigma__d1^4*sigma__v^5 - 52*gamma^5*sigma__d1^2*sigma__d2^2*sigma__v^5 - 48*gamma^5*sigma__d2^4*sigma__v^5 + 32*gamma^3*sigma__d1^2*sigma__v^3 + 160*gamma^3*sigma__d2^2*sigma__v^3)*_Z^3 + (12*gamma^6*sigma__d1^4*sigma__d2^2*sigma__v^6 + 36*gamma^6*sigma__d1^2*sigma__d2^4*sigma__v^6 + 24*gamma^6*sigma__d2^6*sigma__v^6 - 12*gamma^4*sigma__d1^4*sigma__v^4 - 216*gamma^4*sigma__d1^2*sigma__d2^2*sigma__v^4 - 412*gamma^4*sigma__d2^4*sigma__v^4 + 32*gamma^2*sigma__d1^2*sigma__v^2 + 384*gamma^2*sigma__d2^2*sigma__v^2)*_Z^4 + (32*gamma^5*sigma__d1^4*sigma__d2^2*sigma__v^5 + 224*gamma^5*sigma__d1^2*sigma__d2^4*sigma__v^5 + 248*gamma^5*sigma__d2^6*sigma__v^5 - 336*gamma^3*sigma__d1^2*sigma__d2^2*sigma__v^3 - 1392*gamma^3*sigma__d2^4*sigma__v^3 + 384*gamma*sigma__d2^2*sigma__v)*_Z^5 + (4*gamma^6*sigma__d1^6*sigma__d2^2*sigma__v^6 + 16*gamma^6*sigma__d1^4*sigma__d2^4*sigma__v^6 + 16*gamma^6*sigma__d1^2*sigma__d2^6*sigma__v^6 + 4*gamma^6*sigma__d2^8*sigma__v^6 + 16*gamma^4*sigma__d1^4*sigma__d2^2*sigma__v^4 + 512*gamma^4*sigma__d1^2*sigma__d2^4*sigma__v^4 + 1072*gamma^4*sigma__d2^6*sigma__v^4 - 176*gamma^2*sigma__d1^2*sigma__d2^2*sigma__v^2 - 2288*gamma^2*sigma__d2^4*sigma__v^2 + 128*sigma__d2^2)*_Z^6 + (48*gamma^5*sigma__d1^4*sigma__d2^4*sigma__v^5 + 96*gamma^5*sigma__d1^2*sigma__d2^6*sigma__v^5 + 32*gamma^5*sigma__d2^8*sigma__v^5 + 512*gamma^3*sigma__d1^2*sigma__d2^4*sigma__v^3 + 2464*gamma^3*sigma__d2^6*sigma__v^3 - 1792*gamma*sigma__d2^4*sigma__v)*_Z^7 + (32*gamma^4*sigma__d1^4*sigma__d2^4*sigma__v^4 + 208*gamma^4*sigma__d1^2*sigma__d2^6*sigma__v^4 + 96*gamma^4*sigma__d2^8*sigma__v^4 + 192*gamma^2*sigma__d1^2*sigma__d2^4*sigma__v^2 + 3136*gamma^2*sigma__d2^6*sigma__v^2 - 512*sigma__d2^4)*_Z^8 + (192*gamma^3*sigma__d1^2*sigma__d2^6*sigma__v^3 + 128*gamma^3*sigma__d2^8*sigma__v^3 + 2048*gamma*sigma__d2^6*sigma__v)*_Z^9):

First write the polynomial with variable lambda, and normalize so one of the terms is 1 (arbitrarily take the first one)

p:=expand(subs(_Z=lambda,op(p3))):
pnorm:=expand(%/op(1,%));

1+4*sigma__d2^2/sigma__d1^2+4*sigma__d2^4/sigma__d1^4+sigma__d2^6/sigma__d1^6+12*lambda*sigma__d2^2/(gamma*sigma__d1^2*sigma__v)+24*lambda*sigma__d2^4/(gamma*sigma__d1^4*sigma__v)+8*lambda*sigma__d2^6/(gamma*sigma__d1^6*sigma__v)+8*lambda^2*sigma__d2^2/(gamma^2*sigma__d1^2*sigma__v^2)+52*lambda^2*sigma__d2^4/(gamma^2*sigma__d1^4*sigma__v^2)+24*lambda^2*sigma__d2^6/(gamma^2*sigma__d1^6*sigma__v^2)+48*lambda^3*sigma__d2^4/(gamma^3*sigma__d1^4*sigma__v^3)+32*lambda^3*sigma__d2^6/(gamma^3*sigma__d1^6*sigma__v^3)+3/(lambda^2*sigma__d1^2)+9*sigma__d2^2/(lambda^2*sigma__d1^4)+6*sigma__d2^4/(lambda^2*sigma__d1^6)+16*lambda^4*sigma__d2^4/(gamma^4*sigma__d1^4*sigma__v^4)+16*lambda^4*sigma__d2^6/(gamma^4*sigma__d1^6*sigma__v^4)+8/(gamma*lambda*sigma__d1^2*sigma__v)+56*sigma__d2^2/(gamma*lambda*sigma__d1^4*sigma__v)+62*sigma__d2^4/(gamma*lambda*sigma__d1^6*sigma__v)+4/(gamma^2*sigma__d1^2*sigma__v^2)+128*sigma__d2^2/(gamma^2*sigma__d1^4*sigma__v^2)+268*sigma__d2^4/(gamma^2*sigma__d1^6*sigma__v^2)+128*lambda*sigma__d2^2/(gamma^3*sigma__d1^4*sigma__v^3)+616*lambda*sigma__d2^4/(gamma^3*sigma__d1^6*sigma__v^3)-(3/4)/(lambda^4*sigma__d1^2*sigma__d2^2)-(3/4)/(lambda^4*sigma__d1^4)+48*lambda^2*sigma__d2^2/(gamma^4*sigma__d1^4*sigma__v^4)+784*lambda^2*sigma__d2^4/(gamma^4*sigma__d1^6*sigma__v^4)-3/(gamma*lambda^3*sigma__d1^2*sigma__d2^2*sigma__v)-13/(gamma*lambda^3*sigma__d1^4*sigma__v)-12*sigma__d2^2/(gamma*lambda^3*sigma__d1^6*sigma__v)+512*lambda^3*sigma__d2^4/(gamma^5*sigma__d1^6*sigma__v^5)-3/(gamma^2*lambda^2*sigma__d1^2*sigma__d2^2*sigma__v^2)-54/(gamma^2*lambda^2*sigma__d1^4*sigma__v^2)-103*sigma__d2^2/(gamma^2*lambda^2*sigma__d1^6*sigma__v^2)+128*lambda^4*sigma__d2^4/(gamma^6*sigma__d1^6*sigma__v^6)-84/(gamma^3*lambda*sigma__d1^4*sigma__v^3)-348*sigma__d2^2/(gamma^3*lambda*sigma__d1^6*sigma__v^3)-(1/2)/(lambda^6*sigma__d1^4*sigma__d2^2)-44/(gamma^4*sigma__d1^4*sigma__v^4)-572*sigma__d2^2/(gamma^4*sigma__d1^6*sigma__v^4)-2/(gamma*lambda^5*sigma__d1^4*sigma__d2^2*sigma__v)-448*lambda*sigma__d2^2/(gamma^5*sigma__d1^6*sigma__v^5)+6/(gamma^2*lambda^4*sigma__d1^6*sigma__v^2)-128*lambda^2*sigma__d2^2/(gamma^6*sigma__d1^6*sigma__v^6)+8/(gamma^3*lambda^3*sigma__d1^4*sigma__d2^2*sigma__v^3)+40/(gamma^3*lambda^3*sigma__d1^6*sigma__v^3)+8/(gamma^4*lambda^2*sigma__d1^4*sigma__d2^2*sigma__v^4)+96/(gamma^4*lambda^2*sigma__d1^6*sigma__v^4)+96/(gamma^5*lambda*sigma__d1^6*sigma__v^5)+32/(gamma^6*sigma__d1^6*sigma__v^6)

So all the other terms must now be dimensionless. We use the automated non-dimensionalization procedure of E. Hubert, G, Labahn, Found Comput Math 13 (2013) 479–516, doi:10.1007/s10208-013-9165-9. See also https://youtu.be/Nl2FBAbU1pE

terms:=[op(2..-1,pnorm)]:
allvars:=[indets(p3,name)[],lambda]; #lambda as a true variable is last. A simpler result may occur if the ones to eliminate are first
nvars:=nops(allvars);
nterms:=nops(terms)

[gamma, sigma__d1, sigma__d2, sigma__v, lambda]

5

51

Matrix K contains the exponents of the variables in each term. It has 5 rows but rank 3, so we can eliminate 5 - 3 = 2 variables.

By default these are the first 2 in allvars, so rearrange allvars if you want a different result.

K := Matrix(nvars, nterms, (i, j) -> diff(terms[j], allvars[i])*allvars[i]/terms[j]);
r1 := Rank(K);

_rtable[36893490145425527140]

3

Hubert's matrix A has full row rank, equal to the number of variables to eliminate

H, U := HermiteForm(K, output = ['H', 'U'], method = 'integer'):
Determinant(U):
A := U[-nvars + r1 .. -1, .. ];
A . K: # should be zero
rA := Rank(A);

Matrix(2, 5, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 1, (1, 5) = 0, (2, 1) = 1, (2, 2) = -1, (2, 3) = -1, (2, 4) = 0, (2, 5) = 1})

2

Calculate matrix W. Procedure ColumnHermiteForm is in startup code

HA, V := ColumnHermiteForm(A):
W := V^(-1):

Vi := V[.., 1 .. rA]:
Vn := V[.., -nvars + rA .. -1]:
Wu := W[1 .. rA, ..]:
Wd := W[-nvars + rA .. -1, ..]:

nondims := table():
nondimvars := table():
for i to nvars - rA do
    nondimvars[i] := cat(allvars[i + rA], __new);
    nondims[i] := nondimvars[i] = mul(`^`~(allvars, V[ .. , i + rA]));
end do:
nondimvars := convert(nondimvars, list);
nondims := convert(nondims, list);

[sigma__d2__new, sigma__v__new, lambda__new]

[sigma__d2__new = sigma__d2/sigma__d1, sigma__v__new = gamma*sigma__d1*sigma__v, lambda__new = sigma__d1*lambda]

rewrites := table():
for i to nvars do
    rewrites[i] := allvars[i] = mul(`^`~(nondimvars, Wd[ .. , i]));
end do:
rewrites := convert(rewrites, list);

[gamma = 1, sigma__d1 = 1, sigma__d2 = sigma__d2__new, sigma__v = sigma__v__new, lambda = lambda__new]

So the non-dimensionalized p is

newsys := subs(rewrites, p);

64*lambda__new^10*sigma__d2__new^8*sigma__v__new^2+128*lambda__new^9*sigma__d2__new^8*sigma__v__new^3+96*lambda__new^8*sigma__d2__new^8*sigma__v__new^4+32*lambda__new^7*sigma__d2__new^8*sigma__v__new^5+4*lambda__new^6*sigma__d2__new^8*sigma__v__new^6+64*lambda__new^10*sigma__d2__new^6*sigma__v__new^2+192*lambda__new^9*sigma__d2__new^6*sigma__v__new^3+208*lambda__new^8*sigma__d2__new^6*sigma__v__new^4+96*lambda__new^7*sigma__d2__new^6*sigma__v__new^5+16*lambda__new^6*sigma__d2__new^6*sigma__v__new^6+512*lambda__new^10*sigma__d2__new^6+2048*lambda__new^9*sigma__d2__new^6*sigma__v__new+3136*lambda__new^8*sigma__d2__new^6*sigma__v__new^2+32*lambda__new^8*sigma__d2__new^4*sigma__v__new^4+2464*lambda__new^7*sigma__d2__new^6*sigma__v__new^3+48*lambda__new^7*sigma__d2__new^4*sigma__v__new^5+1072*lambda__new^6*sigma__d2__new^6*sigma__v__new^4+16*lambda__new^6*sigma__d2__new^4*sigma__v__new^6+248*lambda__new^5*sigma__d2__new^6*sigma__v__new^5+24*lambda__new^4*sigma__d2__new^6*sigma__v__new^6+192*lambda__new^8*sigma__d2__new^4*sigma__v__new^2+512*lambda__new^7*sigma__d2__new^4*sigma__v__new^3+512*lambda__new^6*sigma__d2__new^4*sigma__v__new^4+4*lambda__new^6*sigma__d2__new^2*sigma__v__new^6+224*lambda__new^5*sigma__d2__new^4*sigma__v__new^5+36*lambda__new^4*sigma__d2__new^4*sigma__v__new^6-512*lambda__new^8*sigma__d2__new^4-1792*lambda__new^7*sigma__d2__new^4*sigma__v__new-2288*lambda__new^6*sigma__d2__new^4*sigma__v__new^2+16*lambda__new^6*sigma__d2__new^2*sigma__v__new^4-1392*lambda__new^5*sigma__d2__new^4*sigma__v__new^3+32*lambda__new^5*sigma__d2__new^2*sigma__v__new^5-412*lambda__new^4*sigma__d2__new^4*sigma__v__new^4+12*lambda__new^4*sigma__d2__new^2*sigma__v__new^6-48*lambda__new^3*sigma__d2__new^4*sigma__v__new^5-176*lambda__new^6*sigma__d2__new^2*sigma__v__new^2-336*lambda__new^5*sigma__d2__new^2*sigma__v__new^3-216*lambda__new^4*sigma__d2__new^2*sigma__v__new^4-52*lambda__new^3*sigma__d2__new^2*sigma__v__new^5-3*lambda__new^2*sigma__d2__new^2*sigma__v__new^6+128*lambda__new^6*sigma__d2__new^2+384*lambda__new^5*sigma__d2__new^2*sigma__v__new+384*lambda__new^4*sigma__d2__new^2*sigma__v__new^2-12*lambda__new^4*sigma__v__new^4+160*lambda__new^3*sigma__d2__new^2*sigma__v__new^3-12*lambda__new^3*sigma__v__new^5+24*lambda__new^2*sigma__d2__new^2*sigma__v__new^4-3*lambda__new^2*sigma__v__new^6+32*lambda__new^4*sigma__v__new^2+32*lambda__new^3*sigma__v__new^3-8*lambda__new*sigma__v__new^5-2*sigma__v__new^6

The naming here is not optimal. sigma__v__new is similar to Gamma previously but with sigma__d1 instead of sigma__d, so Gamma or maybe Gamma_1; lambda__new should be Lamba, and the ratio sigma__d2/sigma__d1 could be some other capital Greek letter, say Psi
So the new Lambda is

newnames:={sigma__d2__new = Psi, sigma__v__new=Gamma};
Lambda := RootOf(eval(newsys,newnames),lambda__new);

{sigma__d2__new = Psi, sigma__v__new = Gamma}

RootOf((64*Gamma^2*Psi^8+64*Gamma^2*Psi^6+512*Psi^6)*_Z^10+(128*Gamma^3*Psi^8+192*Gamma^3*Psi^6+2048*Gamma*Psi^6)*_Z^9+(96*Gamma^4*Psi^8+208*Gamma^4*Psi^6+32*Gamma^4*Psi^4+3136*Gamma^2*Psi^6+192*Gamma^2*Psi^4-512*Psi^4)*_Z^8+(32*Gamma^5*Psi^8+96*Gamma^5*Psi^6+48*Gamma^5*Psi^4+2464*Gamma^3*Psi^6+512*Gamma^3*Psi^4-1792*Gamma*Psi^4)*_Z^7+(4*Gamma^6*Psi^8+16*Gamma^6*Psi^6+16*Gamma^6*Psi^4+1072*Gamma^4*Psi^6+4*Gamma^6*Psi^2+512*Gamma^4*Psi^4+16*Gamma^4*Psi^2-2288*Gamma^2*Psi^4-176*Gamma^2*Psi^2+128*Psi^2)*_Z^6+(248*Gamma^5*Psi^6+224*Gamma^5*Psi^4+32*Gamma^5*Psi^2-1392*Gamma^3*Psi^4-336*Gamma^3*Psi^2+384*Gamma*Psi^2)*_Z^5+(24*Gamma^6*Psi^6+36*Gamma^6*Psi^4+12*Gamma^6*Psi^2-412*Gamma^4*Psi^4-216*Gamma^4*Psi^2-12*Gamma^4+384*Gamma^2*Psi^2+32*Gamma^2)*_Z^4+(-48*Gamma^5*Psi^4-52*Gamma^5*Psi^2-12*Gamma^5+160*Gamma^3*Psi^2+32*Gamma^3)*_Z^3+(-3*Gamma^6*Psi^2-3*Gamma^6+24*Gamma^4*Psi^2)*_Z^2-8*Gamma^5*_Z-2*Gamma^6)

Find how to write the old variables in terms of the new ones, and directly use eval as a check. We only need to eval for the variables that we are not eliminating; the others will automatically disappear

map(lhs,remove(x->rhs(x)=1,rewrites)):
OldToNew := eval(solve(nondims,%)[],newnames);
sigma__d1^4*simplify(eval(p, OldToNew)): # may need to scale to get rid of last of the old variables
RootOf(%,lambda__new):
%-Lambda; # should be zero

[sigma__d2 = Psi*sigma__d1, sigma__v = Gamma/(gamma*sigma__d1), lambda = lambda__new/sigma__d1]

0

allvals:=[allvalues(Lambda)]: # just different index values for the RootOf

plot3d(allvals[1],Gamma=0..10,Psi=0..10);

NULL

Download nondim.mw

The default appears to be 500x500 pixels. If you add the plot option size=[1000,1000] then the exported file will be 1000x1000 pixels.

You need to set up a Table of Contents to load into the help database with HelpTools:-TableOfContents:-Load(TOC,helpdbname);

In the TOC, you set a priority for each entry, so make your overview the highest priority

See here on Mapleprimes for an example.

1 2 3 4 5 6 7 Last Page 1 of 62