dharr

Dr. David Harrington

4449 Reputation

20 Badges

18 years, 281 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

@vs140580 I already answered your (modified) question about finding the length of the shortest walk that returns to the same vertex. Please look through the answer for the procedure and explanation.

@Carl Love 

restart:

with(GraphTheory):

C1 := {seq({i, i+1}, i=1..6), {7, 1}}:
C2 := {{1, 8}, {8, 3}, {8, 6}}:
G := Graph(C1 union C2):
DrawGraph(G, style=planar);

cb:=CycleBasis(G);

[[1, 2, 3, 8], [1, 2, 3, 4, 5, 6, 7], [1, 7, 6, 8]]

How to find [3, 4, 5, 6, 8]? Considered as sets of edges, the symmetric difference (or "ringsum") operation combines cycles into other cycles  - it kills common edges.

Actually finding cycles or counting all of them is nontrivial. Here combining all three cycles works.

cycles:=map(x->Edges(InducedSubgraph(G,x)),cb);

[{{1, 2}, {1, 8}, {2, 3}, {3, 8}}, {{1, 2}, {1, 7}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}}, {{1, 7}, {1, 8}, {6, 7}, {6, 8}}]

symmdiff(cycles[1],cycles[3],cycles[2]);

{{3, 4}, {3, 8}, {4, 5}, {5, 6}, {6, 8}}

NULL

Download FindCycles.mw

@Carl Love The last one in particular is nice. This operation is one I need to do fairly often and am always frustrated there isn't some simpler (to the average user) way. The help page for solve,details has an example using eval(r, op(indets(r)) = 3) but you have to know it's the third op. My answer has the same problem. Yours allow _Z11 etc, so again you have to know what possibililities might occur. This is compounded by the fact that extra solve commands lead to higher numbers, so the first one isn't always _Z1. To really to do it properly probably some slicing and dicing of the name is necessary.

@Samir Khan Thanks. Most other software I have, the corner handle resizes keeping the aspect ratio the same and handles on the vertical and horizontal edges change one dimension and therefore the aspect ratio.

@Carl Love But this replaces _Z4 as well, which is not what the OP wanted.

@rasmusgs At least for integrals the answer is simple; just add the option numeric, e.g. 

int(sin(x), x = 0 .. Pi, numeric)

gives 2.0000

(But sometimes, as here the exact value is nice, so maybe you want to try without numeric first.)

@vs140580 Here some minimal changes. I made it a row Vector rather than 1x15 Matrix. To see more entries of a Vector or Matrix use for example interface(rtablesize=20) at the beginning of your worksheet.

The last line of the procedure would normally be the return value. (I could have used A rather than return A.) Then you would print outside the procedure; normally print is not used within procedures to provide results.

To convert entries to floats use evalf(A). 

To work with all edges in a single line of code use map or map2. For an example see my other answer here.

Download Try_Degree_based.mw

@jud If I understand your question, then this isn't possible in general. GaloisGroup only gives the permutations between the roots, which is not enough to reconstruct the roots. Here is a simpler example that is my take on it.

restart

with(GroupTheory)

Polynomial and its first root alpha

p := x^4-x^2+1; alias(alpha = RootOf(p, x))

x^4-x^2+1

alpha

The splitting field has only one algebraic number (alpha), the extension field is Q(alpha), and so all the roots can be written directly in terms of alpha.

PolynomialTools:-Split(p, x, 'K'); K

(alpha^3-alpha+x)*(x+alpha)*(-alpha^3+alpha+x)*(x-alpha)

{alpha}

So now the permutations in the Galois group take one root to another. Transitive means for any pair from 1,2,3,4 there is a group operation (permutation) that takes the first of the pair to the second - easily seen to be true.

G := GaloisGroup(p, x); els := Elements(G); IsTransitive(G)

_m1631021308864

{_m1631019905600, _m1631019906208, _m1631021307648, _m1631021308704}

true

So first of all we do not know which of 1,2,3,4 is alpha. We are lucky in this case that rotating through the numbers gives the same result, so it doesn't matter.

map2(PermApply, Perm([[1, 2, 3, 4]]), els)

{_m1631019905600, _m1631019906208, _m1631021307648, _m1631021308704}

So taking 1 as alpha then the permutation (14)(23) takes alpha to the 4th root. But this doesn't help me find that root.

Only after I calculate the roots and look at them do I see that automorphism (14)(23) is complex conjugation.

Specifically, it takes a+b*alpha = a+b*(sqrt(3)+I)/2 to a + b*(sqrt(3)-I)/2, which leaves a unchanged (a,b in Q).

 (Roots 1,2,3 and 4 are in quadrants 1,2,3, and 4).

rts := map(convert, [seq(RootOf(p, x, index = i), i = 1 .. 4)], radical); plots:-complexplot(rts, style = point, symbol = solidcircle, symbolsize = 15, color = red, size = [300, 300], scaling = constrained)

[(1/2)*3^(1/2)+(1/2)*I, -(1/2)*3^(1/2)+(1/2)*I, -(1/2)*3^(1/2)-(1/2)*I, (1/2)*3^(1/2)-(1/2)*I]

Compare with this example

p2 := x^4-10*x^2+1; alias(beta = RootOf(p2, x))

x^4-10*x^2+1

alpha, beta

Different relationship between the roots

PolynomialTools:-Split(p2, x, 'K2'); K2

(beta^3-10*beta+x)*(x+beta)*(-beta^3+10*beta+x)*(x-beta)

{beta}

But the same Galois group.

G2 := GaloisGroup(p2, x); Elements(G2)

_m1631006904768

{_m1631019905600, _m1631019906208, _m1631021307648, _m1631021308704}

(1,4)(2,3) still exchanges roots (1 and 4) and (2 and 3), but it is no longer complex conjugation as an automorphism. (It is not even a reflection in the complex plane.)

rts2 := map(convert, [seq(RootOf(p2, x, index = i), i = 1 .. 4)], radical); plots:-complexplot(rts2, style = point, symbol = solidcircle, symbolsize = 15, color = red, size = [300, 300], scaling = constrained)

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

NULL

Download Galois.mw

@Sphericalmoments A general method would be

1. Find all cycles in the corresponding undirected graph.
2. Find the 0, 1 or 2 directed cycles corresponding to each cycle.
3. Count them.
The DirectedCycles routine I gave earlier does step 2. Step 1 is the difficult step. It is possible to find all cycles from the cycle basis, but the algorithm is nontrivial. So here I give a very inefficient way of doing it, which is probably OK for small graphs.

[Edit: Updated routine works by finding all potential directed cycles on n vertices (rotationally distinct permutations on 1..n (with smallest integer first)) and then sees if the graph contains these edges. Still very inefficient]

restart

with(GraphTheory)

Finds all directed k-cycles in a directed graph.

DirectedCycles:=proc(G1::GRAPHLN,k::posint)
   local i,j,m,s,perm,perms,c,cedges,n,edges,G;
   uses GraphTheory;
   if not IsDirected(G1) then error "graph not directed" end if;
   n:=NumberOfVertices(G1);
   if k<2 or k>n then error "need number of vertices  >= k > 1" end if;
   G:=RelabelVertices(G1,[$(1..n)]);
   edges:=Edges(G);
   # Find all rotationally unique lists of k vertices
   # out of a total of n (smallest integer first).
   j:=0;
   s:=table();
   for i to n-k+1 do # first (smallest)
     perms:=Iterator:-Permute(i+1..n,k-1); # the rest
     for perm in perms do
       c:=[i,seq(perm[m],m=1..k-1),i];
       cedges:={seq([c[m],c[m+1]],m=1..k)};
       # check if this is a directed cycle in the graph
       if cedges subset edges then
         j:=j+1;
         s[j]:=c[1..k];
       end if;
     end do;
   end do;
   convert(s,set);     
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)

Find all directed 2-cycles, 3-cycles and 4-cycles

DirectedCycles(G, 2); nops(%); DirectedCycles(G, 3); nops(%); DirectedCycles(G, 4); nops(%)

{[1, 2], [1, 3], [2, 3]}

3

{[1, 2, 3], [1, 3, 2], [2, 4, 3]}

3

{[1, 2, 4, 3]}

1

G := CycleGraph(3, directed = true); vp := GetVertexPositions(G); G := AddVertex(G, 4); SetVertexPositions(G, [vp[], [0, 0]]); AddArc(G, {[1, 4], [4, 2], [4, 3]}); DrawGraph(G, size = [250, 250])

Number of directed 3-cycles

DirectedCycles(G, 3); nops(%)

{[1, 2, 3], [1, 4, 3]}

2

Number of directed 4-cycles

DirectedCycles(G, 4); nops(%)

{[1, 4, 2, 3]}

1

NULL

Download Cycles4.mw

@Christian Wolinski I've just starting using the Logic package, and saw that in a long expression that I'd Imported and assumed it was an extra variable that was always true, added to put the expression in some sort of standard form. But here is a smaller example, showing it is a bug. I'll submit an SCR.

restart;

with(Logic):

q:=(x[2]+x[8]+x[9])*x[1]:

q2 := Import(q, form = MOD2);

Logic:-`&and`(Logic:-`&xor`(x[2], x[8], x[9]), x[true])

NULL

Download x_true.mw

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

3 4 5 6 7 8 9 Last Page 5 of 42