Advanced Search

dharr

Dr. David Harrington

8270 Reputation

22 Badges

20 years, 357 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. 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

works in 1-D...

dharr 8270
August 12 2023
1 5

It works in 1-D without the extra parentheses, which are not in your code region versions. 

CP2 := (X, Y) ->local  x,y; {seq(seq([x, y], y = Y), x = X)};

Edit: The code edit region codes and 1-D codes should work identically. 2-D has not kept up with recent additions to the Maple language, so my workaround is just to use 1-D (or in this case leave out the local declaration and ignore the warnings)

vertical plot range...

dharr 8270
August 05 2023
1 0

If you remove the 0.00005..0.00015 from the plot, you see the curves are all close to 1-eta, and there was nothing to see in the vertical range you specified. Of course I can't comment on why the parameters don't give what you were expecting.

no...

dharr 8270
August 04 2023
2 4

There was earlier comment on related issues here and here. The conclusion seems to be high Digits (pehaps 40) are frequently required. If you have a nice matrix (normal or perhaps even just diagonalizable), you can just do the function on the eigenvalues and bypass Maple's routine. But in your case you have a non-diagonalizable matrix (JordonForm not diagonal). Then the general algorithm is much more complicated and involves an interpolating polynomial though the eigenvalues, and I'd guess this leads to numerical issues. In Horn and Johnson's discussion of this (Topics in Matrix Analysis, ch 6), they use Lagrange-Hermite interpolation for their discussion of the analytical case, but that would not be a good way to do it numerically - not sure of the method Maple uses. Both the nature of the function and the nature of the matrix mean this is a difficult problem numerically.

working...

dharr 8270
August 01 2023
3 3

You didn't upload your worksheet (using green up-arrow in the Mapleprimes editor) so it's hard to tell. For me in Maple 2023 it is working. I changed some "." to "*" ( as shown in red), though thay didn't seem to be a problem. The plot looks correct, the Gibbs oscillations always need lots of terms to be less prominent.

Download fourier.mw

Edit - I realised your output is the general formula for any N. Do you really want that? Usually the sum you get when asking for a specific numberof terms, say fourier_f(5), is what is expected. If you don't need the general formula, add rather than sum is better.

> 

restart;
> 

assume(n,posint);
> 

target_f := x -> piecewise(-Pi < x and x < 0, 0, 0 < x and x < Pi, x, Pi < x and x < 2*Pi, 0, 2*Pi < x and x < 3*Pi, x - 2*Pi);
a0 := simplify(int(target_f(x), x = -Pi .. Pi)/(2*Pi));
a_n:=int(target_f(x)*cos(n*x), x = -Pi .. Pi)/Pi;
b_n:=int(target_f(x)*sin(n*x), x = -Pi .. Pi)/Pi;
fourier_f := unapply(a0+'add'(a_n*cos(n*x) + b_n*sin(n*x), n = 1 .. N),N);


 

target_f := proc (x) options operator, arrow; piecewise(-Pi < x and x < 0, 0, 0 < x and x < Pi, x, Pi < x and x < 2*Pi, 0, 2*Pi < x and x < 3*Pi, x-2*Pi) end proc

(1/4)*Pi

(-1+(-1)^n)/(n^2*Pi)

-(-1)^n/n

proc (N) options operator, arrow; (1/4)*Pi+add((-1+(-1)^n)*cos(n*x)/(n^2*Pi)-(-1)^n*sin(n*x)/n, n = 1 .. N) end proc

> 

plot(fourier_f(20),x=-Pi..Pi);

> 

 

Download fourier_with_add.mw

quiz...

dharr 8270
July 29 2023
2 1


 

> 

restart
> 

with(Grading)
> 

Quiz("What is the sum of $B and $C", proc () options operator, arrow; evalb(Quiz:-Get(`$RESPONSE`) = Quiz:-Get(`$B`)+Quiz:-Get(`$C`)) end proc, proc () Quiz:-Set(`$B` = (rand(0 .. 100))(), `$C` = (rand(0 .. 100))()) end proc)
> 

NULL

Download Quiz1.mw

faster...

dharr 8270
July 29 2023
2 2

The BellmanFordAlgorithm can do a whole row of the matrix at a time if you supply a list as the last argument, leading to a much faster result.

[Edit - newer version using topological order added]

> 

restart
> 

randomize(14161147548192)

14161147548192

> 

G := GraphTheory:-RandomGraphs:-RandomNetwork(200, .2, 'acyclic', 'weights' = 0 .. .2)
> 

G__0 := applyop(`-`, -1, G)NULL

GRAPHLN(directed, weighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200], Array(1..200, {(1) = {2, 3, 4, 5, 6, 7, 8, 9, 10, 11}, (2) = {3, 6, 8, 10, 13, 16, 17, 18, 20}, (3) = {4, 6, 7, 8, 10, 13, 17, 18, 20}, (4) = {5, 6, 7, 8, 9, 13, 14, 18, 20}, (5) = {6, 7, 10, 14, 15, 18, 19, 20}, (6) = {7, 9, 10, 14, 15, 16, 17}, (7) = {11, 12, 19, 20}, (8) = {9, 10, 11, 12, 15, 17, 18, 20}, (9) = {11, 14, 15, 16, 18}, (10) = {15, 16, 17, 19}, (11) = {15, 16, 17, 18, 19, 20}, (12) = {13, 15, 18, 19, 21, 24}, (13) = {16, 18, 20, 21, 23, 25, 26, 27, 28}, (14) = {15, 16, 17, 18, 20, 22, 27}, (15) = {16, 19, 20, 22, 24, 26, 27, 28}, (16) = {17, 19, 22, 23, 24, 27}, (17) = {18, 20, 26, 27}, (18) = {21, 23, 24, 25, 26, 28}, (19) = {23, 25}, (20) = {21, 23, 24, 25, 26, 28}, (21) = {22, 24, 25, 26, 29, 31, 33, 34, 35}, (22) = {23, 25, 26, 28, 29, 31}, (23) = {24, 26, 28, 30, 32, 34, 35}, (24) = {25, 28, 29, 34, 35}, (25) = {26, 27, 29, 32, 34}, (26) = {27, 30, 32, 34}, (27) = {28, 30, 31, 32, 35}, (28) = {29, 31}, (29) = {30, 32, 34, 35, 37, 40, 42, 45}, (30) = {31, 34, 35, 36, 43, 45}, (31) = {33, 34, 35, 36, 37, 38, 40, 41, 45}, (32) = {33, 34, 35, 36, 37, 41, 42, 44, 45}, (33) = {34, 35, 36, 37, 38, 39, 40, 41, 42, 43}, (34) = {35, 36, 39, 41, 42, 44}, (35) = {37, 40, 42, 43, 44}, (36) = {37, 38, 39, 40, 43, 44, 45, 46, 47, 48, 49, 51, 52}, (37) = {43, 44, 45, 47, 51}, (38) = {40, 41, 42, 43, 44, 46, 48, 49, 52}, (39) = {41, 42, 44, 46, 51}, (40) = {42, 43, 44, 45, 46, 48, 50, 51, 52}, (41) = {44, 49, 51}, (42) = {43, 45, 46, 47, 48, 49, 50, 51}, (43) = {47, 49, 50, 52}, (44) = {45, 46, 49, 52}, (45) = {46, 48, 49, 52}, (46) = {48, 49, 50, 51}, (47) = {48, 49, 51, 52, 53}, (48) = {51}, (49) = {52}, (50) = {52}, (51) = {52}, (52) = {53}, (53) = {54, 55, 56, 57, 58, 59, 60, 61}, (54) = {55, 56, 57, 58, 60, 61, 64}, (55) = {57, 58, 62, 64}, (56) = {57, 60, 61, 62}, (57) = {58, 61, 63, 64}, (58) = {59, 61, 62}, (59) = {60, 61, 62, 63}, (60) = {61, 63, 64}, (61) = {63, 64}, (62) = {65, 66}, (63) = {65, 66}, (64) = {65, 66, 67}, (65) = {67, 69, 70, 75, 76}, (66) = {68, 71, 72, 74}, (67) = {68, 69, 70, 71, 73, 74, 76}, (68) = {69, 71, 73, 74, 75, 77, 82, 84, 85, 86}, (69) = {72, 73, 75, 78, 79, 80, 82, 83, 87}, (70) = {73, 75, 77, 79, 80, 81, 82, 83, 84, 85, 86, 87}, (71) = {76, 78, 80, 83, 84, 86}, (72) = {74, 76, 77, 78, 79, 80, 81, 87}, (73) = {76, 77, 79, 80, 81, 82, 85, 87}, (74) = {75, 78, 81, 83, 84, 86}, (75) = {76, 77, 79, 82, 83, 84, 86, 87}, (76) = {77, 78, 79, 80, 85, 87}, (77) = {79, 80, 83, 85, 87, 89, 93}, (78) = {80, 83, 84, 85, 90, 92, 93}, (79) = {80, 81, 82, 83, 85, 88}, (80) = {82, 85, 86, 87, 88, 90, 91}, (81) = {83, 85, 87, 88, 93}, (82) = {83, 86, 87, 88, 89, 91, 92, 93}, (83) = {84, 85, 86, 87, 88, 90, 91, 92, 93}, (84) = {85, 86, 87, 88, 90, 91}, (85) = {87, 88, 90, 91, 92, 93}, (86) = {89, 90, 92, 93}, (87) = {88, 89}, (88) = {89, 90, 93, 94}, (89) = {90, 92, 94}, (90) = {94}, (91) = {92, 94}, (92) = {93, 94}, (93) = {94}, (94) = {95}, (95) = {96, 97, 98, 99, 100, 101, 102, 103}, (96) = {97, 101, 104, 105}, (97) = {98, 99, 100, 102, 103, 106}, (98) = {105}, (99) = {105}, (100) = {101, 102, 104, 105, 106}, (101) = {102, 104, 105}, (102) = {106}, (103) = {104}, (104) = {105, 107, 109, 110, 111, 112, 113}, (105) = {107, 108, 109, 110}, (106) = {107, 109, 111, 113}, (107) = {111, 114}, (108) = {109, 110, 114}, (109) = {110, 111, 112}, (110) = {111, 113}, (111) = {114}, (112) = {113}, (113) = {114}, (114) = {115}, (115) = {116, 117}, (116) = {118}, (117) = {118}, (118) = {119, 120, 121, 122, 123}, (119) = {125, 129, 130}, (120) = {122, 123, 125, 126, 128, 129, 130}, (121) = {126, 128, 129}, (122) = {123, 124, 125, 126, 127, 130}, (123) = {125, 126, 127, 128, 130}, (124) = {125, 128, 131, 133}, (125) = {126, 127, 128, 131, 132, 135, 137, 138}, (126) = {127, 129, 130, 132, 138}, (127) = {131, 134, 136}, (128) = {130, 131, 132, 137}, (129) = {130, 131, 132, 133, 134, 136, 137}, (130) = {131, 135, 136, 137}, (131) = {133, 134, 136, 137, 139}, (132) = {133, 134, 135, 136, 139, 140, 142, 143}, (133) = {136, 137, 138, 140, 141, 142, 143}, (134) = {135, 136, 139, 140, 141, 143}, (135) = {136, 137, 138, 139, 140, 141, 142, 143}, (136) = {139, 142}, (137) = {140, 141, 142, 143}, (138) = {139, 140, 141, 142}, (139) = {140, 143, 144}, (140) = {141, 142, 144}, (141) = {144}, (142) = {144}, (143) = {144}, (144) = {145, 146, 147, 148, 149, 150, 151}, (145) = {149, 150, 151, 152}, (146) = {149}, (147) = {150, 151, 152}, (148) = {149, 150, 151}, (149) = {152}, (150) = {152}, (151) = {152}, (152) = {153, 154}, (153) = {154, 155, 156}, (154) = {156, 157}, (155) = {156, 157, 159, 160, 161, 162, 163, 164}, (156) = {158, 160, 161, 162, 163, 164}, (157) = {158, 160, 161, 163}, (158) = {161, 162, 164, 166, 167, 170, 172}, (159) = {160, 161, 162, 163, 164, 167, 170, 171, 172}, (160) = {162, 163, 165, 166, 168, 169, 170, 171}, (161) = {163, 164, 168, 169, 170, 171, 173}, (162) = {163, 166, 167, 169, 170, 171, 172, 173}, (163) = {164, 165, 167, 168}, (164) = {165, 167, 169, 170}, (165) = {167, 168, 169, 173, 174, 175, 176, 181}, (166) = {167, 170, 175, 177, 178, 179, 181, 182, 183, 184}, (167) = {169, 170, 176, 177, 180, 182, 183, 184}, (168) = {169, 170, 171, 173, 174, 176, 178, 179, 180, 181}, (169) = {171, 174, 175, 176, 178, 179, 181, 183}, (170) = {174, 175, 177, 179, 180, 182}, (171) = {172, 175, 176, 180, 181, 182}, (172) = {173, 174, 175, 176, 184}, (173) = {175, 178, 180, 181, 183}, (174) = {175, 176, 178, 181, 182, 183, 184, 185, 186, 187, 188, 191}, (175) = {176, 177, 181, 182, 183, 184, 185, 189}, (176) = {177, 179, 180, 181, 182, 183, 185, 188, 189, 190, 191}, (177) = {179, 180, 182, 183, 185, 187, 188, 189, 191}, (178) = {180, 182, 183, 184, 185, 186, 188, 189, 190}, (179) = {180, 181, 182, 185, 187, 188, 189}, (180) = {184, 185, 187, 188, 189, 190}, (181) = {183, 184, 186, 188, 189, 191}, (182) = {186, 188, 189, 190, 191}, (183) = {186, 187, 190}, (184) = {185, 188, 189, 191}, (185) = {186, 190, 191, 192, 193}, (186) = {188, 189, 190, 191, 192}, (187) = {189, 192, 193}, (188) = {190, 192, 193}, (189) = {191, 192, 193}, (190) = {191, 193}, (191) = {192}, (192) = {193}, (193) = {194}, (194) = {195, 196}, (195) = {196, 197, 198}, (196) = {197, 198}, (197) = {198, 199}, (198) = {199}, (199) = {200}, (200) = {}}), `GRAPHLN/table/1`, )

> 

t, s := combinat:-randcomb(GraphTheory:-Vertices(G__0), 5^2), combinat:-randcomb(GraphTheory:-Vertices(G__0), integermul2exp(5, 2))

[4, 21, 25, 26, 30, 35, 47, 48, 60, 74, 84, 93, 100, 104, 110, 126, 140, 144, 158, 165, 167, 171, 179, 182, 196], [5, 40, 42, 43, 54, 74, 80, 81, 86, 89, 93, 96, 103, 104, 123, 125, 126, 133, 144, 160]

Original version

> 

"DataFrame((`M__1`:=CodeTools:-Usage(Matrix(numelems(s),numelems(t),(i,j)->numelems((GraphTheory:-BellmanFordAlgorithm(`G__0`,s[i],t[j]))[1]),datatype=integer[2]))),'columns'=t,'rows'=s):"

memory used=8.45GiB, alloc change=-2.00MiB, cpu time=4.58m, real time=4.85m, gc time=28.11s

New version BellmanFordAlgorithm doing one row at a time.

> 

"DataFrame((`M__2`:=CodeTools:-Usage(Matrix(numelems(s), numelems(t),[seq(map(x->nops(x[1]), GraphTheory:-BellmanFordAlgorithm(`G__0`,s[i],t) ),i=1..numelems(s))]))),'columns'=t,'rows'=s):"

memory used=348.13MiB, alloc change=0 bytes, cpu time=9.42s, real time=9.42s, gc time=406.25ms

Original version sped up with topological sort

> 

"DataFrame((`M__3`:=CodeTools:-Usage(Matrix(numelems(s),numelems(t),proc(i::posint,j::posint,` $`)::nonnegint; uses ListTools,GraphTheory; local ts::list(posint):=TopologicSort(`G__0`,'output'='permutation'),q::posint:=Search(t[j],ts),p::posint:=Search(s[i],ts); if p>q then 0 elif q=p then 1 else numelems(BellmanFordAlgorithm(`G__0`,s[i],t[j])[1]) fi end,datatype=integer))),':-columns'=t,':-rows'=s):"

memory used=4.75GiB, alloc change=0 bytes, cpu time=2.29m, real time=2.21m, gc time=10.86s

New version sped up with topological sort

> 

ShortestPathsDAG:=proc(G::Graph,s::list,t::list)
  local ts,tge,lthsge,ns,nt,lths,i,j,M;
  ns:=nops(s);
  nt:=nops(t);
  M:=Matrix(ns,nt);
  ts:=GraphTheory:-TopologicSort(G,'output'='permutation');
  for i to ns do  # row i in Matrix
    # only do Bellman Ford on values in t that are later in topo seq than vertex s[i] for this row
    j:=ListTools:-Search(s[i],ts);
    tge:=convert({t[]} intersect {ts[j..-1][]},list);
    lthsge:=map(x->nops(x[1]),GraphTheory:-BellmanFordAlgorithm(G,s[i],tge));
    lths:=table('sparse',tge=~lthsge);
    M[i]:=Vector[row](nt,j->lths[t[j]]);
  end do;
  M;
end proc:
> 

"DataFrame((`M__4`:=CodeTools:-Usage(ShortestPathsDAG(`G__0`,s,t))),'columns'=t,'rows'=s):"

memory used=348.33MiB, alloc change=0 bytes, cpu time=9.38s, real time=9.06s, gc time=718.75ms

> 

andmap(LinearAlgebra:-Equal, [M__2, M__3, M__4], M__1)

true

> 

NULL

Download longest_paths_in_a_DAG.mw

 

some things...

dharr 8270
July 29 2023
1 1

If you upload your worksheet (green up arrow in the Mapleprimes editor) then we can run the code without retyping it in.

In your first procedure, the wiggly line on line 2 tells you there is a syntax error there, which is that it is seeing you want m as a local even though you used it already in the first line. Your first line should be U::Matrix and not U::Matrix(m).

In the second procedure A:=Matrix is not correct and should be removed. In both cases, the Dimension command is in the LinearAlgebra package and so should be LinearAlgebra:-Dimension (and since it returns two things the second version uses it correctly)

You should assign the procedure to a name, for example fillme:=proc(A)

 

something...

dharr 8270
July 22 2023
2 1

Here's a small step in the right direction.

> 

restart;
> 

eq:=erf(x)=erf(Pi);

erf(x) = erf(Pi)

Define inverse function

> 

inverf:=erf@@(-1);

erf@@(-1)

Apply it to both sides of the equation

> 

map(inverf,eq);

x = Pi

> 

solve(inverf(erf(x))=inverf(erf(Pi)),x);

Pi

But inverf isn't a linear function, so we don't expect anything more complicated to work.

> 

inverf(erf(x)-erf(Pi))=inverf(0);

(erf@@(-1))(erf(x)-erf(Pi)) = (erf@@(-1))(0)

> 

NULL

Download inverf.mw

Can be positive semidefinite...

dharr 8270
July 20 2023
3 1

The following parameters make it positive semidefinite:

My procedure was pretty ad-hoc, using Satisfy and just adding conditions until it worked.

positivedefinite.mw

@mmcdara - I think you needed C[1]<=0, but in 2015 it didn't lead to a solution even when I made that change, 2023 just gave one solution when solving the conditions. Since the matrix is symmetric, the discriminants of the quadratics have to be positive anyway, and the signs of the coefficents can be used to test for two positive roots. But I don't really see why your method gave a different result.


@sursumCorda I'd be suspicious of the deprecated linalg routine. All principal minors (of all sizes) need to non-neg, so perhaps systematically working through these would be a better method. 

datasetlabels = contents...

dharr 8270
July 03 2023
1 0

Use datasetlabels = contents for this, but it isn't pretty.

Functorial composition in combstruct...

dharr 8270
July 02 2023
0 0

The short answer is that there is no obvious way to do this. You can certainly make subsets and two-element subsets as you have described if you are using labeled structures (and exponential generating functions). But the OEIS series you gave is for unlabelled graphs, so I don't see how it could arise from labelled structures. The combstruct package takes a grammar and then enumerates for either labeled or unlabelled structures but not for a mixture.

I don't really understand the functorial composition. But as far as I can see, there doesn't seem to be a simple way that its generating function is related to the generating function of the operands, and the basis of the combstruct package is that the grammar translates directly to operations on generating functions.

(If your object is just to calculate the numbers in the OEIS series, that can be done through Polya counting and can be programmed in Maple.)

improved labels...

dharr 8270
June 28 2023
1 0

For consistency in the A and nA labels, use MathML for both, so set A:=`#mo(A)` (or use mi for both A and nA if you want italics). I also added padding=10 to StyleVertex so the lines don't touch the labels.

For the weight locations you could use GetVertexPositions and then programmatically adjust the locations for the colored edges to have the same y value as the corresponding uncolored ones, and then save with SetVertexPostions. But that is a lot of work, and for me the weights look OK already.

Arbrepondéré.mw

 

dchange...

dharr 8270
June 24 2023
3 15

Perhaps could be done in fewer steps, but this works.

Note: I scaled the result arbitrarily to make it simpler, under the assumption that it is was intended to be equal to zero, but you can leave that out.

> 

restart:
> 

local gamma;

gamma

> 

eqn := (diff(theta(x, z, t), x))^2*(K[1]-K[3])*cos(theta(x, z, t))*sin(theta(x, z, t))+(diff(theta(x, z, t), x))*((diff(theta(x, z, t), z))*(-K[1]*cos(2*theta(x, z, t))+K[3]*cos(2*theta(x, z, t)))-(1/2)*gamma[1]*(4*sin(theta(x, z, t))^2*u(x, z, t)+2*u(x, z, t)*cos(2*theta(x, z, t))))+(diff(theta(x, z, t), z))^2*(K[3]-K[1])*cos(theta(x, z, t))*sin(theta(x, z, t))-(1/2)*gamma[1]*(diff(theta(x, z, t), z))*(4*sin(theta(x, z, t))^2*v(x, z, t)+2*v(x, z, t)*cos(2*theta(x, z, t)))+(diff(theta(x, z, t), z, x))*(-2*K[1]+2*K[3])*cos(theta(x, z, t))*sin(theta(x, z, t))-(diff(u(x, z, t), z))*((1/2)*gamma[2]*cos(2*theta(x, z, t))+(1/2)*gamma[1]*(2*sin(theta(x, z, t))^2+cos(2*theta(x, z, t))))-(diff(v(x, z, t), x))*((1/2)*gamma[2]*cos(2*theta(x, z, t))+(1/2)*gamma[1]*(-2*sin(theta(x, z, t))^2-cos(2*theta(x, z, t))))-(1/2)*gamma[1]*(4*sin(theta(x, z, t))^2*(diff(theta(x, z, t), t))+2*(diff(theta(x, z, t), t))*cos(2*theta(x, z, t)))+((diff(u(x, z, t), x))*gamma[2]-(diff(v(x, z, t), z))*gamma[2])*cos(theta(x, z, t))*sin(theta(x, z, t))+f[2](theta(x, z, t))*(diff(theta(x, z, t), x, x))+f[1](theta(x, z, t))*(diff(theta(x, z, t), z, z));
 

(diff(theta(x, z, t), x))^2*(K[1]-K[3])*cos(theta(x, z, t))*sin(theta(x, z, t))+(diff(theta(x, z, t), x))*((diff(theta(x, z, t), z))*(-K[1]*cos(2*theta(x, z, t))+K[3]*cos(2*theta(x, z, t)))-(1/2)*gamma[1]*(4*sin(theta(x, z, t))^2*u(x, z, t)+2*u(x, z, t)*cos(2*theta(x, z, t))))+(diff(theta(x, z, t), z))^2*(K[3]-K[1])*cos(theta(x, z, t))*sin(theta(x, z, t))-(1/2)*gamma[1]*(diff(theta(x, z, t), z))*(4*sin(theta(x, z, t))^2*v(x, z, t)+2*v(x, z, t)*cos(2*theta(x, z, t)))+(diff(diff(theta(x, z, t), x), z))*(-2*K[1]+2*K[3])*cos(theta(x, z, t))*sin(theta(x, z, t))-(diff(u(x, z, t), z))*((1/2)*gamma[2]*cos(2*theta(x, z, t))+(1/2)*gamma[1]*(2*sin(theta(x, z, t))^2+cos(2*theta(x, z, t))))-(diff(v(x, z, t), x))*((1/2)*gamma[2]*cos(2*theta(x, z, t))+(1/2)*gamma[1]*(-2*sin(theta(x, z, t))^2-cos(2*theta(x, z, t))))-(1/2)*gamma[1]*(4*sin(theta(x, z, t))^2*(diff(theta(x, z, t), t))+2*(diff(theta(x, z, t), t))*cos(2*theta(x, z, t)))+((diff(u(x, z, t), x))*gamma[2]-(diff(v(x, z, t), z))*gamma[2])*cos(theta(x, z, t))*sin(theta(x, z, t))+f[2](theta(x, z, t))*(diff(diff(theta(x, z, t), x), x))+f[1](theta(x, z, t))*(diff(diff(theta(x, z, t), z), z))

Change the independent variables first

> 

tr1:={t = T*tau,x = X*h, z = Z*h};
eqn2:=PDEtools:-dchange(tr1, eqn*T*h^2, [tau,X,Z],params=[T,h],simplify);

{t = T*tau, x = X*h, z = Z*h}

-2*T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*(K[1]-K[3])*(diff(diff(theta(tau, X, Z), X), Z))+f[2](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), X), X))*T+f[1](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), Z), Z))*T+T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*(K[1]-K[3])*(diff(theta(tau, X, Z), X))^2-(2*(cos(theta(tau, X, Z))^2-1/2)*(K[1]-K[3])*(diff(theta(tau, X, Z), Z))+u(tau, X, Z)*h*gamma[1])*T*(diff(theta(tau, X, Z), X))-T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*(K[1]-K[3])*(diff(theta(tau, X, Z), Z))^2-(diff(theta(tau, X, Z), Z))*v(tau, X, Z)*T*h*gamma[1]+(-(cos(theta(tau, X, Z))^2*gamma[2]+(1/2)*gamma[1]-(1/2)*gamma[2])*T*(diff(u(tau, X, Z), Z))-T*(cos(theta(tau, X, Z))^2*gamma[2]-(1/2)*gamma[1]-(1/2)*gamma[2])*(diff(v(tau, X, Z), X))-(diff(theta(tau, X, Z), tau))*h*gamma[1]+T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*gamma[2]*(diff(u(tau, X, Z), X)-(diff(v(tau, X, Z), Z))))*h

Now the dependent variables

> 

tr2:={u(tau,X,Z) = xi*h^2*U(tau,X,Z)/alpha[4], v(tau,X,Z) = xi*h^2*V(tau,X,Z)/alpha[4]};
eqn3:=PDEtools:-dchange(tr2, eqn2*2*alpha[4], [U,V],params=[xi,h,alpha[4]],simplify);

{u(tau, X, Z) = xi*h^2*U(tau, X, Z)/alpha[4], v(tau, X, Z) = xi*h^2*V(tau, X, Z)/alpha[4]}

-4*T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*alpha[4]*(K[1]-K[3])*(diff(diff(theta(tau, X, Z), X), Z))+2*f[2](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), X), X))*T*alpha[4]+2*f[1](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), Z), Z))*T*alpha[4]+2*T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*alpha[4]*(K[1]-K[3])*(diff(theta(tau, X, Z), X))^2-2*T*(2*(cos(theta(tau, X, Z))^2-1/2)*alpha[4]*(K[1]-K[3])*(diff(theta(tau, X, Z), Z))+xi*h^3*U(tau, X, Z)*gamma[1])*(diff(theta(tau, X, Z), X))-2*T*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*alpha[4]*(K[1]-K[3])*(diff(theta(tau, X, Z), Z))^2-2*(diff(theta(tau, X, Z), Z))*xi*h^3*V(tau, X, Z)*T*gamma[1]-(T*h*xi*(2*cos(theta(tau, X, Z))^2*gamma[2]+gamma[1]-gamma[2])*(diff(U(tau, X, Z), Z))-T*h*xi*(-2*cos(theta(tau, X, Z))^2*gamma[2]+gamma[1]+gamma[2])*(diff(V(tau, X, Z), X))+2*gamma[1]*(diff(theta(tau, X, Z), tau))*alpha[4]-2*T*gamma[2]*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*h*xi*(diff(U(tau, X, Z), X)-(diff(V(tau, X, Z), Z))))*h^2

Now the rest, assuming f[1](theta(x, z, t)) was really intended as a function and not as multiplying by f[1] (same for f[2])

> 

tr3 := {K[3] = K[1]*k[3], f[1](theta(x, z, t)) = K[1]*F[1](theta(x, z, t)), f[2](theta(x, z, t)) = K[1]*F[2](theta(x, z, t)), gamma[1] = mu*Gamma[1], gamma[2] = mu*Gamma[2]};

eqn4:=PDEtools:-dchange(tr3, eqn3, [k[3], F[1], F[2], Gamma[1], Gamma[2]],params=[K[1],mu],simplify);

{K[3] = K[1]*k[3], gamma[1] = mu*Gamma[1], gamma[2] = mu*Gamma[2], f[1](theta(x, z, t)) = K[1]*F[1](theta(x, z, t)), f[2](theta(x, z, t)) = K[1]*F[2](theta(x, z, t))}

4*T*alpha[4]*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*K[1]*(-1+k[3])*(diff(diff(theta(tau, X, Z), X), Z))+2*K[1]*F[2](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), X), X))*T*alpha[4]+2*K[1]*F[1](theta(tau, X, Z))*(diff(diff(theta(tau, X, Z), Z), Z))*T*alpha[4]-2*T*alpha[4]*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*K[1]*(-1+k[3])*(diff(theta(tau, X, Z), X))^2-2*(-2*(cos(theta(tau, X, Z))^2-1/2)*alpha[4]*K[1]*(-1+k[3])*(diff(theta(tau, X, Z), Z))+xi*h^3*U(tau, X, Z)*mu*Gamma[1])*T*(diff(theta(tau, X, Z), X))+2*T*alpha[4]*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*K[1]*(-1+k[3])*(diff(theta(tau, X, Z), Z))^2-2*(diff(theta(tau, X, Z), Z))*xi*h^3*V(tau, X, Z)*T*mu*Gamma[1]-mu*(T*h*xi*(2*cos(theta(tau, X, Z))^2*Gamma[2]+Gamma[1]-Gamma[2])*(diff(U(tau, X, Z), Z))-T*h*xi*(-2*cos(theta(tau, X, Z))^2*Gamma[2]+Gamma[1]+Gamma[2])*(diff(V(tau, X, Z), X))+2*Gamma[1]*(diff(theta(tau, X, Z), tau))*alpha[4]-2*T*Gamma[2]*cos(theta(tau, X, Z))*sin(theta(tau, X, Z))*h*xi*(diff(U(tau, X, Z), X)-(diff(V(tau, X, Z), Z))))*h^2

> 

NULL

Download dchange.mw

integers but not very random...

dharr 8270
June 22 2023
1 2

Well, all but 2 are random...

> 

restart;
> 

with(LinearAlgebra):

Choose determinant required and Matrix size:

> 

det:=9;n:=4;

9

4

Try to adjust 1,1 and 1,2 entries until it works. Values can be too high to look random - could play with the q values

> 

imax:=10:
for i to imax do
  R:=RandomMatrix(n,n,generator=-5..5):
  R[1,1]:=a:
  R[1,2]:=b:
  eq:=isolve(Determinant(R)=det,q);
until eq<>NULL:
if i>imax then
  print(`too hard`);
else
  print(i);
  R:=eval(R,eval(eq,q=0));
  R:=eval(R,indets(R,name)=~0); # if there is an a or b left, set = 0
  Determinant(%);
end if;

1

Matrix(%id = 36893490371187920396)

9

> 

NULL

Download GivenDet.mw

Zeon algebra...

dharr 8270
June 21 2023
1 0

See also the thread here
How-To-Obtain-All-Cycles--In-A-Graph

> 

restart
> 

with(GraphTheory); with(Bits)

Using the adjacency matrix and zeons to count cycles.

See R. Schott and G.S. Staples, Complexity of counting cycles using zeons, Comp. Math. Appl. 62 (2011) 1828–1837, doi: 10.1016/j.camwa.2011.06.026

Zeons implemented with bitwise operations on Maple integers in the Bits package.

> 

G := AddVertex(PathGraph(5, directed), [6, 7]); AddArc(G, {[1, 3], [3, 7], [4, 6], [6, 7], [7, 1]}); DrawGraph(G, layout = circle, size = [250, 250])

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

> 

n := NumberOfVertices(G)

7

Setting a bit in an integer.

> 

SetBit := proc (p) Bits:-Join(convert(Vector(n, {p = 1}, fill = 0), list)) end proc; SetBits := proc () options operator, arrow; foldl(Bits:-Or, seq(SetBit(i), `in`(i, args))) end proc

Warning, (in SetBits) `i` is implicitly declared local

> 

SetBit(3); String(%); q := SetBits(2, 3); String(%)

4

"001"

6

"011"

Going back

> 

Decode := proc (p) options operator, arrow; zip(proc (x, y) options operator, arrow; if x = 1 then y end if end proc, Bits:-Split(p), [`$`(1 .. n)]) end proc; Decode(q)

[2, 3]

The algorithm starts with two versions of the adjacency matrix. Matrix B has the (i,j) entry as a coded version of i and j that is the beginning of a path or cycle.

Matrix A entry (i,j) has only the coded versionsof j; this is the vertex that will be added to the path at each step by evaluating B.A with a special rule for multiplication.

At any stage the list in entry (i,j) of B contains integers representing possible paths between i and j of a certain length k.

> 

Adj := AdjacencyMatrix(G); A := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then SetBit(j) else 0 end if end proc); B := Matrix(n, n, proc (i, j) options operator, arrow; if Adj[i, j] = 1 then [Or(SetBit(i), SetBit(j))] else [] end if end proc); B, A

Matrix(%id = 36893490887404493628), Matrix(%id = 36893490887404491348)

Consider a list of two paths of length 2 from vertex 1 to vertex 3: 1-2-3 and 1-7-3 that would be in entry (1,3). Part of the next step would be trying to add vertex 7 to each of these paths.

> 

b := [SetBits(1, 2, 3), SetBits(1, 7, 3)]; a := SetBit(7)

[7, 69]

64

The multiplication of two entries in the process of Matrix multiplication adds the a vertex to each of the possibilities in b. If it is a repeated vertex, then the possibility disappears. A duplicate is detected as the And of two values being nonzero. Adding the vertex is done by Or

> 

And(7, 64); And(69, 64); Or(7, 64); Decode(%)

0

64

71

[1, 2, 3, 7]

Note that although we add the vertex, we lose information about the order. So this algorithm cannot output the cycles themselves. The following implements this multiplication.

> 

map((x,y)-> if Bits:-And(x,y)=0 then Bits:-Or(x,y) else NULL end if,b,a);

[71]

A repeatng vertex isn't necessarily a cycle; only if the starting and ending vertex are the same. This will happen if we are multplying B[i,j]*A[j,p] with i=p, i.e., we are calculating a diagonal entry.

To prevent double counting, we only keep the cycle if the smallest vertex is i. The full multiplication rule is

> 

mult:=(x,y)->if y = 0 or x=[] then []
             else
               map((u,w)-> if Bits:-And(u,w)=0 then
                              return Bits:-Or(u,w)
                           else
                              if i=p and i=Bits:-FirstNonzeroBit(u) + 1 then record_cycle end if;
                              return NULL;
                           end if,
                    x,y)
             end if:

Implement the Matrix multiplication with the above routine

> 

MM:=proc(B,A) local C,i,j,p;
C:=Matrix(n,n):
for i to n do # rows of B diving down
  for p to n do # cols of A     
     C[i,p]:=[seq(mult(B[i,j],A[j,p])[],j=1..n)];  #Row(B,i).Column(A,p)
  end do
end do:
C
end proc:
> 

C := MM(B, A)
> 

Show := proc (M) options operator, arrow; map(proc (x) options operator, arrow; if 0 < nops(x) then map(Decode, x) else x end if end proc, M) end proc; Show(C)

Matrix(%id = 36893490887341852180)

Next one

> 

B := C; C := MM(B, A); Show(C)

Matrix(%id = 36893490887341836396)

And again

> 

B := C; C := MM(B, A); Show(C)

Matrix(%id = 36893490887253294844)

So here is a procedure to count cycles in a graph or digraph.

Output is a sparse table, indexed by the number of vertices in the cycles

Use add(eval(CycleCounts(G))) to find the total number of cycles.
This includes 2-cycles, one for each edge, so subtract the number of edges if you do not want 2-cycles in the total.

For undirected graphs, some authors count each cycle twice, one in each orientation, but here each cycle is counted only once.

> 

CycleCounts:=proc(G::Graph)
  uses GraphTheory;
  local k,A,B,C,p,n,i,j,mult,cycles,SetBit;
  cycles:=table('sparse');
  n:=NumberOfVertices(G);
  SetBit := m-> Bits:-Join(convert(Vector(n, {m = 1}, 'fill' = 0), list));
  mult:=(x,y)->if y = 0 or x = [] then []
             else
               map((u,w)-> if Bits:-And(u,w)=0 then
                              return Bits:-Or(u,w)
                           else
                              if i=p and i=Bits:-FirstNonzeroBit(u) + 1 then cycles[k]++ end if;
                              return NULL;
                           end if,
                    x,y)
             end if;
  A:=Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then SetBit(j) else 0 end if);
  B:= Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then {Bits:-Or(SetBit(i),SetBit(j))} else [] end if);
  for k from 2 to n do
    C:=Matrix(n,n):
    for i to n do # rows of B diving down
      for p to n do # cols of A     
        C[i,p]:=[seq(mult(B[i,j],A[j,p])[],j=1..n)];  #Row(B,i).Column(A,p)
      end do
    end do;
    B:=C;
  end do;
  if not IsDirected(G) then # correct for double counting except for 2-cycles
        cycles:=map(i->i/2,cycles);
        if assigned(cycles[2]) then
            if cycles[2]=0 then # edgeless graphs
               evaln(cycles[2])
            else
               cycles[2]:=cycles[2]*2
            end if
        end if
  end if;
  eval(cycles);
end proc:
 
> 

cycles := CycleCounts(G)
> 

"for i,cyc in eval(cycles) do i,cyc end do;add(eval(cycles));"

3, 1

4, 1

5, 1

6, 1

4

> 

"G3:=CompleteGraph(5); cycles:=CycleCounts(G3): for i,cyc in eval(cycles) do i,cyc end do;add(eval(cycles)); "

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

2, 10

3, 10

4, 15

5, 12

47

> 

G2 := RandomGraphs:-RandomGraph(15, .3)

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], Array(1..15, {(1) = {3, 8, 12}, (2) = {5, 7}, (3) = {1, 6, 13, 14}, (4) = {7, 11}, (5) = {2, 8, 14}, (6) = {3, 9, 15}, (7) = {2, 4, 8, 10, 12}, (8) = {1, 5, 7, 10, 11, 12, 13, 15}, (9) = {6, 14, 15}, (10) = {7, 8, 11, 15}, (11) = {4, 8, 10, 12, 15}, (12) = {1, 7, 8, 11, 14}, (13) = {3, 8, 15}, (14) = {3, 5, 9, 12, 15}, (15) = {6, 8, 9, 10, 11, 13, 14}}), `GRAPHLN/table/6`, 0)

> 

cycles := CodeTools:-Usage(CycleCounts(G2))

memory used=138.66MiB, alloc change=0 bytes, cpu time=3.23s, real time=3.36s, gc time=156.25ms

> 

"for i,cyc in eval(cycles) do i,cyc end do;add(eval(cycles)); "

2, 31

3, 11

4, 29

5, 58

6, 112

7, 215

9, 490

8, 349

11, 566

10, 579

13, 254

12, 450

15, 8

14, 78

3230

> 

NULL

Download ZeonsCycleCounts.mw

First 28 29 30 31 32 33 34 Last Page 30 of 82