dharr

Dr. David Harrington

5968 Reputation

21 Badges

19 years, 295 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
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 answers submitted by dharr

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

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;

-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

 

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

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

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

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

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

 

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

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

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

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

I changed the rho's in the solution to rho__v's (as in the Eqs) and it worked.
solution_check.mw

I agree with Carl. Perhaps the following helps to illustrate the point.

restart

NULLN := (-t^2-1)*u^2-2*(t^2+1)^2*u+2*tNULL

(-t^2-1)*u^2-2*(t^2+1)^2*u+2*t

A := u^2*a[2]+u*a[1]+a[0]; B := u^2*b[2]+u*b[1]+b[0]

u^2*a[2]+u*a[1]+a[0]

u^2*b[2]+u*b[1]+b[0]

Q := collect(B*(diff(A, u))-A*(diff(B, u)), u)

(-a[1]*b[2]+a[2]*b[1])*u^2+(-2*a[0]*b[2]+2*a[2]*b[0])*u-a[0]*b[1]+b[0]*a[1]

We want to find when Q = N for any u. Since doubling the a's can be compensated for by halving the b's we can fix the scale by arbitrarily choosing a value for b[0], say 1. Or we can just write the solution for the other variables in terms of b[0]. Since we have 3 equations in 6 unknowns, we can extend this logic to another two variables aside from b[0], i.e., as @Carl Love noted we need to choose three variables to solve for in terms of the other three. For example, the following gives a[1], a[2] and b[2] in terms of the other 3. (Solve/identity is just a shortcut for the 3 equations equating the coefficients.

solve(identity(N = Q, u), {a[1], a[2], b[2]})

{a[1] = (a[0]*b[1]+2*t)/b[0], a[2] = -(1/2)*(t^4*a[0]*b[1]+2*t^5-t^2*a[0]*b[0]+2*t^2*a[0]*b[1]+4*t^3-a[0]*b[0]+a[0]*b[1]+2*t)/(t*b[0]), b[2] = -(1/2)*(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1])/t}

Suppose we give a 4th variable to solve for - now we see the solution contains b[0] = b[0], meaning we can choose b[0] arbitrarily, the others are as above.

solve(identity(N = Q, u), {a[1], a[2], b[0], b[2]})

{a[1] = (a[0]*b[1]+2*t)/b[0], a[2] = -(1/2)*(t^4*a[0]*b[1]+2*t^5-t^2*a[0]*b[0]+2*t^2*a[0]*b[1]+4*t^3-a[0]*b[0]+a[0]*b[1]+2*t)/(t*b[0]), b[0] = b[0], b[2] = -(1/2)*(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1])/t}

In this case, Maple chose to make a[1] arbitrary; we get the same as if we had chosen a[2], b[1] and b[2] to solve for.

solve(identity(N = Q, u), {a[1], a[2], b[1], b[2]})

{a[1] = a[1], a[2] = -(1/2)*(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1])/t, b[1] = -(-a[1]*b[0]+2*t)/a[0], b[2] = (1/2)*(-t^4*a[1]*b[0]+2*t^5+t^2*a[0]*b[0]-2*t^2*a[1]*b[0]+4*t^3+a[0]*b[0]-a[1]*b[0]+2*t)/(t*a[0])}

NULL

Download Huang.mw

I made the opposite assumption to @acer (that you started with the equations), and kept implicitplot butdid them all in one, which is simple here despite the inefficiency.

restart

Load the necessary packages

with(plots); with(LinearAlgebra)

Define the equations to solve and plot

eq1 := 2*x-y+z = -5; eq2 := y+3*z = 7; eq3 := z = 2

Convert to the corresponding Matrix of coefficients and Vector (because we shouldn't have to enter equivalent information twice, and less likely to make a mistake)

A, b := GenerateMatrix([eq1, eq2, eq3], [x, y, z])

Matrix(%id = 36893490370101429060), Vector[column](%id = 36893490370101428940)

sol := LinearSolve(A, b)

Vector[column](%id = 36893490370101427364)

We can do all equations in one implicitplot. (I like the brighter colors, so I left the quotes off.)

implicitplot3d([eq1, eq2, eq3], x = -10 .. 10, y = -10 .. 10, z = -10 .. 10, color = [red, blue, green], style = surface, title = "3 D Plot of Linear System")

``

Download linear_systems.mw

1 2 3 4 5 6 7 Last Page 1 of 63