Advanced Search

dharr

Dr. David Harrington

8200 Reputation

22 Badges

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

No...

dharr 8200
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 8200
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 8200
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 8200
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 8200
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 8200
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 8200
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 8200
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 8200
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

rho and rho__v...

dharr 8200
May 08 2024
1 1

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

Solve for three variables....

dharr 8200
April 28 2024
2 7

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

from eqns to Matrix...

dharr 8200
April 27 2024
2 0

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

Fourier-Bessel...

dharr 8200
April 27 2024
1 0

Depends a bit on how you want to enter it.

> 

restart
> 

A := m -> int(xi*(1 - xi^2)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric)/
          int(xi*BesselJ(0, BesselJZeros(0, m)*xi)^2, xi = 0 .. 1, numeric);

proc (m) options operator, arrow; (int(xi*(-xi^2+1)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric))/(int(xi*BesselJ(0, BesselJZeros(0, m)*xi)^2, xi = 0 .. 1, numeric)) end proc

> 

A(3)

0.4547647069e-1

> 

A := proc (m) options operator, arrow; evalf((Int(xi*(-xi^2+1)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1))/(Int(xi*BesselJ(0, BesselJZeros(0, m)*xi)^2, xi = 0 .. 1))) end proc

proc (m) options operator, arrow; evalf((Int(xi*(-xi^2+1)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1))/(Int(xi*BesselJ(0, BesselJZeros(0, m)*xi)^2, xi = 0 .. 1))) end proc

> 

A(3)

0.4547647069e-1

Save an integral - see DLMF

> 

A := m -> 2*int(xi*(1 - xi^2)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric)/
          evalf(BesselJ(1, BesselJZeros(0, m))^2);

proc (m) options operator, arrow; 2*(int(xi*(-xi^2+1)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric))/evalf(BesselJ(1, BesselJZeros(0, m))^2) end proc

> 

A(3)

0.4547647070e-1

> 

NULL

Download Bessel.mw

BesselJZeros...

dharr 8200
April 27 2024
4 1 
evalf(BesselJZeros(0, 1..5));

gives: 2.404825558, 5.520078110, 8.653727913, 11.79153444, 14.93091771

Minpoly...

dharr 8200
April 26 2024
0 4

Here's how I would do it. Not sure why I get the two different answers.

> 

restart
> 

interface(version)

`Standard Worksheet Interface, Maple 2024.0, Windows 11, March 01 2024 Build ID 1794891`

> 

Digits := 40

40

> 

use alpha=1-(1/2)/(1-(RootOf(16*_Z*(_Z*(2*_Z*(_Z*(8*_Z*(_Z*(_Z*(_Z*(32*_Z*(8*_Z-33)+1513)-812)-13)+267)-1469)-330)+811)+279)+345,index=2)-1/2)**2) in
        expr:=(1+alpha)*sqrt(1-alpha**2)+(3+4*alpha)/12*sqrt(3-4*alpha**2)+2*(1+alpha)/3*sqrt(2*(1+alpha)*(1-2*alpha))+(1+2*alpha)/6*sqrt(2*((1-alpha)**2-3*alpha**2))
end:
> 

evalf(expr)

2.912910414540209166042800910083320363837

Looks correct

> 

PolynomialTools:-MinimalPolynomial(expr); fsolve(%)[2]

847155746562*_X^2-2718665563253*_X+731072836523

2.912910414540209166042800910083320363822

> 

alias(r__2 = RootOf(16*_Z*(_Z*(2*_Z*(_Z*(8*_Z*(_Z*(_Z*(_Z*(32*_Z*(8*_Z-33)+1513)-812)-13)+267)-1469)-330)+811)+279)+345, index = 2))

r__2

> 

expr2 := factor(evala(expr))

-(1/427776)*(-4*r__2^2+4*r__2+2)^(1/2)*(6307840*r__2^9+25065472*r__2^8-115353024*r__2^7+132632192*r__2^6-64357200*r__2^5+10339760*r__2^4+18097724*r__2^3-8262480*r__2^2-4524419*r__2-370681)

Square roots problematic, so convert to RootOf

> 

expr3 := convert(expr2, RootOf)

-(1/427776)*RootOf(_Z^2+4*r__2^2-4*r__2-2, index = 1)*(6307840*r__2^9+25065472*r__2^8-115353024*r__2^7+132632192*r__2^6-64357200*r__2^5+10339760*r__2^4+18097724*r__2^3-8262480*r__2^2-4524419*r__2-370681)

We want the minimum polynomial. Over the rationals - not sure why this is not the same as above

> 

evala(Minpoly(expr3, x)); fsolve(%)[-1]

x^20-(1099/72)*x^18+(1575661/20736)*x^16-(22165835/124416)*x^14+(40761112873/214990848)*x^12-(1002499649777/15479341056)*x^10+(11498404615715/1486016741376)*x^8-(3141147731417/17832200896512)*x^6-(86762801956691/184884258895036416)*x^4+(460470261157/164341563462254592)*x^2+33252616609/584325558976905216

2.912910414540209166042800910083320363836

With r__2as a field extension it is a quadratic

> 

evala(Minpoly(expr3, x, r__2)); fsolve(%)[2]

(8159503/855552)*r__2-(236960/1671)*r__2^9+(2085340/5013)*r__2^8-(13276063/40104)*r__2^7+(414451/20052)*r__2^6+(10073755/160416)*r__2^5-(13077487/160416)*r__2^4+(1471997/213888)*r__2^3+(5380211/160416)*r__2^2+1770613/2566656+x^2

2.912910414540209166042800910083320363854

> 

NULL

Download minpoly3.mw

First 17 18 19 20 21 22 23 Last Page 19 of 81