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

complexplot3d...

dharr 8200
August 01 2024
1 1

I've added missing information with random values. This gives a magnitude plot colored by phase.

ps - if providing code, please don't include the output.

> 

restart;
> 

local gamma;

gamma

Added random values for the parameters theta and gamma

> 

params := {alpha = 2.5, k = 3, w = 2, beta[3] = 3, beta[4] = 1.7, theta = 0, gamma = 0.2};
xi := beta[3]*(2*alpha*k*t + x)*sqrt(1/(36*alpha*beta[4] + 36*gamma*beta[4]));

{alpha = 2.5, gamma = .2, k = 3, theta = 0, w = 2, beta[3] = 3, beta[4] = 1.7}

beta[3]*(2*alpha*k*t+x)*(1/(36*alpha*beta[4]+36*gamma*beta[4]))^(1/2)

We perhaps want to plot the magnitude and color it with the phase - so make [mag, phase]

> 

sol1 := [U(xi), (-sqrt(1/(36*alpha*beta[4] + 36*gamma*beta[4]))*x + w*t + theta)];

[U(beta[3]*(2*alpha*k*t+x)*(1/(36*alpha*beta[4]+36*gamma*beta[4]))^(1/2)), -(1/(36*alpha*beta[4]+36*gamma*beta[4]))^(1/2)*x+w*t+theta]

Assume some functional form for U(xi)

> 

sol2 := eval(sol1, U(xi) = xi^2);

[beta[3]^2*(2*alpha*k*t+x)^2/(36*alpha*beta[4]+36*gamma*beta[4]), -(1/(36*alpha*beta[4]+36*gamma*beta[4]))^(1/2)*x+w*t+theta]

insert numerical values

> 

solnum :=eval(sol2, params);

[0.5446623094e-1*(15.0*t+x)^2, -0.7779333800e-1*x+2*t]

> 

plots:-complexplot3d(solnum, x = 0.. 5, t = 0..5);

> 

NULL

Download complexplot3d.mw

 

guess...

dharr 8200
July 31 2024
2 0

If you try

for y in [0, 0.5, 1, 5, 10] do
  p := plot(y);
  plottools:-getdata(p, rangesonly);
end do;

Then the y range is always reported just y..y, but the actual ranges are

y=0: -1..1

other y:   y-y/2..y+y/2

so perhaps that is a general rule.

The plot,ranges help just says "If the vertical range is unspecified, a suitable range will be deduced from the values computed" :-(

dsolve...

dharr 8200
July 30 2024
1 0

dsolve works here, and finds a simple correct solution. I'm guessing something about the problem is not correctly formulated here.

> 

restart
> 

with(Physics)
> 

Setup(realobjects = {a, g, v, K(diff(rho(x), x)), f(diff(rho(x), x))})

[realobjects = {a, g, v, K(`rho'`), f(`rho'`)}]

Original divided by nu^2

> 

y_1 := -`rho'`^2*(diff(K(`rho'`), `rho'`, `rho'`))-2*`rho'`*(diff(K(`rho'`), `rho'`))-K(`rho'`)*(K(`rho'`)^2*a*`rho'`^2*v^2+(1/2)*f(`rho'`)^2+(-2*`rho'`-1)*f(`rho'`)-a*v^2*`rho'`^2+2*`rho'`-3/2)

-`rho'`^2*(diff(diff(K(`rho'`), `rho'`), `rho'`))-2*(diff(K(`rho'`), `rho'`))*`rho'`-K(`rho'`)*(K(`rho'`)^2*a*`rho'`^2*v^2+(1/2)*f(`rho'`)^2+(-2*`rho'`-1)*f(`rho'`)-a*v^2*`rho'`^2+2*`rho'`-3/2)

> 

y_2 := numer(normal(2*(diff(f(`rho'`), `rho'`, `rho'`))/g^2+(-4*f(`rho'`)^3+24*f(`rho'`)^2+(-`rho'`^2*v^2*g^2*K(`rho'`)^2-44)*f(`rho'`)+24+(2*(`rho'`+1/2))*g^2*`rho'`^2*v^2*K(`rho'`)^2)/(2*`rho'`^2*g^2)))

-f(`rho'`)*K(`rho'`)^2*g^2*`rho'`^2*v^2+2*`rho'`^3*v^2*g^2*K(`rho'`)^2+`rho'`^2*v^2*g^2*K(`rho'`)^2+4*(diff(diff(f(`rho'`), `rho'`), `rho'`))*`rho'`^2-4*f(`rho'`)^3+24*f(`rho'`)^2-44*f(`rho'`)+24

> 

ics := f(0) = 1, K(0) = 0, (D(f))(inf) = 0, (D(K))(inf) = 0

f(0) = 1, K(0) = 0, (D(f))(inf) = 0, (D(K))(inf) = 0

> 

params := {a = 2, g = 1.5, inf = 20, v = 3}

{a = 2, g = 1.5, inf = 20, v = 3}

> 

sol := dsolve(eval({ics, y_1, y_2}, params), numeric)

Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead

> 

sol := dsolve(eval({ics, y_1, y_2}, params), numeric, method = bvp[midrich])

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = 2.8571428571428505, (3) = 5.714285714285708, (4) = 8.57142857142857, (5) = 11.428571428571427, (6) = 14.285714285714283, (7) = 17.142857142857142, (8) = 20.0}, datatype = float[8], order = C_order); Y := Matrix(8, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 1.0, (1, 4) = 0.17931463541076747e-18, (2, 1) = .0, (2, 2) = .0, (2, 3) = 1.0, (2, 4) = -0.5379439062323024e-18, (3, 1) = .0, (3, 2) = .0, (3, 3) = 1.0, (3, 4) = -0.1050271435977354e-17, (4, 1) = .0, (4, 2) = .0, (4, 3) = 1.0, (4, 4) = -0.15715090097179726e-17, (5, 1) = .0, (5, 2) = .0, (5, 3) = 1.0, (5, 4) = -0.20940736112877165e-17, (6, 1) = .0, (6, 2) = .0, (6, 3) = 1.0, (6, 4) = -0.64954072162389115e-18, (7, 1) = .0, (7, 2) = .0, (7, 3) = 1.0, (7, 4) = 0.26576829106897595e-18, (8, 1) = .0, (8, 2) = .0, (8, 3) = 1.0, (8, 4) = .0}, datatype = float[8], order = C_order); YP := Matrix(8, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = -.0, (2, 3) = -0.5379439062323024e-18, (2, 4) = .0, (3, 1) = .0, (3, 2) = -.0, (3, 3) = -0.1050271435977354e-17, (3, 4) = .0, (4, 1) = .0, (4, 2) = -.0, (4, 3) = -0.15715090097179726e-17, (4, 4) = .0, (5, 1) = .0, (5, 2) = -.0, (5, 3) = -0.20940736112877165e-17, (5, 4) = .0, (6, 1) = .0, (6, 2) = -.0, (6, 3) = -0.64954072162389115e-18, (6, 4) = .0, (7, 1) = .0, (7, 2) = -.0, (7, 3) = 0.26576829106897595e-18, (7, 4) = .0, (8, 1) = .0, (8, 2) = -.0, (8, 3) = .0, (8, 4) = .0}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = 2.8571428571428505, (3) = 5.714285714285708, (4) = 8.57142857142857, (5) = 11.428571428571427, (6) = 14.285714285714283, (7) = 17.142857142857142, (8) = 20.0}, datatype = float[8], order = C_order); Y := Matrix(8, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = -0.2390861805476899e-18, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = 0.7172585416430703e-18, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = 0.14003619146364716e-17, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = 0.2095345346290631e-17, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = 0.2792098148383623e-17, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = 0.86605429549852e-18, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = -0.3543577214253028e-18, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(2.792098148383623e-18) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 8, [K(`rho'`), diff(K(`rho'`), `rho'`), f(`rho'`), diff(f(`rho'`), `rho'`)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 4, X, Y, outpoint, yout, L, V) end if; [`rho'` = outpoint, seq('[K(`rho'`), diff(K(`rho'`), `rho'`), f(`rho'`), diff(f(`rho'`), `rho'`)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(2.792098148383623e-18) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (true), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (true), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(8, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [`rho'`, K(`rho'`), diff(K(`rho'`), `rho'`), f(`rho'`), diff(f(`rho'`), `rho'`)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [`rho'` = res[1], seq('[K(`rho'`), diff(K(`rho'`), `rho'`), f(`rho'`), diff(f(`rho'`), `rho'`)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc

> 

plots:-odeplot(sol, [[diff(rho(x), x), K(`rho'`)], [diff(rho(x), x), f(`rho'`)]], 0 .. eval(inf, params), color = [red, blue])

Indeed K=0 and f=1 is a solution.

> 

odetest([K(diff(rho(x), x)) = 0, f(diff(rho(x), x)) = 1], [y_1, y_2])

[0, 0]

> 

NULL

Download dsolve.mw

some suggestions...

dharr 8200
July 29 2024
0 0

You have a large polynomial system that is going to be slow. Among other changes, removing the explicit option saves about 25%. I have also added general suggestions that may be useful in the future.

new1.mw

select...

dharr 8200
July 29 2024
1 4

Your worksheet does not run so I copied the output.

> 

restart
> 

sol := {k = k, w = w, A[0] = 0, A[1] = A[1], B[1] = 0, a[1] = a[1], a[2] = a[2], a[3] = a[3], a[4] = a[4], a[5] = a[5]}, {k = k, w = -4*A[0]*a[5], A[0] = A[0], A[1] = A[1], B[1] = B[1], a[1] = 0, a[2] = -4*a[5], a[3] = 0, a[4] = 0, a[5] = a[5]}, {k = k, w = (1/2)*A[0]*(3*k^2*A[0]^2*a[4]+2*k^2*A[0]*a[3]+k^2*a[2]+4*k^2*a[5]+2*A[0]^2*a[4]+2*A[0]*a[3]+2*a[2]), A[0] = A[0], A[1] = 0, B[1] = 0, a[1] = -(1/2)*A[0]*(3*A[0]^2*a[4]+2*A[0]*a[3]+a[2]+4*a[5]), a[2] = a[2], a[3] = a[3], a[4] = a[4], a[5] = a[5]}, {k = k, w = w, A[0] = A[0], A[1] = 0, B[1] = 0, a[1] = a[1], a[2] = (-A[0]^3*a[4]+k^2*a[1]-A[0]^2*a[3]+w)/A[0], a[3] = a[3], a[4] = a[4], a[5] = a[5]}, {k = k, w = 4*A[1]*a[5]+4*B[1]*a[5], A[0] = -A[1]-B[1], A[1] = A[1], B[1] = B[1], a[1] = 0, a[2] = -4*a[5], a[3] = 0, a[4] = 0, a[5] = a[5]}, {k = k, w = -k^2*a[1]-4*A[0]*a[5]+a[1], A[0] = A[0], A[1] = A[0]^2/(4*B[1]), B[1] = B[1], a[1] = a[1], a[2] = -4*a[5], a[3] = 0, a[4] = 0, a[5] = a[5]}, {k = k, w = w, A[0] = 2*B[1], A[1] = B[1], B[1] = B[1], a[1] = a[1], a[2] = (k^2*a[1]+w-a[1])/(2*B[1]), a[3] = 0, a[4] = 0, a[5] = -(k^2*a[1]+w-a[1])/(8*B[1])}, {k = k, w = w, A[0] = A[0], A[1] = B[1], B[1] = B[1], a[1] = 0, a[2] = w/A[0], a[3] = 0, a[4] = 0, a[5] = -w/(4*A[0])}, {k = k, w = 0, A[0] = 0, A[1] = B[1], B[1] = B[1], a[1] = 0, a[2] = a[2], a[3] = 0, a[4] = 0, a[5] = -(1/4)*a[2]}

You can select by number

> 

sol[3]

{k = k, w = (1/2)*A[0]*(3*k^2*A[0]^2*a[4]+2*k^2*A[0]*a[3]+k^2*a[2]+4*k^2*a[5]+2*A[0]^2*a[4]+2*A[0]*a[3]+2*a[2]), A[0] = A[0], A[1] = 0, B[1] = 0, a[1] = -(1/2)*A[0]*(3*A[0]^2*a[4]+2*A[0]*a[3]+a[2]+4*a[5]), a[2] = a[2], a[3] = a[3], a[4] = a[4], a[5] = a[5]}

select(has, .., is useful here. Note the [] around sol.

> 

select(has, [sol], a[1] = a[1])

[{k = k, w = w, A[0] = 0, A[1] = A[1], B[1] = 0, a[1] = a[1], a[2] = a[2], a[3] = a[3], a[4] = a[4], a[5] = a[5]}, {k = k, w = w, A[0] = A[0], A[1] = 0, B[1] = 0, a[1] = a[1], a[2] = (-A[0]^3*a[4]+k^2*a[1]-A[0]^2*a[3]+w)/A[0], a[3] = a[3], a[4] = a[4], a[5] = a[5]}, {k = k, w = -k^2*a[1]-4*A[0]*a[5]+a[1], A[0] = A[0], A[1] = (1/4)*A[0]^2/B[1], B[1] = B[1], a[1] = a[1], a[2] = -4*a[5], a[3] = 0, a[4] = 0, a[5] = a[5]}, {k = k, w = w, A[0] = 2*B[1], A[1] = B[1], B[1] = B[1], a[1] = a[1], a[2] = (1/2)*(k^2*a[1]+w-a[1])/B[1], a[3] = 0, a[4] = 0, a[5] = -(1/8)*(k^2*a[1]+w-a[1])/B[1]}]

remove or selectremove may also be useful.

> 

NULL

Download choosecase.mw

convert, polynom...

dharr 8200
July 26 2024
1 1

The problem is the O(h^7) term. If q is your expand(SUMY-T) then

coeff(convert(q, polynom), h, 6)

(actually you should remove it from SUMY)

difficult...

dharr 8200
July 26 2024
1 0

I don't see any way in the general case. But Z2 is simple enough to do. The comments suggest why it will be difficult in the general case. Zp is the cyclic group of order p.

> 

restart
> 

with(GroupTheory); with(GraphTheory)
> 

p := 11

11

Construct direct product Zp x Zp. Goes to Z x Z in limit of large p.

> 

Zp2 := `<|>`(`<,>`(a, b), `<,>`(a^p = 1, b^p = 1, 1/b.(1/a).b.a))

_m2184784207904

> 

cg := CayleyGraph(Zp2)

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121], Array(1..121, {(1) = {2, 3}, (2) = {5, 6}, (3) = {6, 7}, (4) = {1, 8}, (5) = {10, 11}, (6) = {11, 12}, (7) = {12, 13}, (8) = {2, 14}, (9) = {4, 15}, (10) = {16, 17}, (11) = {17, 18}, (12) = {18, 19}, (13) = {19, 20}, (14) = {5, 21}, (15) = {1, 22}, (16) = {23, 24}, (17) = {24, 25}, (18) = {25, 26}, (19) = {26, 27}, (20) = {27, 28}, (21) = {10, 29}, (22) = {3, 30}, (23) = {31, 32}, (24) = {32, 33}, (25) = {33, 34}, (26) = {34, 35}, (27) = {35, 36}, (28) = {36, 37}, (29) = {16, 38}, (30) = {7, 39}, (31) = {40, 41}, (32) = {41, 42}, (33) = {42, 43}, (34) = {43, 44}, (35) = {44, 45}, (36) = {45, 46}, (37) = {46, 47}, (38) = {23, 48}, (39) = {13, 49}, (40) = {50, 51}, (41) = {51, 52}, (42) = {52, 53}, (43) = {53, 54}, (44) = {54, 55}, (45) = {55, 56}, (46) = {56, 57}, (47) = {57, 58}, (48) = {31, 59}, (49) = {20, 60}, (50) = {61, 62}, (51) = {62, 63}, (52) = {63, 64}, (53) = {64, 65}, (54) = {65, 66}, (55) = {66, 67}, (56) = {67, 68}, (57) = {68, 69}, (58) = {69, 70}, (59) = {40, 71}, (60) = {28, 72}, (61) = {15, 73}, (62) = {73, 74}, (63) = {74, 75}, (64) = {75, 76}, (65) = {76, 77}, (66) = {77, 78}, (67) = {78, 79}, (68) = {79, 80}, (69) = {80, 81}, (70) = {4, 81}, (71) = {50, 83}, (72) = {37, 84}, (73) = {22, 85}, (74) = {85, 86}, (75) = {86, 87}, (76) = {87, 88}, (77) = {88, 89}, (78) = {89, 90}, (79) = {90, 91}, (80) = {91, 92}, (81) = {8, 92}, (82) = {9, 70}, (83) = {9, 61}, (84) = {47, 94}, (85) = {30, 95}, (86) = {95, 96}, (87) = {96, 97}, (88) = {97, 98}, (89) = {98, 99}, (90) = {99, 100}, (91) = {100, 101}, (92) = {14, 101}, (93) = {82, 83}, (94) = {58, 82}, (95) = {39, 102}, (96) = {102, 103}, (97) = {103, 104}, (98) = {104, 105}, (99) = {105, 106}, (100) = {106, 107}, (101) = {21, 107}, (102) = {49, 108}, (103) = {108, 109}, (104) = {109, 110}, (105) = {110, 111}, (106) = {111, 112}, (107) = {29, 112}, (108) = {60, 113}, (109) = {113, 114}, (110) = {114, 115}, (111) = {115, 116}, (112) = {38, 116}, (113) = {72, 117}, (114) = {117, 118}, (115) = {118, 119}, (116) = {48, 119}, (117) = {84, 120}, (118) = {120, 121}, (119) = {59, 121}, (120) = {93, 94}, (121) = {71, 93}}), `GRAPHLN/table/1`, 0)

Torus as expected

> 

DrawGraph(cg, dimension = 3, orientation = [34, 50, 35])

The problem is that the order of the elements is not useful, and the representation is not unique, e.g., b^2 can be b^3.b^-1. So we need some manipulations that are easy here but wont be in the general case.

> 

els := [Elements(Zp2)[]]

Find the orders of a and b

> 

degs := proc (els) local a, b, prods; a, b := indets(els)[]; prods := map(mul, els); map(proc (x) options operator, arrow; [`mod`(degree(x, a), p), `mod`(degree(x, b), p)] end proc, prods) end proc
> 

d := degs(els)

[[0, 0], [1, 0], [0, 1], [0, 10], [2, 0], [1, 1], [0, 2], [1, 10], [10, 10], [3, 0], [2, 1], [1, 2], [0, 3], [2, 10], [10, 0], [4, 0], [3, 1], [2, 2], [1, 3], [0, 4], [3, 10], [10, 1], [5, 0], [4, 1], [3, 2], [2, 3], [1, 4], [0, 5], [4, 10], [10, 2], [6, 0], [5, 1], [4, 2], [3, 3], [2, 4], [1, 5], [0, 6], [5, 10], [10, 3], [7, 0], [6, 1], [5, 2], [4, 3], [3, 4], [2, 5], [1, 6], [0, 7], [6, 10], [10, 4], [8, 0], [7, 1], [6, 2], [5, 3], [4, 4], [3, 5], [2, 6], [1, 7], [0, 8], [7, 10], [10, 5], [9, 0], [8, 1], [7, 2], [6, 3], [5, 4], [4, 5], [3, 6], [2, 7], [1, 8], [0, 9], [8, 10], [10, 6], [9, 1], [8, 2], [7, 3], [6, 4], [5, 5], [4, 6], [3, 7], [2, 8], [1, 9], [10, 9], [9, 10], [10, 7], [9, 2], [8, 3], [7, 4], [6, 5], [5, 6], [4, 7], [3, 8], [2, 9], [9, 9], [10, 8], [9, 3], [8, 4], [7, 5], [6, 6], [5, 7], [4, 8], [3, 9], [9, 4], [8, 5], [7, 6], [6, 7], [5, 8], [4, 9], [9, 5], [8, 6], [7, 7], [6, 8], [5, 9], [9, 6], [8, 7], [7, 8], [6, 9], [9, 7], [8, 8], [7, 9], [9, 8], [8, 9]]

We want to split between degrees 5 and 6. Make an Array of the element numbers indexed by the degrees

> 

M := Array(0 .. p-1, 0 .. p-1, {map(proc (x) options operator, arrow; op(lhs(x)) = rhs(x) end proc, `~`[`=`](d, [`$`(1 .. nops(d))]))[]})

_rtable[36893490332336398692]

Delete arcs between 5 and 6. In this case we know the arc directions, but in the more general case this won't be true.

> 

start := (p-1)*(1/2); fin := start+1

5

6

> 

DeleteArc(cg, {convert(M[() .. (), start .. fin], listlist)[]}); DeleteArc(cg, {convert(M*LinearAlgebra:-Transpose([start .. fin, () .. ()]), listlist)[]})

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121], Array(1..121, {(1) = {2, 3}, (2) = {5, 6}, (3) = {6, 7}, (4) = {1, 8}, (5) = {10, 11}, (6) = {11, 12}, (7) = {12, 13}, (8) = {2, 14}, (9) = {4, 15}, (10) = {16, 17}, (11) = {17, 18}, (12) = {18, 19}, (13) = {19, 20}, (14) = {5, 21}, (15) = {1, 22}, (16) = {23, 24}, (17) = {24, 25}, (18) = {25, 26}, (19) = {26, 27}, (20) = {27, 28}, (21) = {10, 29}, (22) = {3, 30}, (23) = {31, 32}, (24) = {32, 33}, (25) = {33, 34}, (26) = {34, 35}, (27) = {35, 36}, (28) = {36}, (29) = {16, 38}, (30) = {7, 39}, (31) = {40, 41}, (32) = {41, 42}, (33) = {42, 43}, (34) = {43, 44}, (35) = {44, 45}, (36) = {45}, (37) = {46, 47}, (38) = {23, 48}, (39) = {13, 49}, (40) = {50, 51}, (41) = {51, 52}, (42) = {52, 53}, (43) = {53, 54}, (44) = {54, 55}, (45) = {55}, (46) = {56, 57}, (47) = {57, 58}, (48) = {31, 59}, (49) = {20, 60}, (50) = {61, 62}, (51) = {62, 63}, (52) = {63, 64}, (53) = {64, 65}, (54) = {65, 66}, (55) = {66}, (56) = {67, 68}, (57) = {68, 69}, (58) = {69, 70}, (59) = {40, 71}, (60) = {28}, (61) = {15, 73}, (62) = {73, 74}, (63) = {74, 75}, (64) = {75, 76}, (65) = {76, 77}, (66) = {77}, (67) = {78, 79}, (68) = {79, 80}, (69) = {80, 81}, (70) = {4, 81}, (71) = {50, 83}, (72) = {37, 84}, (73) = {22, 85}, (74) = {85, 86}, (75) = {86, 87}, (76) = {87, 88}, (77) = {88}, (78) = {89, 90}, (79) = {90, 91}, (80) = {91, 92}, (81) = {8, 92}, (82) = {9, 70}, (83) = {9, 61}, (84) = {47, 94}, (85) = {30, 95}, (86) = {95, 96}, (87) = {96, 97}, (88) = {97}, (89) = {98, 99}, (90) = {99, 100}, (91) = {100, 101}, (92) = {14, 101}, (93) = {82, 83}, (94) = {58, 82}, (95) = {39, 102}, (96) = {102, 103}, (97) = {103}, (98) = {104, 105}, (99) = {105, 106}, (100) = {106, 107}, (101) = {21, 107}, (102) = {49, 108}, (103) = {108}, (104) = {109, 110}, (105) = {110, 111}, (106) = {111, 112}, (107) = {29, 112}, (108) = {60}, (109) = {113, 114}, (110) = {114, 115}, (111) = {115, 116}, (112) = {38, 116}, (113) = {72, 117}, (114) = {117, 118}, (115) = {118, 119}, (116) = {48, 119}, (117) = {84, 120}, (118) = {120, 121}, (119) = {59, 121}, (120) = {93, 94}, (121) = {71, 93}}), `GRAPHLN/table/1`, 0)

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121], Array(1..121, {(1) = {2, 3}, (2) = {5, 6}, (3) = {6, 7}, (4) = {1, 8}, (5) = {10, 11}, (6) = {11, 12}, (7) = {12, 13}, (8) = {2, 14}, (9) = {4, 15}, (10) = {16, 17}, (11) = {17, 18}, (12) = {18, 19}, (13) = {19, 20}, (14) = {5, 21}, (15) = {1, 22}, (16) = {23, 24}, (17) = {24, 25}, (18) = {25, 26}, (19) = {26, 27}, (20) = {27, 28}, (21) = {10, 29}, (22) = {3, 30}, (23) = {32}, (24) = {32, 33}, (25) = {33, 34}, (26) = {34, 35}, (27) = {35, 36}, (28) = {36}, (29) = {16, 38}, (30) = {7, 39}, (31) = {40, 41}, (32) = {42}, (33) = {42, 43}, (34) = {43, 44}, (35) = {44, 45}, (36) = {45}, (37) = {46, 47}, (38) = {23}, (39) = {13, 49}, (40) = {50, 51}, (41) = {51, 52}, (42) = {53}, (43) = {53, 54}, (44) = {54, 55}, (45) = {55}, (46) = {56, 57}, (47) = {57, 58}, (48) = {31, 59}, (49) = {20, 60}, (50) = {61, 62}, (51) = {62, 63}, (52) = {63, 64}, (53) = {65}, (54) = {65, 66}, (55) = {66}, (56) = {67, 68}, (57) = {68, 69}, (58) = {69, 70}, (59) = {40, 71}, (60) = {28}, (61) = {15, 73}, (62) = {73, 74}, (63) = {74, 75}, (64) = {75, 76}, (65) = {77}, (66) = {77}, (67) = {78, 79}, (68) = {79, 80}, (69) = {80, 81}, (70) = {4, 81}, (71) = {50, 83}, (72) = {37, 84}, (73) = {22, 85}, (74) = {85, 86}, (75) = {86, 87}, (76) = {87, 88}, (77) = {}, (78) = {89, 90}, (79) = {90, 91}, (80) = {91, 92}, (81) = {8, 92}, (82) = {9, 70}, (83) = {9, 61}, (84) = {47, 94}, (85) = {30, 95}, (86) = {95, 96}, (87) = {96, 97}, (88) = {97}, (89) = {99}, (90) = {99, 100}, (91) = {100, 101}, (92) = {14, 101}, (93) = {82, 83}, (94) = {58, 82}, (95) = {39, 102}, (96) = {102, 103}, (97) = {103}, (98) = {104, 105}, (99) = {106}, (100) = {106, 107}, (101) = {21, 107}, (102) = {49, 108}, (103) = {108}, (104) = {109, 110}, (105) = {110, 111}, (106) = {112}, (107) = {29, 112}, (108) = {60}, (109) = {113, 114}, (110) = {114, 115}, (111) = {115, 116}, (112) = {38}, (113) = {72, 117}, (114) = {117, 118}, (115) = {118, 119}, (116) = {48, 119}, (117) = {84, 120}, (118) = {120, 121}, (119) = {59, 121}, (120) = {93, 94}, (121) = {71, 93}}), `GRAPHLN/table/1`, 0)

And we get part of an infinite grid, as expected.

> 

DrawGraph(cg, layout = spring)

> 

NULL

Download CayleyGraphZ2.mw

showstat...

dharr 8200
July 26 2024
2 0

Try

showstat(Units:-TestDimensions);

 

assumptions with invlaplace...

dharr 8200
July 24 2024
1 3

invlaplace very often needs assumptions to make progress. Even though s is complex, assuming everything is positive often seems to work. I just worked out lap3, since I am out of time, but this method is probably enough to get the rest done.

> 

NULL
> 

restart
> 

with(inttrans)
> 

pde := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = t^2*x+x

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = t^2*x+x

Work with the three parts separately

> 

pt1, pt2, pt3 := diff(u(x, t), t), u(x, t)*(diff(u(x, t), x)), t^2*x+x

diff(u(x, t), t), u(x, t)*(diff(u(x, t), x)), t^2*x+x

part 1 by hand. L1 is invlaplace of (laplace divided by s^alpha)

> 

L1 := u(x, t)

u(x, t)

> 

L2 := invlaplace(laplace(pt2, t, s)/s^alpha, s, t)

invlaplace(laplace(u(x, t)*(diff(u(x, t), x)), t, s)/s^alpha, s, t)

Part 3 requires some assumptions that work even if they aren't correct

> 

L3 := `assuming`([invlaplace(laplace(pt3, t, s)/s^alpha, s, t)], [positive])

x*t^alpha*(1/GAMMA(1+alpha)+2*t^2/GAMMA(3+alpha))

> 

lap3 := L1 = L3-L2

u(x, t) = x*t^alpha*(1/GAMMA(1+alpha)+2*t^2/GAMMA(3+alpha))-invlaplace(laplace(u(x, t)*(diff(u(x, t), x)), t, s)/s^alpha, s, t)

> 

``

Download lap3.mw

unknown function...

dharr 8200
July 24 2024
2 3

Your proposed solution (9) contains the unknown function R(xi), so I'm not sure what you expect here. You can use

simplify(eval(SO, f))

to see the resulting ode in R(xi). I'm not sure if this is the question you are asking.

Minimize...

dharr 8200
July 24 2024
2 1

Since you are using floating point numbers, an alternative is

Optimization:-Minimize(G(x,y), x = 0..1, y = 0..1)

In principle it might find a local rather than global minimum, but a plot shows it is fine in this case.

Array, then Matrix...

dharr 8200
July 22 2024
3 2

@Scot Gould 's solution is nice (Vote up), but it is possible to put the Vectors directly in the Matrix if you first make an Array.

> 

restart
> 

vectors1 := [seq(`<,>`(i, i+1, i+2), i = 1 .. 3)]; vectors2 := [seq(`<,>`(i, 2*i, 3*i), i = 1 .. 3)]

[Vector(3, {(1) = 1, (2) = 2, (3) = 3}), Vector(3, {(1) = 2, (2) = 3, (3) = 4}), Vector(3, {(1) = 3, (2) = 4, (3) = 5})]

[Vector[column](%id = 36893491208837358164), Vector[column](%id = 36893491208837358284), Vector[column](%id = 36893491208837358404)]

Arrays accept Vectors as entries

> 

arrowList := Array(1 .. 3, 1 .. 2, [seq([vectors1[n], vectors2[n]], n = 1 .. 3)])

Array(%id = 36893491208837351772)

And we can convert that to a Matrix

> 

arrowList := Matrix(arrowList)

Matrix(%id = 36893491208837333588)

> 

arrowList[3, 2]

Vector[column](%id = 36893491208837358404)

> 

arrowList[3, 2][1]

3

> 

NULL

Download MaplePrimes_Matrix_of_Vectors.mw

psi assigned...

dharr 8200
July 22 2024
2 5

The main error here was that psi(x,t) was already assigned when you did pdetest. I don't know the assumptions about the signs of things, but to make progress you will need to make them explicit. I think the presence of the absolute values in pde is going to get you into trouble if you are using complex functions. You had an earlier ode in U(xi) that you derived from the pde. So I would make sure that odetest on that works before attempting the full pde. 

You had kt which I changed to k*t but in your notes it is just k. Isn't all the time dependence in the exponential?

pde-solve.mw

solve,identity...

dharr 8200
July 20 2024
2 1

solve, identity will automatically do all the work for you

solve(identity(E1, xi), {k, lambda, w, A[0], A[1], B[1]})

See the attached worksheet for the results.

solve_identity.mw

pseudo-inverse...

dharr 8200
July 18 2024
3 1

The pseudo inverse of A seems to work here, but off-hand I can't justify this.

Groebner.mw

 

First 15 16 17 18 19 20 21 Last Page 17 of 81