dharr

Dr. David Harrington

4791 Reputation

21 Badges

19 years, 33 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 replies submitted by dharr

@lcz I agree, my code actually tests for a cutset, which is a minimal cut, but not all cuts are minimal. Apologies; I sorted all this out before, but forgot about this distinction.

"A cutset S of a connected graph G is a minimal set of edges of G such that removal of S [dis]connects G into exactly two components", p.43, K. Thulasiraman and M.N.S. Swamy, "Graphs, Theory and Algorithms", Wiley 1992. doi:10.1002/9781118033104.

The stackexchange algorithm involving contraction detecting extra edges with loops is not easy to implement in Maple with its Contract or ContractSubgraph commands because it doesn't keep track of where loops came from. But one can just check for extra edges that are in one component before contracting. Here it is: 

restart;

IsCut returns true if the edges are a cut

IsCut:=proc(G1::GRAPHLN,edges::set)
       local G,partition,subgraphs,alledges;
       uses GraphTheory;
       if not IsConnected(G1) then error "Graph must be connected" end if;
       G:=UnderlyingGraph(G1); # make undirected, unweighted, without loops
       alledges:=Edges(G);
       partition:=ConnectedComponents(DeleteEdge(G,edges,'inplace'=false));
       if nops(partition) <= 1 then return false end if;
       # check for edges that have both ends in one component
       subgraphs:=map2(InducedSubgraph,G,partition);
       if not (alledges minus `union`(map(Edges,subgraphs)[]) = edges) then return false end if;
       # contract the components and see if the graph is bipartite
       IsBipartite(foldl(ContractSubgraph,G,subgraphs[]));
end proc:

with(GraphTheory):
with(combinat):
G:=CompleteGraph(3,3):
edge:=choose(Edges(G), 7):
`and`(seq(IsCutSet(G,s),s in edge)); # any 7 edges cuts the graph into 2 or more components

true

`or`(seq(IsCut(G,s),s in edge)); # no set of 7 edges is a cut

false

DrawGraph(G,size=[250,250]);

edges:={{1,4},{1,5},{1,6},{3,6}}; # {3,6} edge is extra
IsCut(G,edges);
 

{{1, 4}, {1, 5}, {1, 6}, {3, 6}}

false

edges:={{1,4},{1,5},{1,6}};
IsCut(G,edges);

{{1, 4}, {1, 5}, {1, 6}}

true

edges:={{1,4},{1,5},{1,6},{3,6},{3,5},{3,4}}; #3 components after cut; contracted graph is bipartite
IsCut(G,edges);

{{1, 4}, {1, 5}, {1, 6}, {3, 4}, {3, 5}, {3, 6}}

true

G2:=CycleGraph(3): #3 components, but contracted graph (same as original) is not bipartite
IsCut(G2,Edges(G2));

false

NULL

Download IsCut.mw

@vv @Christian Wolinski Thanks! Interestingly, I was using a more complicated expression earlier, and I got an error message:

with(Logic):
q:=&or((&not x[2]) &and (&not x[8]),(&not x[2]) &and (&not x[true]),&and(x[3],x[9]
,&not x[true]),&and(&not x[3],&not x[9],&not x[true]),&and(x[2],x[3],x[8],x[9])
,&and(x[2],x[8],&not x[3],&not x[9]),&and(x[3],x[true],&not x[8],&not x[9]),
&and(x[9],x[true],&not x[3],&not x[8])):
SymmetryGroup(q);

gives: Error, (in Logic:-SymmetryGroup) not in conjunctive normal form

so it can test. So then I didn't think about the form, which was a mistake. I'll submit an SCR about this (and Normalize)

@mmcdara Not sure what the OP wants; let's see. I would think a cycle in the underlying (undirected) graph could correspond to 0,1, or 2 directed cycles (neither direction works, clockwise works or counterclockwise works but not both, both directions work). Perhaps like this:

restart

with(GraphTheory)

Takes a list of vertices [v1, ..., vn] and returns [v1,...,vn] if v1->v2...->vn->v1 is a directed cycle and/or [vn,...v1] if vn->....->v1->nv is a directed cycle. Returns NULL if neither of these is a directed cycle (includes not a cycle at all).

DirectedCycles:=proc(G::GRAPHLN,vlist::list)
   local rlist,cycle1,cycle2,subedges;
   uses GraphTheory;
   if not IsDirected(G) then return NULL end if;
   subedges:=Edges(InducedSubgraph(G,vlist));
   if Edges(Graph(Trail(vlist[],vlist[1]),directed))
       subset subedges then cycle1:=vlist else cycle1:=NULL end if;
   rlist:=ListTools:-Reverse(vlist);  
   if Edges(Graph(Trail(rlist[],rlist[1]),directed))
       subset subedges then cycle2:=rlist else cycle2:=NULL end if;
   cycle1,cycle2;     
end proc:

G := Graph({[1, 2], [1, 3], [2, 1], [2, 3], [2, 4], [3, 1], [3, 2], [4, 3]})

GRAPHLN(directed, unweighted, [1, 2, 3, 4], Array(1..4, {(1) = {2, 3}, (2) = {1, 3, 4}, (3) = {1, 2}, (4) = {3}}), `GRAPHLN/table/1`, 0)

DrawGraph(G, size = [250, 250], layout = spring)

DirectedCycles(G, [1, 2, 3])

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

DirectedCycles(G, [3, 2, 4])

[3, 2, 4]

DirectedCycles(G, [2, 1, 3, 4])NULL

[4, 3, 1, 2]

NULLDirectedCycles(G, [2, 1, 4])NULL

Find some of the cycles and check them out.

cycles := CycleBasis(UnderlyingGraph(G)); dcycles := map2(DirectedCycles, G, cycles)

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

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

``

Download Cycles.mw

@sursumCorda I guess that's what the help means by "The isolve command has some limited ability to deal with inequalities." In my experience isolve is limited. In this case, as you probably already tried, you can just inelegantly iterate through all the _Z1 etc values, (Iterator:-Multiseq might be a little faster.)

restart;

local gamma:
alias(alpha=RootOf(_Z^3-4*_Z^2+_Z+1,index=1),beta=RootOf(_Z^3-4*_Z^2+_Z+1,index=2),gamma=RootOf(_Z^3-4*_Z^2+_Z+1,index=3));

alpha, beta, gamma

eqs:=eval~(k[1]*x^3+k[2]*x^2*y+k[5]*x^2*z+k[3]*x*y^2+k[6]*x*y*z+k[8]*x*z^2+k[4]*y^3+k[7]*y^2*z+k[9]*y*z^2+k[10]*z^3,{{x=1,y=1,z=1},{x=RootOf(_Z^3-4*_Z^2+_Z+1,index=1),y=RootOf(_Z^3-4*_Z^2+_Z+1,index=2),z=RootOf(_Z^3-4*_Z^2+_Z+1,index=3)},{x=RootOf(_Z^3-4*_Z^2+_Z+1,index=2),y=RootOf(_Z^3-4*_Z^2+_Z+1,index=3),z=RootOf(_Z^3-4*_Z^2+_Z+1,index=1)},{x=RootOf(_Z^3-4*_Z^2+_Z+1,index=3),y=RootOf(_Z^3-4*_Z^2+_Z+1,index=1),z=RootOf(_Z^3-4*_Z^2+_Z+1,index=2)}}):

eqs2:=solve(evala(eqs)):

eqs3:=isolve(eqs2);

{k[1] = _Z1, k[2] = 10*_Z1-_Z2+2*_Z3+_Z4-_Z5-5*_Z6, k[3] = -4*_Z1-4*_Z2-_Z3-_Z5-4*_Z6, k[4] = _Z2, k[5] = _Z3, k[6] = -4*_Z1+3*_Z2-_Z3-3*_Z4+4*_Z5+22*_Z6, k[7] = -3*_Z1+_Z2-_Z3+_Z4-3*_Z5-14*_Z6, k[8] = _Z4, k[9] = _Z5, k[10] = _Z6}

mn:=-5;mx:=5;i:=0:
for Z1 from mn to mx do for Z2 from mn to mx do for Z3 from mn to mx do
 for Z4 from mn to mx do for Z5 from mn to mx do for Z6 from mn to mx do
 s:=eval(eqs3,{_Z1=Z1,_Z2=Z2,_Z3=Z3,_Z4=Z4,_Z5=Z5,_Z6=Z6});
 if andmap(x->evalb(rhs(x)>=mn and rhs(x)<=mx),s) then
   i:=i+1;
   sol[i]:=s;
 end if;
od od od od od od:
i;

-5

5

1407

NULL

Download algeqns.mw

@Christian Wolinski Definitely unexpected to get a simpler RootOf here - documentation for convert,radical says "The conversion can fail if Maple cannot find radical expressions for the roots or if the correct radical expression cannot be selected. If the conversion fails, the RootOf remains unchanged." [my bold].

@sursumCorda If you know which RootOf you want the answer in terms of, then Algfield will help.

restart;

a:=[RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3), -4/5*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3)^2 + 19/5*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3) + 3/5, 1];
b:=[RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3), RootOf(_Z^3 + 10*_Z^2 + 3*_Z - 1, index = 3), 1];
c:=[RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 1)/RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 3), RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 2)/RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 3), 1];

[RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3), -(4/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3)^2+(19/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3)+3/5, 1]

[RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3), RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3), 1]

[RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3), RootOf(_Z^3-4*_Z^2+_Z+1, index = 2)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3), 1]

evalf(a)evalf(b)evalf(c)

[-1.924983964, -9.679389670, 1.]

[-1.924983964, -9.679389672, 1.]

[-1.924983964, -9.679389671, 1.]

From a[2] to b[2]

qb := evala(Algfield(a[2], {RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3)})); evala(eval(a[2], qb[1]))

RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3)

From a to c

qc := evala(Algfield(a, {RootOf(_Z^3-4*_Z^2+_Z+1, index = 1), RootOf(_Z^3-4*_Z^2+_Z+1, index = 2), RootOf(_Z^3-4*_Z^2+_Z+1, index = 3)})); evala(eval(a, qc[1]))

[RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-2*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-1, 3*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-10*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-4, 1]

I think Maple will always prefer a polynomial over a rational function. Notice here that since they can all be written in terms of one index value, so that is simpler. The splitting field for this polynomial has only one algebraic number (by default index=1)

PolynomialTools:-Split(_Z^3-4*_Z^2+_Z+1, _Z, 'K'); K

{RootOf(_Z^3-4*_Z^2+_Z+1)}

But you can check that c[1] is the same as the above result.

evala(RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3))

RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-2*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-1

NULL

Download Algfield.mw

@Zeineb The Routh-Hurwitz 'array'  treatment is equivalent so will lead to the same answer. If you got the polynomial from a matrix, then working from the matrix (or even the original differential equations) can sometimes be easier. For example a symmetric matrix proves real roots, but that is harder to tell from the polynomial.

@Axel Vogt Creative approach! (gfun is now a standard package within Maple.)

@C_R I branched it off to here

@C_R Interesting. Yes, the index=real[4] RootOf does not behave like a typical RootOf in some respects. Normally evalf on a RootOf would give a numerical value but here it returns unevaluated.  This would normally be a way to work toward an exact solution, since convert(RootOf(...,numeric value),radical) could then be used.

@sursumCorda A combination of implicitplot (plot_real_curves is sometimes too slow) and algcurves:-singularities can track down the roots. 

restart

with(algcurves)

expr := a^3*b^3*(12*(a*b+a+b)-28*(a+b+1)^2)+4*a^2*b^2*(a+b+1)*(a^2+b^2+(a+1)*(b+1))*(13*(a+b+1)^2-8*(a*b+a+b))+a*b*(7*(a+b+1)^8-40*(a*b+a+b)*(a+b+1)^6+32*(a*b+a+b)^2*(a+b+1)^4+28*(a*b+a+b)^3*(a+b+1)^2-20*(a*b+a+b)^4)+(a+b+1)*(a*b+a+b)*((a+b+1)^4-6*(a*b+a+b)*(a+b+1)^2+6*(a*b+a+b)^2)^2:

plots:-implicitplot(expr, a = -5 .. 5, b = -5 .. 5)

Curve is symmetric - easy to find solutions with a=b. There are three non-negative ones.

simplify(eval(expr, b = a)); solve(%); evalf(%); solns := {a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)^(1/4), b = (1/2)^(1/4)}

8*(a-1)^2*a*(a^4-1/2)^2

0, 1, 1, (1/2)*2^(3/4), ((1/2)*I)*2^(3/4), -(1/2)*2^(3/4), -((1/2)*I)*2^(3/4), (1/2)*2^(3/4), ((1/2)*I)*2^(3/4), -(1/2)*2^(3/4), -((1/2)*I)*2^(3/4)

0., 1., 1., .8408964155, .8408964155*I, -.8408964155, -.8408964155*I, .8408964155, .8408964155*I, -.8408964155, -.8408964155*I

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}

The plot shows that there are two other solutions along each of the a and b axes - these are the index=1 and index=2 ones

simplify(eval(expr, b = 0)); q := solve(%); evalf(%)

(a+1)*a*(a^4-2*a^3-2*a+1)^2

-1, 0, RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 1), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 2), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 3), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 4), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 1), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 2), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 3), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 4)

-1., 0., .4354205447, 2.296630263, -.3660254038+.9306048591*I, -.3660254038-.9306048591*I, .4354205447, 2.296630263, -.3660254038+.9306048591*I, -.3660254038-.9306048591*I

solns := solns, {a = 0, b = convert(q[3], radical)}, {a = convert(q[3], radical), b = 0}, {a = 0, b = convert(q[4], radical)}, {a = convert(q[4], radical), b = 0}

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}, {a = 0, b = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 0, b = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4), b = 0}

Any other points are not part of the main curve, so are singularities. Here they are in homogeneous coordinates, but since the last coordinate is 1, we can just take the first two coordinates. We already had the a=1,b=1 solution and the a=(1/2)^(1/4),b=(1/2)^(1/4) one

s := singularities(expr, a, b)

{[[1, 1, 1], 2, 1, 2], [[1, RootOf(_Z^4-2), 1], 2, 1, 2], [[RootOf(_Z^4-2), 1, 1], 2, 1, 2], [[RootOf(2*_Z^4-1), RootOf(2*_Z^4-1), 1], 2, 1, 2]}

simplify(eval(expr, {a = 1, b = 2^(1/4)}))

0

So the complete set is

solns := solns, {a = 1, b = 2^(1/4)}, {a = 2^(1/4), b = 1}

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}, {a = 0, b = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 0, b = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 1, b = 2^(1/4)}, {a = 2^(1/4), b = 1}

NULL

Download algcurves.mw

For the second one, there are only singularities, and I'll let you break them down into all the subcases.

algcurves2.mw

@C_R  Interesting about index = real[4]. I've not seen that before in many years using Maple, and the stuff on the forum isn't clear about it being a bug or not. Perhaps it is a bug in the documentation. There do seem to be just six solutions, from plotting the curve with algcurves.

restart

with(algcurves)

expr := 36*a^3*b^3+8*a^2*b^2*(9*(a+b+1)^2+4*(a*b+a+b))*(a+b+1)+((a-b)^2-2*(a+b)+1)^2*(a+b+1)^5+a*b*((a-b)^2-2*(a+b)+1)*(17*(a+b+1)^2+4*(a*b+a+b))*(a+b+1)^2:

Curve is symmetric - easy to find solutions with a=b

solve(eval(expr, a = b)); evalf(%)

1, 1, 1/2-(1/2)*3^(1/2), 1/2+(1/2)*3^(1/2), 1/2-(1/2)*3^(1/2), 1/2+(1/2)*3^(1/2)

1., 1., -.3660254040, 1.366025404, -.3660254040, 1.366025404

The positive points here points seem to be isolated from the main curve - and there are 6 points in the closed first quadrant.

plot_real_curve(expr, a, b, abserr = 1.*10^(-3), view = [-5 .. 5, -5 .. 5])

 

Download algcurves.mw

@Christian Wolinski I sumitted an SCR.

Use the green up-arrow to upload the worksheet.

@tomleslie @mmcdara Well, @nm asked for "a more elegent way to remove any entry in piecewise which has undefined in it". I didn't really answer that for the general case, but did for the case in question with a Heaviside function. As always, Maple considers the general case and one needs to be careful in special cases. I'm sure @nm will consider that. However, in the context of laplace transforms outputting the unit step function is fine because integrals are involved. In fact taking the Laplace transform of the different outputs gives the same answer:

restart;

with(inttrans):
expr:=exp(-s)/s;
t1:=invlaplace(expr,s,t);
p1:=convert(t1,piecewise);
_EnvUseHeavisideAsUnitStep:=true;
p2:=convert(t1,piecewise);
laplace(t1,t,s);
laplace(p1,t,s);
laplace(p2,t,s);

expr := exp(-s)/s

t1 := Heaviside(t-1)

p1 := piecewise(t < 1, 0, t = 1, undefined, 1 < t, 1)

_EnvUseHeavisideAsUnitStep := true

piecewise(t < 1, 0, 1 <= t, 1)

exp(-s)/s

exp(-s)/s

exp(-s)/s

 

Download Heaviside.mw

For the general case, Maple straddles the two extremes of treating Heaviside as a function or a distribution, and so for other applications I agree the unit step function might not be appropriate.

First 7 8 9 10 11 12 13 Last Page 9 of 45