## dharr

Dr. David Harrington

## 5136 Reputation

19 years, 102 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.

## MaplePrimes Activity

### These are replies submitted by dharr

@emersondiaz I assume you are talking about the dsolve(DE,...) solution. Your worksheet doesn't reproduce the error you mention (I instead get a message about a singulaity) so I don't know how you got that. So here are some comments:

1. The default method automatically adjust the stepsize, so you should not use the stepsize option.

2. If your equations occur with many diffferent timescales (or rscales here), then you may need a stiff solver - add stiff=true (I routinely use this).

3. Terms in the equations and the numbers you are interested in for r and R are many orders of magnitude away from 1 (10^6 and 10^45 appear with 10^(-31). For example (r^2 - 1.690735000*10^37*r + 0.81). So the equations need to be rescaled (non-dimensionalize would be typical, or use a different unit solution) to avoid these numbers.

Point 3 is the most important point.

## with output...

Like this: ## 2-region PDE...

@NIMA112 You want now to solve a PDE on two regions with matching conditions at the common boundary, Eqs 6.30-6.33. As far as I know, Maple's numerical PDE solver only handles single rectangular regions. A Tool like COMSOL is probably better suted to this. Or you could follow the procedure given to find the complicated series solution; I didn't try to understand that in detail.

## not obvious...

@mmcdara Thanks. Of course I first tried the obvious elimination of all of v,a,b,R,T, which doesn't work. v turns into Z to make the cubic, so I tried that first and made progress. Then since A is a nondimensionalized version of a, and similarly for b I added them and the R and T disappeared because they "had to" as they weren't in the final result.

I use solve more than I use eliminate, but here I would have got a RootOf solving for Z: solve({eqA, eqB, eqP, eqZ}, {Z, a, b, v}).

## nice analytical solution method...

@mmcdara You analysis appeared as I was working on the numerical solution. Vote up.

## remembering...

@Carl Love Actually, I wasn't relying on that in the procedure. If it does remember (probably most times), then the loop executes faster, but otherwise it is just slower.

## convert, polynom...

@ider I think @sursumCorda's point is that by using a truncation order of zero you only get the principal part and the O() term, which you can remove by convert/polynom. So the following should give the principal part:

`convert(series(BesselK(4, x), x = 0, 0),polynom);`

You suggested it sometimes stops earlier, which I didn't understand. Do you have an example of that?

(The documentation for laurent suggests some subtlety when there are an infinite number of negative powers, which I didn't follow.)

## VF2++...

@dharr Now the full VF2++ algorithm is implemented. Good for relatively sparse graphs but dense ones still difficult. Perhaps VF3...

Notes: It's currently a memory hog. Information about states that won't be visited could be deleted; T1 etc sets could be incrementally generated as needed. But for dense graphs, the VF3 algorithm may be better, since it (dynamically) uses information about G2 to prune the search tree; VF2++'s preprocessing isonly for G1.

Find if G1 is isomorphic to an induced subgraph of G2, using the VF2++ algorithm.
Graphs must be undirected without self-loops.

See VF2++ paper doi: 10.1016/j.dam.2018.02.018 for notation and description.
Code is in startup region (and is reproduced at the end of this worksheet)

 > restart;
 > with(GraphTheory):
 > claw:=CompleteGraph(1,3): g:=CompleteGraph(2,2,2,2): #DrawGraph(claw,size=[200,200]); #DrawGraph(g,size=[200,200]); IsSubgraphInducedIsomorphic(claw,g);   > IsSubgraphInducedIsomorphic(G1,G2,matching=true); Snake in the box - see Wikipedia

 > G1:=PathGraph(14); # longest for HyperGraph(5); G2:=SpecialGraphs:-HypercubeGraph(5); #DrawGraph(G1); #DrawGraph(G2);  > t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true); IsIsomorphic(G1,InducedSubgraph(G2,map(rhs,m)));  Coil in a box

 > G1:=CycleGraph(14); # longest for HyperGraph(5); G2:=SpecialGraphs:-HypercubeGraph(5); #DrawGraph(G1); #DrawGraph(G2);  > t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true); IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));  But it is much slower to decide there are no matches when the cycle (or path) is one unit longer.

 > CodeTools:-Usage(IsSubgraphInducedIsomorphic(CycleGraph(15),G2,matching=true));

memory used=16.95GiB, alloc change=0.53GiB, cpu time=5.50m, real time=3.23m, gc time=2.75m Random low density graph (length of output message here is a bug - SCR reported)

 > randomize(740472561204): G2 := RandomGraphs:-RandomGraph(10000, 0.1); NumberOfEdges(G2); p:=combinat:-randcomb(Vertices(G2),30): G1 := RelabelVertices(InducedSubgraph(G2,p),[\$1..nops(p)]);   > t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true)); IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

memory used=483.58MiB, alloc change=-114.43MiB, cpu time=3.09s, real time=2.59s, gc time=718.75ms  Random smaller higher density graph is much harder

 > randomize(22807444796519): G2 := RandomGraphs:-RandomGraph(50, 0.5); p:=combinat:-randcomb(Vertices(G2),11): G1 := RelabelVertices(InducedSubgraph(G2,p),[\$1..nops(p)]);  > t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true)); IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

memory used=5.47GiB, alloc change=0 bytes, cpu time=54.72s, real time=47.40s, gc time=8.72s  > Code from startup region reproduced here

 > IsSubgraphInducedIsomorphic := proc(G1::Graph, G2::Graph, {matching::truefalse:=false})         description "Determines if undirected loopless G1 is isomorphic to an induced subgraph of G2",                         "and if so optionally returns the matching function for the vertices";         local n1,n2,D1,W1,V1,A1orig,A1,D2,W2,V2,A2,perm,MM,Morder,r,Vd,V1remaining,                 depth,allV1,allV2,Mout,Mstate,expandmatching,Premaining,                 M,M1,M2,T1,T1sorted,T2,T1t,T2t,P,pair,consistent,cut,i,j;         D1,W1,V1,A1orig := op(1..4,G1); # directedness,weightedness,vertices,neighbors         D2,W2,V2,A2     := op(1..4,G2);         if D1 <> 'undirected' or D2 <> 'undirected' then error "both graphs must be undirected" end if;         n1 := nops(V1):                # number of vertices in G1         n2 := nops(V2):         if orseq(i in A1orig[i], i = 1..n1) or orseq(i in A2[i], i=1..n2) then error "graph(s) have self-loops" end if;         allV1 := {\$1..n1};                # all vertices in G1         allV2 := {\$1..n2};         # simple rejections based on number of vertices and max degrees           if n1 > n2 then return false end if;         if max(map(nops,A1orig)) > max(map(nops,A2)) then return false end if;         # First presort vertices of G1 in special BFS order         # We are not using labels so can presort in decreasing degree order outside outermost while loop         V1remaining := sort([\$1..n1], (a,b) -> nops(A1orig[a]) >= nops(A1orig[b]));         Morder:=Array([]);         MM:={};         while numelems(V1remaining) <> 0 do # loop for each connected component                 r := V1remaining; # root for BFS tree for this component                 Vd := [r]; # depth = 0                 while nops(Vd) > 0 do # for each depth                         # sort according to decreasing connections to M and                         # otherwise decreasing degree                         if nops(Vd) > 1 then                                 Vd:=sort(Vd, proc(a,b)                                                         local na := nops(A1orig[a] intersect MM),                                                                  nb := nops(A1orig[b] intersect MM);                                                         na > nb  or na = nb and nops(A1orig[a]) >= nops(A1orig[b])                                                    end proc)                         end if;                         # add these to the order and update                         Morder ,= Vd[];                         MM := MM union {Vd[]};                         V1remaining := remove(has, V1remaining, Vd);                         # go one level deeper by finding new neighbors of Vd;                         Vd:=convert(`union`(seq(A1orig[Vd])) minus MM, list);                 end do;         end do;                 perm:=[seq(Morder)]; # order of vertices we want.         # We now reorder and relabel the graph so the chosen order becomes the natural order         A1 := A1orig[perm];         A1 := eval(A1, perm =~ [\$1..n1]);         # Start of main part of algorithm         # current matching M is represented as a set of equations         # that map G1 vertices (lhs) to G2 vertices (rhs), e.g., (2=4,3=5}         #         # expandmatching works out M1,M2,T1,T1sorted,T2,T1t,T2t for a given M         expandmatching := proc(M) option remember;                 local                         M1 := map(lhs,M), # set of verts in G1 in matching                         M2 := map(rhs,M), # set of verts in G2 in matching                         T1 := {seq(op(A1[i]), i in M1)} minus M1, # set of neighbors of M1 not in M1                         # sort T1 in order of decreasing numbers of neighbors in M1                         T1sorted := sort([T1[]],(a,b)->nops(A1[a] intersect M1) >= nops(A1[b] intersect M1)),                         T2 := {seq(op(A2[i]), i in M2)} minus M2, # set of neighbors of M2 not in M2                         # T1t = T1tilde are further away vertices                         T1t := allV1 minus M1 minus T1,                         T2t := allV2 minus M2 minus T2;                 M1,M2,T1,T1sorted,T2,T1t,T2t         end proc;         # Premaining is table indexed by a matching that gives         # the remaining addition possibilities (P) not yet tried for this matching         Premaining := table('sparse');         M := {};         Mout := 0;         depth := 0;         # main loop - each iteration assesses a candidate state M         # in a depth-first search         while depth >= 0 do                 M1,M2,T1,T1sorted,T2,T1t,T2t := expandmatching(M);                 if M1 = allV1 then # success                         Mout := M;                         break                 end if;                 # find remaining candidate pairs for addition to this state                 if Premaining[M] = 0 then # new M                         if nops(T1) <> 0 then                                 P := {seq([T1sorted,j],j in T2)};                         elif         nops(allV1 minus M1) > 0 then                                         P := {seq([(allV1 minus M1),j], j in (allV2 minus M2))};                         else                                 P:={}                         end if;                         Premaining[M] := P;                 else                         P := Premaining[M];                 end if;                 # check out possibilities and go deeper if one is found                 while numelems(P) > 0 do                         pair := P;                         P := P minus {pair};                         Premaining[M] := P;                         i,j := pair[];                         # consistency check for [i,j]                         # neighbors of i in M1 must correspond to neighbors of j in M2                         consistent:=evalb(eval(A1[i] intersect M1, M) = (A2[j] intersect M2));                         # cut rule                         cut := evalb(nops(A2[j] intersect T2) < nops(A1[i] intersect T1)                                 or nops(A2[j] intersect T2t) < nops(A1[i] intersect T1t));                         if consistent and not cut then                                 Mstate[depth] := M; # save state                                 depth := depth + 1;                                   M := M union {i = j};                                 next 2                         end if                 end do;                 # no possibilities, backtrack                 depth := depth - 1;                 M := Mstate[depth];         end do;         forget(expandmatching);         # return values         if Mout <> 0 then                    if matching then                         true, map(eq->V1[perm[lhs(eq)]]=V2[rhs(eq)], Mout)                 else                         true                     end if;           else                 false         end if; end proc:
 > ## workbook files...

@mmcdara @jalal The "this:///..." notation refers to a file within the workbook (not a .mw file), where "this" means within the current workbook; see

this help

So I think the original problem was, as @mmcdara pointed out, that it was not a string. But you don't wand to append this to currentdir(). I think the Read needs only the "this///...", but until you upload your workbook it is hard for someone to check this. As @acer pointed out, you need to compress your workbook to a .zip file before it can be uploaded to Mapleprimes - I see you have now done that.

## timimgs...

@sursumCorda Yes, it is very fast. Adding iterations=100 to the CodeTools:-Usage shows they are both about the same time, so the apparently much faster time is a fluctuation. Vectors of integers also can be Trimmed.

## like this...  > > >  >  >  >  ## version...

@mmcdara I agree with the OPs observations of errors in Maple 2023. I found them hard to track down, and got confused by the gammas. But please check, you have both gamma_i and gamma__i in the gap procedure. .I usually define the ode and ic set outside of the procedure to be passed to minimize, and do eval on them within the procedure.

## exact solution...

@mmcdara Your

`solution := tanh((1/4)*sqrt(2)*x)`

is I believe the solution the OP referred to. It is an exact solution for z(0)=0, z(infinity)=1 (as verified by taking the limit), but is only an approximate solution for z(0)=0, z(15)=1.

## confusion...

@mmcdara It is a confusing situation - while there isn't a help page for, say, `[]`, asking for help does take you to a page that descibes the user-level equivalent. And yet there is a help page for `?[]`.

## `<,>` and `<|>`...

@mmcdara I think the OP wanted documentation for their usage in the form `<,>`(1,2,3) and `<|>`(1,2,3).

﻿