Advanced Search

dharr

Dr. David Harrington

8482 Reputation

22 Badges

21 years, 33 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

one way...

dharr 8482
March 05 2023
1 0 
mul(map(mul,ListTools:-Collect(A)))

one way...

dharr 8482
March 05 2023
0 0

CN and N1 are 

CN := (G, u) -> add(map2(GraphTheory:-Degree, G, GraphTheory:-Neighborhood(G, u, 'closed')));
N1 := G -> add(map(e->(CN(G,e[1])+CN(G,e[2]))/2, convert(GraphTheory:-Edges(G),list)));

The others are straightforward modifications.

step by step...

dharr 8482
March 02 2023
3 1

This gives the answer in term of four algebraic numbers, each of which can be written in terms of the first root, which I think is more in the spirit of Galois theory than radicals. I also show the minimal polynomials that show the Galois group is order 72, But I'm not using the Galois group for much of this, as the OP asked for. Sorry if this is already obvious to the OP; perhaps we can have some clarification if the answers here are off track.

> 

restart
> 

_local(gamma); poly := x^6-3*x+3

gamma

x^6-3*x+3

Doesn't factor in Q

> 

irreduc(poly)

true

Let alpha be the first root (index=1 by default)

> 

alias(alpha = RootOf(poly, x))

alpha

alpha is the only root that can be written in terms of alpha. Notice that the rest factors, i.e. the quintic is reducible in  Q(alpha)

> 

p6alpha := evala(Factor(poly, {alpha}))

(x^3+((6/17)*alpha-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+15/17)*x^2+((6/17)*alpha^2-33/17+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3-(3/17)*alpha)*x+(5/17)*alpha^3+18/17+(7/17)*alpha^4+(6/17)*alpha^2+(3/17)*alpha^5-(3/17)*alpha)*((11/17)*alpha^5+(3/17)*alpha^4+(7/17)*alpha^3+(5/17)*alpha^2+(6/17)*alpha-36/17+((6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17)*x+x^2)*(x-alpha)

Collect the factors

> 

p3alpha, p2alpha, p1alpha := op(p6alpha)

x^3+((6/17)*alpha-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+15/17)*x^2+((6/17)*alpha^2-33/17+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3-(3/17)*alpha)*x+(5/17)*alpha^3+18/17+(7/17)*alpha^4+(6/17)*alpha^2+(3/17)*alpha^5-(3/17)*alpha, (11/17)*alpha^5+(3/17)*alpha^4+(7/17)*alpha^3+(5/17)*alpha^2+(6/17)*alpha-36/17+((6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17)*x+x^2, x-alpha

Let beta be the first root of the quadratic; the other can be expressed in terms of alpha and beta

> 

p2alpha; alias(beta = RootOf(p2alpha, x))

(11/17)*alpha^5+(3/17)*alpha^4+(7/17)*alpha^3+(5/17)*alpha^2+(6/17)*alpha-36/17+((6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17)*x+x^2

alpha, beta

> 

p2alphabeta := evala(Factor(p2alpha, {alpha, beta}))

(x+(6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17+beta)*(x-beta)

The cubic can't be done in terms of alpha and beta, so let gamma be its first root.

> 

p3alpha; irreduc(p3alpha, {alpha, beta}); alias(gamma = RootOf(p3alpha, x))

x^3+((6/17)*alpha-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+15/17)*x^2+((6/17)*alpha^2-33/17+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3-(3/17)*alpha)*x+(5/17)*alpha^3+18/17+(7/17)*alpha^4+(6/17)*alpha^2+(3/17)*alpha^5-(3/17)*alpha

true

alpha, beta, gamma

It factors into a quadratic and the linear term, so we need another algebraic number as the first root of the quadratic

> 

p3alphabetagamma := evala(Factor(p3alpha, {alpha, beta, gamma}))

(x^2+(-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+(6/17)*alpha+15/17+gamma)*x+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3+(6/17)*alpha^2-(3/17)*alpha-33/17-(6/17)*gamma*alpha^5-(14/17)*gamma*alpha^4-(10/17)*gamma*alpha^3-(12/17)*gamma*alpha^2+(6/17)*gamma*alpha+(15/17)*gamma+gamma^2)*(x-gamma)

> 

p2alphabetagamma, p1alphabetagamma := op(p3alphabetagamma)

x^2+(-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+(6/17)*alpha+15/17+gamma)*x+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3+(6/17)*alpha^2-(3/17)*alpha-33/17-(6/17)*gamma*alpha^5-(14/17)*gamma*alpha^4-(10/17)*gamma*alpha^3-(12/17)*gamma*alpha^2+(6/17)*gamma*alpha+(15/17)*gamma+gamma^2, x-gamma

> 

alias(delta = RootOf(p2alphabetagamma, x))

alpha, beta, gamma, delta

> 

evala(Factor(p2alphabetagamma, {alpha, beta, delta, gamma}))

(x-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+(6/17)*alpha+15/17+gamma+delta)*(x-delta)

Can put it all together

> 

evala(AFactor(poly))

(x-delta)*(x+(6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17+beta)*(x-beta)*(x-alpha)*(x-gamma)*(x-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+(6/17)*alpha+15/17+gamma+delta)

This can be done all in one step with PolynomialTools:-Split - the splitting field K with four algebraic numbers is returned (but figuring out the nested RootOf's is painful).

> 

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

(x-delta)*(x+(6/17)*alpha^5+(14/17)*alpha^4+(10/17)*alpha^3+(12/17)*alpha^2+(11/17)*alpha-15/17+beta)*(x-beta)*(x-alpha)*(x-gamma)*(x-(6/17)*alpha^5-(14/17)*alpha^4-(10/17)*alpha^3-(12/17)*alpha^2+(6/17)*alpha+15/17+gamma+delta)

{alpha, beta, delta, gamma}

> 

evala(Minpoly(alpha, x)); dalpha := degree(%, x)

x^6-3*x+3

6

> 

evala(Minpoly(beta, x, {alpha})); dbeta := degree(%, x)

x^2+(6/17)*x*alpha^5+(14/17)*alpha^4*x+(10/17)*alpha^3*x+(12/17)*alpha^2*x+(11/17)*x*alpha-(15/17)*x+(11/17)*alpha^5+(3/17)*alpha^4+(7/17)*alpha^3+(5/17)*alpha^2+(6/17)*alpha-36/17

2

> 

evala(Minpoly(gamma, x, {alpha, beta})); dgamma := degree(%, x)

x^3-(6/17)*alpha^5*x^2-(14/17)*alpha^4*x^2-(10/17)*alpha^3*x^2-(12/17)*alpha^2*x^2+(6/17)*alpha*x^2+(15/17)*x^2+(3/17)*x*alpha^5+(7/17)*alpha^4*x+(5/17)*alpha^3*x+(6/17)*alpha^2*x-(3/17)*x*alpha-(33/17)*x+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3+(6/17)*alpha^2-(3/17)*alpha+18/17

3

> 

evala(Minpoly(delta, x, {alpha, beta, gamma})); ddelta := degree(%, x)

x^2-(6/17)*x*alpha^5-(14/17)*alpha^4*x-(10/17)*alpha^3*x-(12/17)*alpha^2*x+(6/17)*x*alpha+(15/17)*x+x*gamma+(3/17)*alpha^5+(7/17)*alpha^4+(5/17)*alpha^3+(6/17)*alpha^2-(3/17)*alpha-33/17-(6/17)*gamma*alpha^5-(14/17)*gamma*alpha^4-(10/17)*gamma*alpha^3-(12/17)*gamma*alpha^2+(6/17)*gamma*alpha+(15/17)*gamma+gamma^2

2

Order of Galois group.

> 

dalpha*dbeta*dgamma*ddelta

72

> 

NULL

Download Splits.mw

precision and algorithms...

dharr 8482
March 01 2023
3 2

That is a surprising result, and sufficiently far enough from the correct answer that I would call it a bug, even though it is just a numerical issue. Certainly since the matrix is "nice" (symmetric, diagonally dominant). As @vv says it depends on the precision, but certainly you don't need to go to 40 to get OK results. At Digits=14, M_MP.M_MP - C is recognizably close to zero and at Digits=16, the entries are all about 1e-6 or 1e-7.

In the general case, the algorithm for matrix functions is nontrivial, but Maple often has special cases for shape=symmetric or shape = hermitian, so then I think just applying the function to the eigenvalues as you did is straightforward.

@vv's comment about the determinant might be true (I agree), but it seems to imply that the reason that <0,1;0,0> doesn't have a square root is because its determinant is zero. That has more to do with the fact that it is non-diagonalizable; for example <1,0;0,0> (positive semidefinite) and <1,1;0,0> both have square roots and have determinant zero. 

Detecting edge cuts...

dharr 8482
February 25 2023
2 5

For the edge separator, Maple's IsCutSet can be used - it just checks for an increase in the number of components. Checking for a Cut is more complicated since just checking for 2 components can fail to detect extra removed edges that are in one of the components. Here is one way to do it.

[Edit - this actually tests for a cutset, which is a minimal Cut, so I've renamed it IsMinCut; see below for the test for a Cut]

> 

restart;

IsMinCut returns true if the edges are a minimal cut (cutset), i.e., removal of the edges cuts the graph into two components, and no edge can be added back without reconnecting the graph.

This version doesn't handle loops for efficiency, but could be modified to use the underlying graph and ignore loops (provided edges doesn't contain a loop).

> 

IsMinCut:=proc(G::GRAPHLN,edges::set)
       local partition,G1,G2;
       uses GraphTheory;
       if not IsConnected(G) or IsDirected(G) then error "Graph must be connected and undirected" end if;
       partition:=ConnectedComponents(DeleteEdge(G,edges,'inplace'=false));
       if nops(partition) <> 2 then return false end if;
       G1:=InducedSubgraph(G,partition[1]);
       G2:=InducedSubgraph(G,partition[2]);
       if (Edges(G) minus Edges(G1) minus Edges(G2)) = edges then
          true
       else
          false
       end if;
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(IsMinCut(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
IsMinCut(G,edges);
 

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

false

> 

edges:={{1,4},{1,5},{1,6}};
IsMinCut(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
IsMinCut(G,edges);

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

false

> 

NULL

Download IsMinCut.mw

a guess...

dharr 8482
February 24 2023
1 0

It looks as though when there are more unkowns than equations that Mathematica sets the extra ones first in order to be zero. So something like this achieves that result.

> 

restart;
> 

eqs := [2*c2 + c0 = 1, 6*c3 + 2*c1 = 2, 3*c2 + 12*c4 = 1];
oldvars:={c0,c1,c2,c3,c4};
zerovars:=oldvars[1..(nops(oldvars)-nops(eqs))];
vars:=oldvars minus zerovars;
sol1:=solve(eqs, vars);
zeroes:=zerovars=~0;
sol:=eval(sol1,zeroes) union zeroes;

[2*c2+c0 = 1, 6*c3+2*c1 = 2, 3*c2+12*c4 = 1]

{c0, c1, c2, c3, c4}

{c0, c1}

{c2, c3, c4}

{c2 = -(1/2)*c0+1/2, c3 = -(1/3)*c1+1/3, c4 = -1/24+(1/8)*c0}

{c0 = 0, c1 = 0}

{c0 = 0, c1 = 0, c2 = 1/2, c3 = 1/3, c4 = -1/24}

> 

NULL

Download solve.mw

ListTools:-Occurrences...

dharr 8482
February 24 2023
0 5

There are many ways; here is one. (Note that CycleBasis doesn't give all cycles in the graph)

> 

restart
> 

with(GraphTheory); G := SpecialGraphs:-OctahedronGraph()

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

> 

Cycles := CycleBasis(G)

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

> 

ListTools:-Occurrences(3, map(numelems, Cycles))

4

> 

ListTools:-Occurrences(4, map(numelems, Cycles))

3

> 

NULL

Download countcycles.mw

variation...

dharr 8482
February 23 2023
2 2

A variation on @acer 's and @Kitonum's responses that bypasses radicals or floating points.

> 

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)}});

{k[1]+k[2]+k[3]+k[4]+k[5]+k[6]+k[7]+k[8]+k[9]+k[10], k[1]*alpha^3+k[2]*alpha^2*beta+k[5]*alpha^2*gamma+k[3]*alpha*beta^2+k[6]*alpha*beta*gamma+k[8]*alpha*gamma^2+k[4]*beta^3+k[7]*beta^2*gamma+k[9]*beta*gamma^2+k[10]*gamma^3, k[1]*beta^3+k[2]*beta^2*gamma+k[5]*beta^2*alpha+k[3]*beta*gamma^2+k[6]*alpha*beta*gamma+k[8]*beta*alpha^2+k[4]*gamma^3+k[7]*gamma^2*alpha+k[9]*gamma*alpha^2+k[10]*alpha^3, k[1]*gamma^3+k[2]*gamma^2*alpha+k[5]*gamma^2*beta+k[3]*gamma*alpha^2+k[6]*alpha*beta*gamma+k[8]*gamma*beta^2+k[4]*alpha^3+k[7]*alpha^2*beta+k[9]*alpha*beta^2+k[10]*beta^3}

> 

eqs2:=solve(evala(eqs));

{k[1] = (13/24)*k[4]-(7/24)*k[5]-(1/8)*k[6]-(1/6)*k[7]-(5/24)*k[8]+(5/12)*k[10], k[2] = (37/8)*k[4]-(7/8)*k[5]-(11/8)*k[6]-(3/2)*k[7]-(13/8)*k[8]+(17/4)*k[10], k[3] = -(143/24)*k[4]+(5/24)*k[5]+(3/8)*k[6]+(5/6)*k[7]+(7/24)*k[8]-(7/12)*k[10], k[4] = k[4], k[5] = k[5], k[6] = k[6], k[7] = k[7], k[8] = k[8], k[9] = -(5/24)*k[4]-(1/24)*k[5]+(1/8)*k[6]-(1/6)*k[7]+(13/24)*k[8]-(61/12)*k[10], k[10] = k[10]}

> 

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}

> 

NULL

Download algeqns.mw

Condition for stable roots...

dharr 8482
February 20 2023
1 2

You can use PolynomialTools:-Hurwitz to find the conditions for stable roots, i.e for roots with negative real parts.

> 

restart
> 

a1 := 1; g[2] := -S0*eta[2]+b+r; a[1] := c+epsilon+g[2]; a[2] := epsilon*g[2]+c*(epsilon+g[2])-S0*b*eta[3]; a[3] := -S0*Pi*b*eta[1]+c*(-S0*b*eta[3]+epsilon*g[2])

1

> 

eq := lambda^3+lambda^2*a[1]+lambda*a[2]+a[3]

lambda^3+lambda^2*(-S0*eta[2]+b+c+epsilon+r)+lambda*(epsilon*(-S0*eta[2]+b+r)+c*(-S0*eta[2]+b+epsilon+r)-S0*b*eta[3])-S0*Pi*b*eta[1]+c*(-S0*b*eta[3]+epsilon*(-S0*eta[2]+b+r))

> 

PolynomialTools:-Hurwitz(eq, lambda, 's', 'g')

FAIL

> 

s

[-lambda/(S0*eta[2]-b-c-epsilon-r), (S0*eta[2]-b-c-epsilon-r)^2*lambda/(S0^2*b*eta[2]*eta[3]+S0^2*c*eta[2]^2+S0^2*epsilon*eta[2]^2+Pi*S0*b*eta[1]-S0*b^2*eta[3]-2*S0*b*c*eta[2]-2*S0*b*epsilon*eta[2]-S0*b*epsilon*eta[3]-S0*b*r*eta[3]-S0*c^2*eta[2]-2*S0*c*epsilon*eta[2]-2*S0*c*r*eta[2]-S0*epsilon^2*eta[2]-2*S0*epsilon*r*eta[2]+b^2*c+b^2*epsilon+b*c^2+2*b*c*epsilon+2*b*c*r+b*epsilon^2+2*b*epsilon*r+c^2*epsilon+c^2*r+c*epsilon^2+2*c*epsilon*r+c*r^2+epsilon^2*r+epsilon*r^2), (S0^2*b*eta[2]*eta[3]+S0^2*c*eta[2]^2+S0^2*epsilon*eta[2]^2+Pi*S0*b*eta[1]-S0*b^2*eta[3]-2*S0*b*c*eta[2]-2*S0*b*epsilon*eta[2]-S0*b*epsilon*eta[3]-S0*b*r*eta[3]-S0*c^2*eta[2]-2*S0*c*epsilon*eta[2]-2*S0*c*r*eta[2]-S0*epsilon^2*eta[2]-2*S0*epsilon*r*eta[2]+b^2*c+b^2*epsilon+b*c^2+2*b*c*epsilon+2*b*c*r+b*epsilon^2+2*b*epsilon*r+c^2*epsilon+c^2*r+c*epsilon^2+2*c*epsilon*r+c*r^2+epsilon^2*r+epsilon*r^2)*lambda/((S0*eta[2]-b-c-epsilon-r)*(Pi*S0*b*eta[1]+S0*b*c*eta[3]+S0*c*epsilon*eta[2]-b*c*epsilon-c*epsilon*r))]

These must be positive:

> 

conds := map(coeff, s, lambda)

[-1/(S0*eta[2]-b-c-epsilon-r), (S0*eta[2]-b-c-epsilon-r)^2/(S0^2*b*eta[2]*eta[3]+S0^2*c*eta[2]^2+S0^2*epsilon*eta[2]^2+Pi*S0*b*eta[1]-S0*b^2*eta[3]-2*S0*b*c*eta[2]-2*S0*b*epsilon*eta[2]-S0*b*epsilon*eta[3]-S0*b*r*eta[3]-S0*c^2*eta[2]-2*S0*c*epsilon*eta[2]-2*S0*c*r*eta[2]-S0*epsilon^2*eta[2]-2*S0*epsilon*r*eta[2]+b^2*c+b^2*epsilon+b*c^2+2*b*c*epsilon+2*b*c*r+b*epsilon^2+2*b*epsilon*r+c^2*epsilon+c^2*r+c*epsilon^2+2*c*epsilon*r+c*r^2+epsilon^2*r+epsilon*r^2), (S0^2*b*eta[2]*eta[3]+S0^2*c*eta[2]^2+S0^2*epsilon*eta[2]^2+Pi*S0*b*eta[1]-S0*b^2*eta[3]-2*S0*b*c*eta[2]-2*S0*b*epsilon*eta[2]-S0*b*epsilon*eta[3]-S0*b*r*eta[3]-S0*c^2*eta[2]-2*S0*c*epsilon*eta[2]-2*S0*c*r*eta[2]-S0*epsilon^2*eta[2]-2*S0*epsilon*r*eta[2]+b^2*c+b^2*epsilon+b*c^2+2*b*c*epsilon+2*b*c*r+b*epsilon^2+2*b*epsilon*r+c^2*epsilon+c^2*r+c*epsilon^2+2*c*epsilon*r+c*r^2+epsilon^2*r+epsilon*r^2)/((S0*eta[2]-b-c-epsilon-r)*(Pi*S0*b*eta[1]+S0*b*c*eta[3]+S0*c*epsilon*eta[2]-b*c*epsilon-c*epsilon*r))]

And these must have real parts equal to zero, which they do.

> 

map(coeff, s, lambda, 0)

[0, 0, 0]

We can analyse the conditions a little further. All of the following need to be positive

> 

c1 := 1/conds[1]; c2 := factor(normal(c1^2/conds[2])); c3 := factor(simplify(c2/(c1*conds[3])))

-S0*eta[2]+b+c+epsilon+r

S0^2*b*eta[2]*eta[3]+S0^2*c*eta[2]^2+S0^2*epsilon*eta[2]^2+Pi*S0*b*eta[1]-S0*b^2*eta[3]-2*S0*b*c*eta[2]-2*S0*b*epsilon*eta[2]-S0*b*epsilon*eta[3]-S0*b*r*eta[3]-S0*c^2*eta[2]-2*S0*c*epsilon*eta[2]-2*S0*c*r*eta[2]-S0*epsilon^2*eta[2]-2*S0*epsilon*r*eta[2]+b^2*c+b^2*epsilon+b*c^2+2*b*c*epsilon+2*b*c*r+b*epsilon^2+2*b*epsilon*r+c^2*epsilon+c^2*r+c*epsilon^2+2*c*epsilon*r+c*r^2+epsilon^2*r+epsilon*r^2

-Pi*S0*b*eta[1]-S0*b*c*eta[3]-S0*c*epsilon*eta[2]+b*c*epsilon+c*epsilon*r

NULL

NULL

Download Hurwitz.mw

not easy to do...

dharr 8482
February 20 2023
3 2

I also would like an easy way to do this. I think evala(Reduce()) and evala(Normal()) want to have the the answer in the same algebraic field extension so don't do this. It can be done somewhat awkwardly through Minpoly; perhaps you had another objection to this. 

> 

restart
> 

q := -4*RootOf(_Z^3-3*_Z^2-10*_Z-1)^2*(1/5)+19*RootOf(_Z^3-3*_Z^2-10*_Z-1)*(1/5)+3/5; simplify(fnormal(evalf([allvalues(q)])))

-(4/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1)^2+(19/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1)+3/5

[-.519484859, .198874518, -9.679389674]

> 

RootOf(evala(Minpoly(q, x))); simplify(fnormal(evalf([allvalues(%)])))

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

[.198874521, -.519484847, -9.679389673]

> 

q2:=-45658*RootOf(37*_Z^6 - 382*_Z^5 + 1388*_Z^4 - 2188*_Z^3 + 1475*_Z^2 - 406*_Z + 37, index = 6)^5 + 417257*RootOf(37*_Z^6 - 382*_Z^5 + 1388*_Z^4 - 2188*_Z^3 + 1475*_Z^2 - 406*_Z + 37, index = 6)^4 - 1252087*RootOf(37*_Z^6 - 382*_Z^5 + 1388*_Z^4 - 2188*_Z^3 + 1475*_Z^2 - 406*_Z + 37, index = 6)^3 + 1463384*RootOf(37*_Z^6 - 382*_Z^5 + 1388*_Z^4 - 2188*_Z^3 + 1475*_Z^2 - 406*_Z + 37, index = 6)^2 - 558475*RootOf(37*_Z^6 - 382*_Z^5 + 1388*_Z^4 - 2188*_Z^3 + 1475*_Z^2 - 406*_Z + 37, index = 6) + 69230;

-45658*RootOf(37*_Z^6-382*_Z^5+1388*_Z^4-2188*_Z^3+1475*_Z^2-406*_Z+37, index = 6)^5+417257*RootOf(37*_Z^6-382*_Z^5+1388*_Z^4-2188*_Z^3+1475*_Z^2-406*_Z+37, index = 6)^4-1252087*RootOf(37*_Z^6-382*_Z^5+1388*_Z^4-2188*_Z^3+1475*_Z^2-406*_Z+37, index = 6)^3+1463384*RootOf(37*_Z^6-382*_Z^5+1388*_Z^4-2188*_Z^3+1475*_Z^2-406*_Z+37, index = 6)^2-558475*RootOf(37*_Z^6-382*_Z^5+1388*_Z^4-2188*_Z^3+1475*_Z^2-406*_Z+37, index = 6)+69230

> 

num := evalf(q2)

26532.730

> 

convert(RootOf(evala(Minpoly(q2, x)), num), RootOf, form = index)

RootOf(37*_Z^6-7304346*_Z^5+477422219475*_Z^4-12741284944716948*_Z^3+145415493111187762668*_Z^2-720012242195396824623282*_Z+1254681647187128843079859317, index = 4)

> 

NULL

Download RootOf.mw

shape=symmetric...

dharr 8482
February 19 2023
0 0


Download Matrix.mw

See here for simplifying the entry (and also the coversion to the matrix).

bad RootOf...

dharr 8482
February 18 2023
2 12

The index=real[4] in the RootOf looks wrong. [Edit: I submitted an SCR].

> 

restart
> 

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:

You can use the explicit option so you don't get the RootOf, but then you get only four solutions.

> 

sol1 := solve({expr, a >= 0, b >= 0}, {a, b}, allsolutions, explicit)

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

There is something wrong with the RootOf solution, which has index=real[4] instead of index = positive integer
Under the assumption that it was supposed to be an integer, we try various integers and with index=1 one of the RootOfs is resolved, but the other is now not correct (should be 1, but evalf does not give this)

> 

sol2 := solve({expr, a >= 0, b >= 0}, {a, b}, allsolutions)
> 

sol3 := convert(subs(real[4] = 1, [sol2]), radical)

[{a = 1, b = 0}, {a = RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 3), b = 3^(1/2)-1}, {a = 0, b = 1}, {a = 3^(1/2)-1, b = 1}, {a = 1, b = 1}, {a = 1/2+(1/2)*3^(1/2), b = 1/2+(1/2)*3^(1/2)}]

Let's go back and just remove the index from the RootOf

> 

ro := RootOf(op(1, eval(a, sol2[2])))

and the correct answer is the first of these values.

> 

allvalues(ro)

1, RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 1), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 2), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 3), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 4), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 5), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 6), RootOf(_Z^7+(2+3^(1/2))*_Z^6+(11*3^(1/2)-18)*_Z^5+(-5*3^(1/2)+4)*_Z^4+(-144*3^(1/2)+251)*_Z^3+(39*3^(1/2)-66)*_Z^2+(981*3^(1/2)-1698)*_Z+873*3^(1/2)-1512, index = 7)

> 

NULL

Download solve.mw

_EnvUseHeavisideAsUnitStep:=true;...

dharr 8482
February 16 2023
3 1
> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2022.2, Windows 10, October 23 2022 Build ID 1657361`

> 

_EnvUseHeavisideAsUnitStep:=true;
expr:=(s+exp(-Pi*s)-exp(-2*Pi*s))/(s*(s^2+2*s+2));
inttrans:-invlaplace(expr,s,t);
Y:=convert(%,piecewise);
#remove all entries in piecwise which has undefined

_EnvUseHeavisideAsUnitStep := true

expr := (s+exp(-Pi*s)-exp(-2*Pi*s))/(s*(s^2+2*s+2))

exp(-t)*sin(t)+(1/2)*(-1+exp(-t+2*Pi)*(cos(t)+sin(t)))*Heaviside(t-2*Pi)+(1/2)*(1+exp(-t+Pi)*(cos(t)+sin(t)))*Heaviside(t-Pi)

piecewise(t < Pi, exp(-t)*sin(t), t < 2*Pi, exp(-t)*sin(t)+1/2+(1/2)*exp(-t+Pi)*(cos(t)+sin(t)), 2*Pi <= t, exp(-t)*sin(t)+(1/2)*exp(-t+2*Pi)*(cos(t)+sin(t))+(1/2)*exp(-t+Pi)*(cos(t)+sin(t)))

> 

NULL

Download remove_undefined.mw

numeric values required...

dharr 8482
February 11 2023
0 0

If I assign some values such as h_n:=[5,4], etc (for np=2) then it works.

points.mw

or a variation more suitable to a large number of points.

points2.mw

some things...

dharr 8482
February 11 2023
1 0

Your variable s__v_0 actually had an unseen multiplication in. I replaced it with just sv0. Some subscripts were made differently in different places, for example v[1] or v__1; I changed them all to the double-underline form. In your differential equation mu was strange (from a font perhaps instead of typing mu and using escape), and different from the mu in the parameters list. You had two alphas in the parameter list.

You now need to resolve the fact that yor equations have Q(t) but there is no differential equation for Q(t). You have some parameters in your parameter list that are not in your equations - I would have removed them but you give them values later. So next step is to resolve these inconsistencies.

Covid19_Simulation2.mw

First 34 35 36 37 38 39 40 Last Page 36 of 83