Advanced Search

dharr

Dr. David Harrington

6416 Reputation

21 Badges

20 years, 23 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

solution...

dharr 6416
June 23 2024
3 1

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;
 

a*b+c = 2020

b*c+a = 2021

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

> 

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

[{b = -(c-2020)/a}, {a^2-c^2-2021*a+2020*c}]

{b = -(c-2020)/a}, {a^2-c^2-2021*a+2020*c}

> 

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

[{a = 0, c = 0}, {a = 0, c = 2020}, {a = 673, c = 674}, {a = 673, c = 1346}, {a = 896, c = 900}, {a = 896, c = 1120}, {a = 1125, c = 900}, {a = 1125, c = 1120}, {a = 1348, c = 674}, {a = 1348, c = 1346}, {a = 2021, c = 0}, {a = 2021, c = 2020}]

[{a = 673, c = 674}, {a = 673, c = 1346}, {a = 896, c = 900}, {a = 896, c = 1120}, {a = 1125, c = 900}, {a = 1125, c = 1120}, {a = 1348, c = 674}, {a = 1348, c = 1346}, {a = 2021, c = 0}, {a = 2021, c = 2020}]

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,%);

[{a = 673, b = 2, c = 674}, {a = 673, b = 674/673, c = 1346}, {a = 896, b = 5/4, c = 900}, {a = 896, b = 225/224, c = 1120}, {a = 1125, b = 224/225, c = 900}, {a = 1125, b = 4/5, c = 1120}, {a = 1348, b = 673/674, c = 674}, {a = 1348, b = 1/2, c = 1346}, {a = 2021, b = 2020/2021, c = 0}, {a = 2021, b = 0, c = 2020}]

[{a = 673, b = 2, c = 674}, {a = 2021, b = 0, c = 2020}]

> 

NULL

Download isolve.mw

local variable...

dharr 6416
June 22 2024
3 1

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:
 

"C:\Users\dharr\OneDrive - University of Victoria\Desktop"

> 

foo();

"TMP.m do not exist, first time saving to it"

"A=", 10, A

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 10

"Now A is ", 20, " will now save this new value"

"A=", 20, 10

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 20

"Now A is ", 30, " will now save this new value"

"A=", 30, 20

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 30

"Now A is ", 40, " will now save this new value"

"A=", 40, 30

> 

 


 

Download T.mw

if statement...

dharr 6416
June 19 2024
0 3

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.

> 

restart
> 

with(numtheory)
> 

sigma(12)*sigma(6)

336

> 

sigma(12*6)

195

> 

p := sigma(12)-2*12

4

> 

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

abundant

> 

``

Download abundant_edge.mw

variation...

dharr 6416
June 15 2024
2 4

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<n+2,
                                  add(f[i-m+j-1]*binomial(i-1,m)*(-1)^m,
                                      m=0..i-1),
                                  0));
end do;

Vector(1, {(1) = f[1]})

Matrix(2, 2, {(1, 1) = f[1], (1, 2) = f[2], (2, 1) = f[2]-f[1], (2, 2) = 0})

Matrix(3, 3, {(1, 1) = f[1], (1, 2) = f[2], (1, 3) = f[3], (2, 1) = f[2]-f[1], (2, 2) = f[3]-f[2], (2, 3) = 0, (3, 1) = f[3]-2*f[2]+f[1], (3, 2) = 0, (3, 3) = 0})

Matrix(%id = 36893628080372195684)

Download binomial.mw

bug reporting...

dharr 6416
June 15 2024
1 0

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.

comments...

dharr 6416
June 04 2024
2 3

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;

-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];

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

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

> 

m:-SamplePoints;

[[Gamma = 1, rho = -2], [Gamma = 1, rho = 0], [Gamma = 2, rho = 0], [Gamma = 1, rho = 2], [Gamma = 2, rho = 2]]

Value of Lambda at sample point 2

> 

SampleSolutions(m,2);

[[Lambda = .5622851345]]

> 

SampleSolutions(m,3);

[[Lambda = .5161075747]]

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

> 

NULL

Download cells.mw

 

No...

dharr 6416
May 27 2024
1 5

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.

> 

restart
> 

with(LinearAlgebra)
> 

L := Matrix([[(1+I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k), I*rho*(exp((1/2)*k)-exp(-(1/2)*k))/lambda], [I*rho*(exp((-k)*(1/2))-exp((1/2)*k))/lambda, (1-I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp((1/2)*k)]])

Matrix(%id = 36893490028717500884)

> 

X := Matrix([[-I*lambda*(1/2)-(1/2)*w, -rho], [rho, I*lambda*(1/2)+(1/2)*w]])

Matrix(%id = 36893490028738504332)

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

> 

F := L*lambda/(I*(exp(-(1/2)*k)-exp((1/2)*k)))

Matrix(%id = 36893490028738496500)

> 

F__2 := `assuming`([simplify(map(convert, F, exp))], [k > 0, k < (1/2)*Pi])

Matrix(%id = 36893490028738516380)

Diagonals have different sums, so cannot be the same.

> 

simplify(Trace(F__2)); Trace(X)

(w*(exp(k)-1)*tan((1/2)*k)+I*(exp(k)+1)*lambda)/(exp(k)-1)

0

> 

NULL

Download XintoL2.mw

PolynomialSystem...

dharr 6416
May 26 2024
0 4

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...

dharr 6416
May 21 2024
3 6

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

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

Take first one from each class

> 

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

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

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

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

> 

NULL

Download graphlist.mw

missing derivatives...

dharr 6416
May 19 2024
1 0

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...

dharr 6416
May 18 2024
2 2

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:

Download mmcdara4.mw

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)[]];

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

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]];

12344

Total counts in all cases

> 

add(eval(counts));
N*nops(cs);

2500000

2500000

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

1012

> 

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]];

3078

> 

add(eval(counts));

101200000

> 

NULL

Download mmcdara8.mw

2015 compatible version:
mmcdara8a.mw

 

analytical solution...

dharr 6416
May 15 2024
2 0

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.

Download Set_of_Linear_DEs.mw

gfun...

dharr 6416
May 14 2024
2 1

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

x^2-(1/6)*x^6+(1/120)*x^10-(1/5040)*x^14+(1/362880)*x^18-(1/39916800)*x^22+(1/6227020800)*x^26-(1/1307674368000)*x^30+(1/355687428096000)*x^34-(1/121645100408832000)*x^38

[0, 0, 1, 0, 0, 0, -1/6, 0, 0, 0, 1/120, 0, 0, 0, -1/5040, 0, 0, 0, 1/362880, 0, 0, 0, -1/39916800, 0, 0, 0, 1/6227020800, 0, 0, 0, -1/1307674368000, 0, 0, 0, 1/355687428096000, 0, 0, 0, -1/121645100408832000, 0]

> 

gfun:-guessgf(cofs,x);

[sin(x^2), ogf]

> 

``

Download gfun.mw

combine...

dharr 6416
May 13 2024
3 1
> 

restart
> 

s := add(cat(a, i)*cos(x)^i, i = 0 .. 4)

a0+a1*cos(x)+a2*cos(x)^2+a3*cos(x)^3+a4*cos(x)^4

> 

combine(s)collect(%, cos)

a0+(1/2)*a2+(3/8)*a4+a1*cos(x)+(1/2)*a2*cos(2*x)+(3/4)*a3*cos(x)+(1/4)*a3*cos(3*x)+(1/8)*a4*cos(4*x)+(1/2)*a4*cos(2*x)

(a1+(3/4)*a3)*cos(x)+((1/2)*a2+(1/2)*a4)*cos(2*x)+a0+(1/2)*a2+(3/8)*a4+(1/4)*a3*cos(3*x)+(1/8)*a4*cos(4*x)

> 

NULL

Download combine.mw

sequentially...

dharr 6416
May 09 2024
1 0

Yes, you need to do them sequentially - the help page warns it won't work doing them together in one alias command.

alias_in_alias.mw

3 4 5 6 7 8 9 Last Page 5 of 67