## dharr

Dr. David Harrington

## 6416 Reputation

20 years, 23 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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.

## solution...

I agree with @nm; isolve is frequently weak, especially here where the equations are nonlinear. Here is another ad-hoc way that works.

 > restart;
 > eq1 := a*b + c = 2020; eq2 := b*c + a = 2021;

Eliminate b, then isolve can handle the remaining equation in a and c

 > eliminate({eq1, eq2}, {b}); beq,eq3:=%[];

 > [isolve(eq3)]; eq3solns:=remove(has,%,a=0); #avoid division by zero

Substitute back into the equation for b and select the solutions where b is an integer

 > map(x->eval(beq,x) union x,eq3solns); select(x->eval(b,x)::integer,%);

 >

## local variable...

If you change local A:=0,B to global A:=0,B it works as expected. When you print out A, you are printing out the local A, and then you use that value to determine the value to be saved.

The help for read says "If the file is in Maple language format, the statements in the file are read and executed as if they were being entered into Maple interactively", which I take to mean that A := 6 assigns a value to global A.

You could adjust by being more explicit about the locals and globals:

 > restart;
 > currentdir(); try  #remove TMP.m first time     FileTools:-Remove("TMP.m"); catch:     NULL; end try; #this function is called many times. foo:=proc()  local A:=0,B;  A:=A+10;  if FileTools:-Exists("TMP.m") then      print("TMP.m exists, will read its content");      B:=A; #save copy      read "TMP.m";      print("Read old value of A from file",:-A);      A:= :-A + B; # update running total      print("Now A is ",A," will now save this new value");           save A,"TMP.m";   else      print("TMP.m do not exist, first time saving to it");      save A,"TMP.m";   fi;   print("A=",A, :-A); end proc:

 > foo();

 > foo();

 > foo();

 > foo();

 >

## if statement...

Note that sigma(a)*sigma(b) is not the same as sigma(a*b), which it seemed to me you were assuming. I think this answers your question about the if statement.

 >
 >
 >

 >

 >

 > if p>0 then   print(abundant) elif p=0 then   print(perfect) else   print(deficient) end if;

 >

## variation...

Here's a variation on @Ronan's answer.

 > restart;
 > for n to 4 do    A[n]:=Matrix(n,n,(i,j)->local m;                            ifelse(i+j

## bug reporting...

I agree with your analysis, and that it is a bug. To report it, you can (logged into Mapleprimes) choose the "more" tab, and then "submit software change request". I usually copy the mapleprimes link and paste that in the box where you are asked for any further information.

I don't know the specific answer to question 1. If you extend to Gamma negative you don't see a region. Anyway, as you say, it gives the correct number of solutions in the different regions shown. Other explanations below.

 > restart;

Here @acer's analysis

 > Eq := -8*(rho + 1)^4*Lambda^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*Lambda^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*Lambda^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*Lambda + Gamma^2*(rho + 1)*(rho - 1)^2;

 > with(RootFinding:-Parametric):
 > m := CellDecomposition([Eq=0, Lambda>0, Gamma>0, rho>-1, rho<1], [Lambda]):
 > CellPlot(m, samplepoints, view=[0..10,-2..2], size=[400,300]);

Number of solutions in regions 1..5

 > NumberOfSolutions(m)[1..5];

The chosen sample points within the region (as shown on the plot)

 > m:-SamplePoints;

Value of Lambda at sample point 2

 > SampleSolutions(m,2);

 > SampleSolutions(m,3);

The plot3d that I think you are worried about. For the region Lambda>0, there is everywhere a solution (the top sheet), in agreement with the analysis above that shows 1 real solution in regions 2and 3

 > plot3d([allvalues(RootOf(Eq,Lambda))],rho=-1..1,Gamma=0..5, grid = [20,20]);

 >

## No...

If you mean multiply by some function of the variables then the answer is no. The answer is somehow "closer" if you change the red k/2 to -k/2. What are the assumption on k? If one is allowed to multiply by matrices, then the answer might be different. (Matrices not showing correctly on Mapleprimes.

 >
 >
 >

 >

To get the off diagonals of L to be those of X, we need to multiply by

 >

 >

Diagonals have different sums, so cannot be the same.

 >

 >

## PolynomialSystem...

SolveTools:-PolynomialSystem with engine=triade can solve this. PolynomialSystems allows choice of different engines, one of which is groebner, so you can experiment if it gets stuck.

solve_test.mw

## ListTools:-Categorize...

To find the non-somorphic ones I use ListTools:-Categorize. I think IsIsomorphic already uses canonical labelling internally anyway, so is not required explicitly.

 > restart;
 > with(GraphTheory):
 > graph_list := [GraphTheory:-CompleteGraph(3), GraphTheory:-PathGraph(3),Graph({{a,b},{b,c},{c,a}})]:
 > classes:=[ListTools:-Categorize(IsIsomorphic,graph_list)];

Take first one from each class

 > nonisomorphic:=map2(op,1,classes);

If you would like them generically labelled (vertices 1,2,etc) then (or this could be done on the original list)

 > map(G->Graph(op(4,G)),nonisomorphic);

 >

## missing derivatives...

It's hard to know what form linjeelementer expects for an ODE, but if I assume it's the standard Maple form, then the error message is clear - you have given a differential equation without a derivative. Probably what you have given is the growth rate diff(y(t),t) so then you want

`ODE := diff(y(t), t) = -0.000032229*y(t)^2 + 0.065843*y(t) - 15.103`

## complete algorithm...

Here's a complete algorithm, assuming I've understood your objective correctly. Your C and C' can just be combined into a single combined C that I deal with - the counts for C={1}, C'={2,3,4} are the same as for C={1,2}, C'={3,4}, so I just count {1,2,3,4}.

Each row and column-choice combination is a pattern for which a count is recorded, so there are no unsuccessful searches.

Older version:

Edit: new much faster version takes advantage of the fact that there are many duplicate rows. Now L=10 can be done in about a second on my machine:

 > restart;
 > N := 10^5: L := 5: M := LinearAlgebra:-RandomMatrix(N, L, generator=0..1):

All subsets that can be C union C' with C, C' disjoint and nonempty

 > Lset:={\$1..L}: cs:=[map(convert,combinat:-powerset(L) minus {{},Lset,seq({i},i=1..L)},list)[]];

Procedure that counts the possibilities for the matrix M and the column choices in cs.

 > count := proc(M, cs);         local counts, C, n, N, P, row, MT, nrows;         N := upperbound(M)[1];         # prescan Matrix for duplicate rows and count them         MT := table('sparse'):         for n to N do                 row := [seq(M[n])];                 MT[row]++;         end do:         counts := table('sparse');         for row, nrows in eval(MT) do                 for C in cs do                         P := row[C];                         counts[P, C] += nrows;                 end do;         end do;         eval(counts); end proc:
 > counts := CodeTools:-Usage(count(M, cs)):

memory used=56.16MiB, alloc change=5.00MiB, cpu time=47.00ms, real time=145.00ms, gc time=15.62ms

How many times does the pattern [1,1,0] occur in columns [1, 2, 4]?

 > counts[[1,1,0], [1,2,4]];

Total counts in all cases

Try L=10;

 > N := 10^5: L := 10: M := LinearAlgebra:-RandomMatrix(N, L, generator=0..1):
 > Lset := {\$1..L}: cs := [map(convert,combinat:-powerset(L) minus {{},Lset,seq({i},i=1..L)},list)[]]: nops(cs);

 > counts := CodeTools:-Usage(count(M, cs)):

memory used=168.70MiB, alloc change=24.18MiB, cpu time=765.00ms, real time=1.10s, gc time=78.12ms

 > counts[[1,0,1,0,1], [1,4,5,6,9]];

 >

2015 compatible version:
mmcdara8a.mw

## analytical solution...

An accurate analytical solution for approximate (floating point) coefficients varying over many orders of magnitude is unlikely, so it would be better if you entered the system using rationals (from which the floating values presumably came). Maple begins by converting the floating point values to rationals, but an analytical solution can then be obtained using method = laplace.

## gfun...

Here is your sin(x^2) example, using guessgf directly when no odes are involved.

 > restart;
 > convert(series(sin(x^2),x,40),polynom); cofs:=[seq(coeff(%,x,i),i=0..39)];

 > gfun:-guessgf(cofs,x);

 >

 >
 >

 >

 >