## dharr

Dr. David Harrington

## 4449 Reputation

18 years, 281 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com 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.

## other methods...

@arashghgood You can try stiff methods other than the default resenbrock one, e.g. lsode. You could put Digits higher for more accuracy, but it will get very slow. In combination with the transformation method to add a small imaginary part you may get closer to what you want. I don't really understand what you are expecting - near a sigularity inaccuracy seems inevitable.

## Stiff method...

@arashghgood Stiff method seems to work here, without transformation. (just plotting real and imm parts of Q)

Edit: at least for the second case. For the first, it varies depending on the exact values of the ic.

## small mod...

@arashghgood  The dchange can be done more simply by

`ode2:=PDEtools:-dchange({K=x+a},ode,[x],params=[a]):`

and still using Q. Result is the same. Hard to understand any discrepancy.

## reformulate...

@arashghgood Perhaps you can reformulate, e.g., K=x+a with x real (including changing the derivatives appropriately). PDEtools:-dchange can do a change of variables - see attached. But ignoring singularities seems risky.

## Z/nZ...

@Thomas Dean Here's my interpretation. The wikipedia page you cited makes clear is the notation for a ring, modulo n. So for the additive group modulo n, the order is n, but for the multiplicative group, 0 is excluded and the order is n-1. So your ord=12, mod 13, is Z/13Z.

This is how you used it in your original post.

## zero and order...

@Thomas Dean Your original table didn't have 0 in the multiplicative case, so I emulated that without thinking too much about it. The Wikipedia page talks about integers 0..n-1 mod n. So that may be part of it. But notice that you are doing ord=12, which is mod 13, which is prime. If you change ord to 11, so mod 12, you get an error message, since the default is to check that it is a group. Try other values, some work and some don't.

## coded elements...

 > > Matrix has to have the entries coded as 1..4

 >   >  >  >  >  > ## FileTools...

I tried to find this out a while ago, but concluded it wasn't possible. You can use FileTools:-ListDirectory to find a list of all *.mw files in the interface(worksheetdir) directory, but there seems not to be a simple way to find which one is running. I suppose there must be some system-dependent system calls to figure out which programs are running.

## vote up...

@mmcdara Much, much simpler than mine!

## algebra...

@gharouce So, to be clear, you never want the actual numbers out, but just to evaluate some special algebra in which e^2=1, e^3=e etc for one specific variable e in an expression? That is easy to do if your expression is a polynomial in e, but requires more work if you want it to work for, say, e/sin(e^2). Here's the polynomial case:

 > > evalpm:=proc(expr,var)   local p,sum,power,term,terms;   p:=collect(expr,var);   if not type(p,polynom(anything,var)) then error "%1 not polynomial in %2",expr,var end if;   if nops(p)=1 then terms:=[p] else terms:=[op(p)] end if;   sum:=0;   for term in terms do     power:=diff(term,var)/term*var;     sum:=sum+ifelse(power::even,coeff(term,var,power),var*coeff(term,var,power));   end do;   sum; end proc:
 >   >  >  > For the more general case

```evalpm:=proc(expr,var)
local p,pwr;
p:=collect(expr,var);
evalindets(p,`^`,
proc(pwr) local res;
if op(1,pwr)=var then
if op(2,pwr)::even then
res:=1
elif op(2,pwr)::odd then
res:=var
else
res:=pwr
end if
else
res:=pwr;
end if;
res;
end proc
);
end proc:```

## eval...

@gharouce You can only set X to one value at a time. But you can do something like this:

```s:=X^3*Y + X^2*Z + X;
eval(s,X=1);
eval(s,X=-1);```

or you can do something like this:

```solve({s=X^3*Y + X^2*Z + X,X^2=1},{s,X});
```

which gives the two possibilities @Ronan I originally had ModuleUnload remove the types, thinking that when I did unwith() the module would unload and the types would be removed. But the types still existed - unwith only gets rid of the global names (exports) the package created. The ModuleUnload procedure "is called when the module is destroyed. A module is destroyed either when it is no longer accessible and is garbage collected, or when Maple exits".

So my conclusion is that it isn't really useful.

## Example...

@Ronan Here's a working example:

 > restart
 > packagedir:=interface(worksheetdir); > MyThings:=module()   option package;   export mysqr,mysqrt;   local ModuleLoad;   uses TT=TypeTools;   ModuleLoad:=proc()    TT:-AddType(mytype,{identical(x),posint});   end proc;   mysqr:=proc(x::mytype);     x^2   end proc;   mysqrt:=proc(x::mytype);     sqrt(2)   end proc;   end module:
 > libfile:=cat(packagedir, "/MyThings.mla"); march( 'create', libfile, 20 ); savelib('MyThings',libfile); >
 > > restart;
 > packagedir:=interface(worksheetdir); > libname:=libname,packagedir:
 > with(MyThings); > mysqr(y);

Error, invalid input: MyThings:-mysqr expects its 1st argument, x, to be of type mytype, but received y

 > mysqrt(2.5);

Error, invalid input: MyThings:-mysqrt expects its 1st argument, x, to be of type mytype, but received 2.5

 > mysqr(x); > ## New version...

Deciding whether or not a given cut is minimal (breaks the graph only into two components) is the slow part for the above algorithm; the test uses ConnectedComponents that needs to perform an algorithm that searches through the graph. In the description of the above algorithm two cases were mentioned where this hard test could be avoided. One was that the fundamental cuts were minimal by construction. The second was for ring sums between two fundamental cuts: if the ring sum was between two edge sets with no common edges, it was not minimal and the hard test was not required.

I speculated that this might be possible to simplify for ringsums of more than two cuts. The following algorithm accumulates the ringsums incrementally, and at each stage tests for lack of common edges.

We need to test all 2^r subsets of the set of r fundamental cuts (except the empty set). Let's code the subsets by a list of integers that index the fundamental sets. For example the index [1,2,4,5] will mean evaluating the ring sum of the first, second, fourth and fifth subsets. The subsets iterator in the combinat package can iterate through all subsets of integers in the list [1,2,3,...,r], thereby successively finding all subsets. More importantly, it does it in an order that means we can use previously evaluated ringsums to find new ones. The order is all subsets with zero elements, then one element, then two elements, etc., and in numerical order within those cases, e.g., for r=3 this is [], , , , [1,2], [1,3], [2,3], [1,2,3]. So as an example (for r>4), when we come to evaluate the ringsum of subsets [1,2,3,7,9] we find the ringsum of subsets [1,2,3,7] has already been done, and so we can evaluate the ringsum between the [1,2,3,7] result and fundamental cut , and store the result in cutsets[1,2,3,7,9] for later use. If the [1,2,3,7] and  edge sets have no common edges, this will not be minimal and we can avoid the hard test in this case.

This was implemented using a Maple table cutsets that was indexed with the lists of integers. All 2^r cutsets need to be stored, so there is a larger memory requirement. Profiling the algorithm on @Icz's graph with 20 vertices  and 72 edges and 314415 minimal cuts showed that the time saved by incrementally doing the ringsums was approximately offset by the extra complexity of the algorithm. For this graph at least, although many more hard tests were avoided, this was still a small fraction of the total.

Old         New
fundamental             19          19
easy (no common edges) 104       25063
hard tests          524164      499205
total = 2^19-1      524287      524287

The total times were appromimately the same. 342 out of the 370 seconds (92%) was spent on hard tests, so those are still limiting here. The number of hard tests was reduced by only 5%.

Edit: This analysis suggests that testing all pairs out of 2^r may be the way to go, or perhaps avoid a brute-force method using the much more complicated algorithm of @Icz's ref .

Update: testing all pairs is much, much slower, so the mixed strategy here or above seems to be optimum for a brute force method.

 > > MinimalEdgeCuts:=proc(G1::GRAPHLN)   local mincutsets,edges,treeedges,G,treeG,     i,j,r,vertexpartition,partitionedges,nf,fc,     fsets,ref,fcutsets,cutsetedges,cutsets,iter,f1,f2,f1intf2,f1ringf2,s;   uses GraphTheory;   if not IsConnected(G1) then error "graph is not connected" end if;   G:=UnderlyingGraph(G1); # make undirected, unweighted, without loops   r:=NumberOfVertices(G)-1;   edges:=Edges(G);   cutsets:=table();     # all 2^r cuts   mincutsets:=table();  # first r ones are the fundamental ones   iter:=combinat:-subsets([\$1..r]); #iterate through all 2^r   iter[nextvalue](); #skip empty list   # choose a spanning tree   treeG:=SpanningTree(G);   treeedges:=Edges(treeG);   # find r corresponding fundamental cutsets, one for each tree edge   for j to r do     iter[nextvalue]();       vertexpartition:=ConnectedComponents(DeleteEdge(treeG,treeedges[j],'inplace'=false));     partitionedges:=`union`(map(Edges,map2(InducedSubgraph,G,vertexpartition))[]);     mincutsets[j]:=edges minus partitionedges;     cutsets[j]:=mincutsets[j];   end do;   j:=r;   # now find ringsums of all subsets of the set of fundamental cutsets   while not iter[finished] do     s:=iter[nextvalue]();     f1:=cutsets[op(1..-2,s)];     f2:=cutsets[s[-1]];     f1intf2:=f1 intersect f2;     f1ringf2:=(f1 union f2) minus f1intf2;     cutsets[op(s)]:=f1ringf2;     if f1intf2 = {} then       next     else       # hard test       vertexpartition:=ConnectedComponents(DeleteEdge(G,f1ringf2,inplace=false));       if nops(vertexpartition)=2 then # cutset is minimal         j:=j+1;         mincutsets[j]:=f1ringf2       end if     end if;   end do;   Vector(j,mincutsets); end proc:
 > with(CodeTools:-Profiling):
 > Profile(MinimalEdgeCuts);
 > g := "S~tIID@OI?{@n~V?goYEDOWd?qI?sJ?[C"; G := GraphTheory:-ConvertGraph(g); #GraphTheory:-DrawGraph(G);  > mincutsets:=CodeTools:-Usage(MinimalEdgeCuts(G)): numelems(mincutsets);

memory used=22.76GiB, alloc change=447.40MiB, cpu time=6.16m, real time=3.88m, gc time=2.83m > PrintProfiles(MinimalEdgeCuts);

MinimalEdgeCuts
MinimalEdgeCuts := proc(G1::GRAPHLN)
local mincutsets, edges, treeedges, G, treeG, i, j, r, vertexpartition,
partitionedges, nf, fc, fsets, ref, fcutsets, cutsetedges, cutsets, iter, f1
, f2, f1intf2, f1ringf2, s;
|Calls Seconds  Words|
PROC |    1 369.781 3054875702|
1 |    1   0.000    520| if not GraphTheory:-IsConnected(G1) then
2 |    0   0.000      0|     error "graph is not connected"
end if;
3 |    1   0.000  90679| G := GraphTheory:-UnderlyingGraph(G1);
4 |    1   0.000     10| r := GraphTheory:-NumberOfVertices(G)-1;
5 |    1   0.000    837| edges := GraphTheory:-Edges(G);
6 |    1   0.000    263| cutsets := table();
7 |    1   0.000    263| mincutsets := table();
8 |    1   0.000   3528| iter := combinat:-subsets([\$ 1 .. r]);
9 |    1   0.000     72| iter[nextvalue]();
10 |    1   0.016   2570| treeG := GraphTheory:-SpanningTree(G);
11 |    1   0.000    192| treeedges := GraphTheory:-Edges(treeG);
12 |    1   0.000      0| for j to r do
13 |   19   0.000   1432|     iter[nextvalue]();
14 |   19   0.000  22390|     vertexpartition := GraphTheory:-
ConnectedComponents(GraphTheory:-DeleteEdge(
treeG,treeedges[j],('inplace') = false));
15 |   19   0.000  24057|     partitionedges := `union`(map(GraphTheory:-
Edges,map2(GraphTheory:-InducedSubgraph,G,
vertexpartition))[]);
16 |   19   0.000   1710|     mincutsets[j] := edges minus partitionedges;
17 |   19   0.000    266|     cutsets[j] := mincutsets[j]
end do;
18 |    1   0.000      0| j := r;
19 |    1   3.621      4| while not iter[finished] do
20 |524268  10.529 58981573|     s := iter[nextvalue]();
21 |524268   1.023 6028911|     f1 := cutsets[op(1 .. -2,s)];
22 |524268   0.548 1048536|     f2 := cutsets[s[-1]];
23 |524268   3.079 9436865|     f1intf2 := f1 intersect f2;
24 |524268   4.557 50878477|     f1ringf2 := f1 union f2 minus f1intf2;
25 |524268   0.875 4414435|     cutsets[op(s)] := f1ringf2;
26 |524268   1.793 1572856|     if f1intf2 = {} then
27 |25063   0.000      0|         next
else
28 |499205 341.972 2917786560|         vertexpartition := GraphTheory:-
ConnectedComponents(GraphTheory:-
DeleteEdge(G,f1ringf2,inplace = false));
29 |499205   1.013 1497615|         if nops(vertexpartition) = 2 then
30 |314396   0.267      0|             j := j+1;
31 |314396   0.473 2762501|             mincutsets[j] := f1ringf2
end if
end if
end do;
32 |    1   0.015 318580| Vector(j,mincutsets)
end proc

 > 