Advanced Search

dharr

Dr. David Harrington

8195 Reputation

22 Badges

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

union is for sets...

dharr 8195
May 02 2025
2 1

The error message is clear: "union received [...". Union applies to sets, so changing case1 from a list to a set solves your problem.

ode-not-run.mw

Get a better book...

dharr 8195
May 01 2025
2 5

Some of these work for restricted regions, but others are wrong. I didn't work out all the form 2 ones; similar procedures apply.

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

> 

_local(gamma)
> 

declare(G(xi))

G(xi)*`will now be displayed as`*G

> 

NULL

Eq7.

> 

ode := diff(G(xi), xi) = G(xi)*sqrt(A+B*G(xi)^2)

diff(G(xi), xi) = G(xi)*(A+B*G(xi)^2)^(1/2)

Eq 8 for A,B>0. Changed xi[0] to xi__0

> 

S1 := G(xi) = -sqrt(A/B)*csch(sqrt(A)*(xi+`ξ__0`))

G(xi) = -(A/B)^(1/2)*csch(A^(1/2)*(xi+xi__0))

This one isn't zero, so solution is not correct.

> 

res := `assuming`([simplify(odetest(S1, ode))], [A > 0, B > 0, sqrt(A)*(xi+`ξ__0`) > 0])

2*A*csch(A^(1/2)*(xi+xi__0))*coth(A^(1/2)*(xi+xi__0))/B^(1/2)

Only zero at infinity

> 

plot(csch(x)*coth(x))

Eq 9 for A<0,B>0

> 

S2 := G(xi) = sqrt(-A/B)*sec(sqrt(-A)*(xi+`&xi;__0`))

G(xi) = (-A/B)^(1/2)*sec((-A)^(1/2)*(xi+xi__0))

So this almost works.

> 

res := `assuming`([simplify(odetest(S2, ode))], [A < 0, B > 0])

A*sec((-A)^(1/2)*(xi+xi__0))*tan((-A)^(1/2)*(xi+xi__0))*(csgn(tan((-A)^(1/2)*(xi+xi__0)))-1)/B^(1/2)

One possibility is assuming tan(sqrt(-A)*(xi+`&xi;__0`)) > 0or just choose a range that makes that true

> 

res := `assuming`([simplify(odetest(S2, ode))], [A < 0, B > 0, sqrt(-A)*(xi+`&xi;__0`) > 0, sqrt(-A)*(xi+`&xi;__0`) < (1/2)*Pi])

0

So can we convert this to see what Eq 8 should have been?

> 

`assuming`([convert(S2, csch)], [A > 0, B > 0])

G(xi) = I*(-A/B)^(1/2)*csch(I*((-A)^(1/2)*(xi+xi__0)+(1/2)*Pi))

Eq 10 for A>0,B<0

> 

S3 := G(xi) = sqrt(-A/B)*sech(sqrt(A)*(xi+`&xi;__0`))

G(xi) = (-A/B)^(1/2)*sech(A^(1/2)*(xi+xi__0))

Find need sqrt(A)*(xi+`&xi;__0`) < 0- not realistic?

> 

res := `assuming`([simplify(odetest(S3, ode))], [A > 0, B < 0, sqrt(A)*(xi+`&xi;__0`) < 0])

0

Eq 11 for A<0,B>0

> 

S4 := G(xi) = sqrt(-A/B)*csc(sqrt(-A)*(xi+`&xi;__0`))

G(xi) = (-A/B)^(1/2)*csc((-A)^(1/2)*(xi+xi__0))

> 

res := `assuming`([simplify(odetest(S4, ode))], [A < 0, B > 0])

A*csc((-A)^(1/2)*(xi+xi__0))*cot((-A)^(1/2)*(xi+xi__0))*(csgn(cot((-A)^(1/2)*(xi+xi__0)))+1)/B^(1/2)

So

> 

res := `assuming`([simplify(odetest(S4, ode))], [A < 0, B > 0, cot(sqrt(-A)*(xi+`&xi;__0`)) < 0])

0

Eq 12 for A<0, B>0. Eq 12 has an I in.

> 

S5 := G(xi) = cos(sqrt(-A)*(xi+`&xi;__0`))+I*sin(sqrt(-A)*(xi+`&xi;__0`))

G(xi) = cos((-A)^(1/2)*(xi+xi__0))+I*sin((-A)^(1/2)*(xi+xi__0))

This one is wrong

> 

res := `assuming`([simplify(odetest(S5, ode))], [A < 0, B > 0]); eval(res, {`&xi;__0` = 1., A = -1, B = -2, xi = 1})

(cos((-A)^(1/2)*(xi+xi__0))+I*sin((-A)^(1/2)*(xi+xi__0)))*(-(I*B*sin(2*(-A)^(1/2)*(xi+xi__0))+B*cos(2*(-A)^(1/2)*(xi+xi__0))+A)^(1/2)+I*(-A)^(1/2))

.2063024051-.9637957172*I

Eq 13. for A=0,B>0

> 

S6 := G(xi) = 1/(sqrt(B)*(xi+`&xi;__0`))

G(xi) = 1/(B^(1/2)*(xi+xi__0))

This one is wrong.

> 

`assuming`([odetest(S6, subs(A = 0, ode))], [B > 0, xi+`&xi;__0` > 0])

-2/(B^(1/2)*(xi+xi__0)^2)

Eq 14. A=0, B<0

> 

S7 := G(xi) = 1/(sqrt(-B)*(xi+`&xi;__0`))

G(xi) = 1/((-B)^(1/2)*(xi+xi__0))

Also wrong

> 

`assuming`([odetest(S7, subs(A = 0, ode))], [B < 0, xi+`&xi;__0` > 0])

-1/((-B)^(1/2)*(xi+xi__0)^2)-I/((-B)^(1/2)*(xi+xi__0)^2)

> 

``

Download Z1.mw

Corrected...

dharr 8195
May 01 2025
3 1
> 

restart
> 

with(PDEtools)
> 

declare(P(mu))

P(mu)*`will now be displayed as`*P

> 

Corrected version A^2+r -> A^2

> 

Psol := P(mu) = 2*A*lambda/(A^2*exp(lambda*rho*mu)+r*exp(-lambda*rho*mu))

P(mu) = 2*A*lambda/(A^2*exp(lambda*rho*mu)+r*exp(-lambda*rho*mu))

Corrected for sign to correspond to above solution

> 

ode := (diff(P(mu), mu))^2-rho^2*P(mu)^2*(1-r*P(mu)^2) = 0

(diff(P(mu), mu))^2-rho^2*P(mu)^2*(1-r*P(mu)^2) = 0

Since lambda = +/-1 the following is zero

> 

res := simplify(odetest(Psol, ode)); eval(res, lambda = 1); eval(res, lambda = -1)

0

0

 

> 

P_hyper := P(mu) = 2*A*lambda/((A^2+r)*cosh(rho*mu)+(A^2-r)*sinh(rho*mu))

P(mu) = 2*A*lambda/((A^2+r)*cosh(rho*mu)+(A^2-r)*sinh(rho*mu))

> 

res_hyper := simplify(odetest(P_hyper, ode)); eval(res_hyper, lambda = 1); eval(res_hyper, lambda = -1)

16*rho^2*(lambda^4-lambda^2)*r*A^4/((A^2+r)*cosh(rho*mu)+(A^2-r)*sinh(rho*mu))^4

0

0

 

NULL

Download ode.m

subs...

dharr 8195
April 29 2025
1 7

subs is best for substituting name=expression. When you don't have a name on the left, it can still work if the left-hand side is an "indet" of the target. In this case it is not. The lhs has the sqrt, but the target has 1/sqrt and they are not the same.

Algsubs can sometimes solve this problem, but not when there are square roots. The simplest solution here is probably to take the reciprocal, do the substitution, and take the reciprocal of the answer. When doing substitutions involving derivatives, dsubs is preferred, though here it give the same answer.

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

> 

NULL
> 

S := (diff(G(xi), xi))^2-r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2) = 0

(diff(G(xi), xi))^2-r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2) = 0

> 

SS := diff(G(xi), xi) = sqrt(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))

diff(G(xi), xi) = (r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2)

> 

Se := sqrt(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2)) = diff(G(xi), xi)

(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2) = diff(G(xi), xi)

> 

dub := diff(SS, xi)

diff(diff(G(xi), xi), xi) = (1/2)*(2*r^2*G(xi)*(a+b*G(xi)+l*G(xi)^2)*(diff(G(xi), xi))+r^2*G(xi)^2*(b*(diff(G(xi), xi))+2*l*G(xi)*(diff(G(xi), xi))))/(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2)

> 

Dubl2 := simplify(diff(diff(G(xi), xi), xi) = (1/2)*(2*r^2*G(xi)*(a+b*G(xi)+l*G(xi)^2)*(diff(G(xi), xi))+r^2*G(xi)^2*(b*(diff(G(xi), xi))+2*l*G(xi)*(diff(G(xi), xi))))/(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2))

diff(diff(G(xi), xi), xi) = (1/2)*r^2*G(xi)*(diff(G(xi), xi))*(4*l*G(xi)^2+3*b*G(xi)+2*a)/(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2)

> 

subs(Se, Dubl2)

diff(diff(G(xi), xi), xi) = (1/2)*r^2*G(xi)*(diff(G(xi), xi))*(4*l*G(xi)^2+3*b*G(xi)+2*a)/(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2)

> 

indets(Se); indets(Dubl2)

{a, b, l, r, xi, (r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2), G(xi), diff(G(xi), xi)}

{a, b, l, r, xi, 1/(r^2*G(xi)^2*(a+b*G(xi)+l*G(xi)^2))^(1/2), G(xi), diff(G(xi), xi), diff(diff(G(xi), xi), xi)}

> 

1/subs(Se, 1/rhs(Dubl2))

(1/2)*r^2*G(xi)*(4*l*G(xi)^2+3*b*G(xi)+2*a)

> 

1/dsubs(Se, 1/rhs(Dubl2))

(1/2)*r^2*G(xi)*(4*l*G(xi)^2+3*b*G(xi)+2*a)

> 

NULL

Download subs.mw

Comments...

dharr 8195
April 29 2025
1 5

In Maple, you cannot incorporate two values into one variable by using +/-. You can of course do this manually, e.g. A+pm*B and later set pm to +1 or -1.

To get the exp form, use simplify. (But you will need to set a value of c__1 first to get the form in the paper.) I do not know how to get the sech/tanh form. For these sorts of things I usually enter it manually and then simplify the difference to check it is zero and you haven't entered it incorrectly.

Eigenvectors command...

dharr 8195
April 29 2025
0 0

See attached file (not displaying in Mapleprimes right now).

Download eigenvectors.mw

table of coefficients...

dharr 8195
April 28 2025
1 4

Here I make a table of the coefficients for each monomial and you can pick out which ones you want. The list in your question has repeats, but you can change to what you really want. I'm not sure it makes sense to pick out some coefficients and not others. In the other case the only monomials in the expression were those in the paper.

Download S1.mw

M[sun]...

dharr 8195
April 27 2025
0 4

It works with M[sun] rather than M__sun.

> 

restart;
> 

with(ScientificConstants):    
> 

GetConstant(M__Sun)         

Error, (in ScientificConstants:-GetConstant) `M__Sun` is not a known Constant

> 

GetConstant(M[Sun])  

mass_of_Sun, symbol = M[Sun], value = 0.199e31, uncertainty = 0.5e28, units = kg

> 

NULL

Download Msun.mw

Loop...

dharr 8195
April 27 2025
1 5

You can use a for loop to do this. I didn't wait for it to finish.

Trail-pdetest.mw

Solutions...

dharr 8195
April 26 2025
1 8

Once again you had an extra =0. But your treatment doesn't give the correct Eq. (26) in the paper. I guessed how they might have obtained it, but then there is a sign wrong in that equation. So you could try changing the sign, and double-checking the terms since I modified it by hand (some strange thing with PDETools:-declare got in the way of a cut and paste). I didn't sort through the answers from solve, but I assume the correct solutions are in there somewhere, or will be with the sign change.

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

_local(gamma)
> 

declare(u(x, y, z, t))

u(x, y, z, t)*`will now be displayed as`*u

> 

declare(f(x, y, z, t))

f(x, y, z, t)*`will now be displayed as`*f

> 

pde := 9*(diff(u(x, y, z, t), t, x))+diff(u(x, y, z, t), `$`(x, 6))-5*(diff(u(x, y, z, t), `$`(x, 3), y)+diff(u(x, y, z, t), `$`(y, 2)))+15*((diff(u(x, y, z, t), `$`(x, 2)))*(diff(u(x, y, z, t), `$`(x, 3)))+(diff(u(x, y, z, t), x))*(diff(u(x, y, z, t), `$`(x, 4)))-(diff(u(x, y, z, t), x))*(diff(u(x, y, z, t), x, y))-(diff(u(x, y, z, t), `$`(x, 2)))*(diff(u(x, y, z, t), y)))+45*(diff(u(x, y, z, t), x))^2*(diff(u(x, y, z, t), `$`(x, 2)))+alpha*(diff(u(x, y, z, t), `$`(x, 2)))+beta*(diff(u(x, y, z, t), x, y))+gamma*(diff(u(x, y, z, t), x, z))

9*(diff(diff(u(x, y, z, t), t), x))+diff(diff(diff(diff(diff(diff(u(x, y, z, t), x), x), x), x), x), x)-5*(diff(diff(diff(diff(u(x, y, z, t), x), x), x), y))-5*(diff(diff(u(x, y, z, t), y), y))+15*(diff(diff(u(x, y, z, t), x), x))*(diff(diff(diff(u(x, y, z, t), x), x), x))+15*(diff(u(x, y, z, t), x))*(diff(diff(diff(diff(u(x, y, z, t), x), x), x), x))-15*(diff(u(x, y, z, t), x))*(diff(diff(u(x, y, z, t), x), y))-15*(diff(diff(u(x, y, z, t), x), x))*(diff(u(x, y, z, t), y))+45*(diff(u(x, y, z, t), x))^2*(diff(diff(u(x, y, z, t), x), x))+alpha*(diff(diff(u(x, y, z, t), x), x))+beta*(diff(diff(u(x, y, z, t), x), y))+gamma*(diff(diff(u(x, y, z, t), x), z))

> 

oppde := [op(expand(pde))]; u_occurrences := map(proc (i) options operator, arrow; numelems(select(has, [op([op(i)])], u)) end proc, oppde); linear_op_indices := ListTools:-SearchAll(1, u_occurrences); pde_linear := add(oppde[[linear_op_indices]]); pde_nonlinear := expand(simplify(expand(pde)-pde_linear))

9*(diff(diff(u(x, y, z, t), t), x))+diff(diff(diff(diff(diff(diff(u(x, y, z, t), x), x), x), x), x), x)-5*(diff(diff(diff(diff(u(x, y, z, t), x), x), x), y))-5*(diff(diff(u(x, y, z, t), y), y))+alpha*(diff(diff(u(x, y, z, t), x), x))+beta*(diff(diff(u(x, y, z, t), x), y))+gamma*(diff(diff(u(x, y, z, t), x), z))

15*(diff(diff(u(x, y, z, t), x), x))*(diff(diff(diff(u(x, y, z, t), x), x), x))+15*(diff(u(x, y, z, t), x))*(diff(diff(diff(diff(u(x, y, z, t), x), x), x), x))-15*(diff(u(x, y, z, t), x))*(diff(diff(u(x, y, z, t), x), y))-15*(diff(diff(u(x, y, z, t), x), x))*(diff(u(x, y, z, t), y))+45*(diff(u(x, y, z, t), x))^2*(diff(diff(u(x, y, z, t), x), x))

> 

H := u(x, y, z, t) = 2*(diff(ln(f(x, y, z, t)), x))

u(x, y, z, t) = 2*(diff(f(x, y, z, t), x))/f(x, y, z, t)

> 

H1 := int(pde_linear, x)

(diff(u(x, y, z, t), z))*gamma+alpha*(diff(u(x, y, z, t), x))+beta*(diff(u(x, y, z, t), y))+9*(diff(u(x, y, z, t), t))-5*(int(diff(diff(u(x, y, z, t), y), y), x))+diff(diff(diff(diff(diff(u(x, y, z, t), x), x), x), x), x)-5*(diff(diff(diff(u(x, y, z, t), x), x), y))

No =0 here

> 

L := eval(H1, H)
> 

F1 := (1/2)*numer(normal(L))

The following does not agree with Eq. (26) in the paper

> 

collect(F1, {alpha, beta, gamma, f(x, y, z, t)})

(f(x, y, z, t)^5*(diff(diff(f(x, y, z, t), x), x))-f(x, y, z, t)^4*(diff(f(x, y, z, t), x))^2)*alpha+(f(x, y, z, t)^5*(diff(diff(f(x, y, z, t), x), y))-f(x, y, z, t)^4*(diff(f(x, y, z, t), y))*(diff(f(x, y, z, t), x)))*beta+((diff(diff(f(x, y, z, t), x), z))*f(x, y, z, t)^5-(diff(f(x, y, z, t), z))*f(x, y, z, t)^4*(diff(f(x, y, z, t), x)))*gamma+(9*(diff(diff(f(x, y, z, t), t), x))+diff(diff(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x), x), x)-5*(diff(diff(diff(diff(f(x, y, z, t), x), x), x), y))-5*(diff(diff(f(x, y, z, t), y), y)))*f(x, y, z, t)^5+(-6*(diff(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x), x))*(diff(f(x, y, z, t), x))-15*(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x))*(diff(diff(f(x, y, z, t), x), x))-10*(diff(diff(diff(f(x, y, z, t), x), x), x))^2+5*(diff(diff(diff(f(x, y, z, t), x), x), x))*(diff(f(x, y, z, t), y))-9*(diff(f(x, y, z, t), x))*(diff(f(x, y, z, t), t))+15*(diff(diff(f(x, y, z, t), x), x))*(diff(diff(f(x, y, z, t), x), y))+5*(diff(f(x, y, z, t), y))^2+15*(diff(diff(diff(f(x, y, z, t), x), x), y))*(diff(f(x, y, z, t), x)))*f(x, y, z, t)^4+(30*(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x))*(diff(f(x, y, z, t), x))^2+120*(diff(diff(diff(f(x, y, z, t), x), x), x))*(diff(f(x, y, z, t), x))*(diff(diff(f(x, y, z, t), x), x))-30*(diff(f(x, y, z, t), x))^2*(diff(diff(f(x, y, z, t), x), y))-30*(diff(diff(f(x, y, z, t), x), x))*(diff(f(x, y, z, t), x))*(diff(f(x, y, z, t), y))+30*(diff(diff(f(x, y, z, t), x), x))^3)*f(x, y, z, t)^3+(-120*(diff(diff(diff(f(x, y, z, t), x), x), x))*(diff(f(x, y, z, t), x))^3+30*(diff(f(x, y, z, t), x))^3*(diff(f(x, y, z, t), y))-270*(diff(diff(f(x, y, z, t), x), x))^2*(diff(f(x, y, z, t), x))^2)*f(x, y, z, t)^2+360*(diff(diff(f(x, y, z, t), x), x))*(diff(f(x, y, z, t), x))^4*f(x, y, z, t)-120*(diff(f(x, y, z, t), x))^6

Take only terms with f^4 or higher and divide them by f^4. Agrees with Eq.  (26) except has +15*fxx*fyy instead of -15*fxx*fyy

> 

F1:=(f(x,y,z,t)*diff(diff(f(x,y,z,t),x),x)-diff(f(x,y,z,t),x)^2)*alpha
+(f(x,y,z,t)*diff(diff(f(x,y,z,t),x),y)-diff(f(x,y,z,t),y)*diff(f(x,y,z,t),x))*beta
+(-diff(f(x,y,z,t),x)*diff(f(x,y,z,t),z)+f(x,y,z,t)*diff(diff(f(x,y,z,t),x),z))*gamma
+(diff(diff(diff(diff(diff(diff(f(x,y,z,t),x),x),x),x),x),x)-5*diff(diff(diff(diff(f(x,y,z,t),x),x),x),y)+9*diff(
diff(f(x,y,z,t),t),x)-5*diff(diff(f(x,y,z,t),y),y))*f(x,y,z,t)
+(-6*diff(diff(
diff(diff(diff(f(x,y,z,t),x),x),x),x),x)*diff(f(x,y,z,t),x)-9*diff(f(x,y,z,t),x
)*diff(f(x,y,z,t),t)+15*diff(diff(diff(f(x,y,z,t),x),x),y)*diff(f(x,y,z,t),x)
+15*diff(diff(f(x,y,z,t),x),x)*diff(diff(f(x,y,z,t),x),y)
-15*diff(diff(diff(diff(f(x,y,z,t),x),x),x),x)*diff(diff(f(x,y,z,t),x),x)+5*diff(f(x,y,z,t),y)^2+5*
diff(diff(diff(f(x,y,z,t),x),x),x)*diff(f(x,y,z,t),y)-10*diff(diff(diff(f(x,y,z
,t),x),x),x)^2);

(f(x, y, z, t)*(diff(diff(f(x, y, z, t), x), x))-(diff(f(x, y, z, t), x))^2)*alpha+(f(x, y, z, t)*(diff(diff(f(x, y, z, t), x), y))-(diff(f(x, y, z, t), x))*(diff(f(x, y, z, t), y)))*beta+(-(diff(f(x, y, z, t), x))*(diff(f(x, y, z, t), z))+f(x, y, z, t)*(diff(diff(f(x, y, z, t), x), z)))*gamma+(9*(diff(diff(f(x, y, z, t), t), x))+diff(diff(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x), x), x)-5*(diff(diff(diff(diff(f(x, y, z, t), x), x), x), y))-5*(diff(diff(f(x, y, z, t), y), y)))*f(x, y, z, t)-6*(diff(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x), x))*(diff(f(x, y, z, t), x))-15*(diff(diff(diff(diff(f(x, y, z, t), x), x), x), x))*(diff(diff(f(x, y, z, t), x), x))-10*(diff(diff(diff(f(x, y, z, t), x), x), x))^2+5*(diff(diff(diff(f(x, y, z, t), x), x), x))*(diff(f(x, y, z, t), y))-9*(diff(f(x, y, z, t), x))*(diff(f(x, y, z, t), t))+15*(diff(diff(f(x, y, z, t), x), x))*(diff(diff(f(x, y, z, t), x), y))+5*(diff(f(x, y, z, t), y))^2+15*(diff(diff(diff(f(x, y, z, t), x), x), y))*(diff(f(x, y, z, t), x))

> 

NULL
> 

T := f(x, y, z, t) = g(x, y, z, t)^2+h(x, y, z, t)^2+a[11]
> 

T1 := g(x, y, z, t) = t*a[4]+x*a[1]+y*a[2]+z*a[3]+a[5]
> 

T2 := h(x, y, z, t) = t*a[9]+x*a[6]+y*a[7]+z*a[8]+a[10]
> 

L2 := subs({T1, T2}, T)

f(x, y, z, t) = (t*a[4]+x*a[1]+y*a[2]+z*a[3]+a[5])^2+(t*a[9]+x*a[6]+y*a[7]+z*a[8]+a[10])^2+a[11]

> 

L3 := eval(F1, L2)
> 

L4 := collect(L3, [x, y, z, t], 'distributed')
> 

L5 := {coeffs(L4, [x, y, z, t], 'monomials')}; nops(%)

15

> 

monomials

1, x, z, z*t, y*t, y*z, x*z, x*t, x*y, y, t, x^2, y^2, z^2, t^2

> 

indets(L5); nops(%)

{alpha, beta, gamma, a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11]}

14

> 

ans := CodeTools:-Usage([solve(L5), explicit]); nops(%)

memory used=2.66GiB, alloc change=102.89MiB, cpu time=16.78s, real time=15.00s, gc time=3.03s

38

> 

ans[1]

{alpha = alpha, beta = beta, gamma = gamma, a[1] = 0, a[2] = RootOf(_Z^2+1)*a[7], a[3] = 0, a[4] = a[9]*(RootOf(_Z^2+1)-1)/(RootOf(_Z^2+1)+1), a[5] = a[5], a[6] = 0, a[7] = a[7], a[8] = 0, a[9] = a[9], a[10] = -a[5]*RootOf(_Z^2+1), a[11] = a[11]}

 

NULL

Download systems.mw

polylog...

dharr 8195
April 24 2025
2 1

Maple 2025 also gives the polylog form, which appears to be real for 0<a<b<Pi.
@Zeineb So although "I" appears, the expression is actually real if you put numbers in.

@Alfred_F So this is a closed form IMO, though not in terms of elementary functions. I tried separating this into real and imaginary parts with evalc but no luck so I think there isn't a simpler form (other than perhaps combining some terms), in agreement with your comment. Nor does making assumptions about a,b simplify it any further.

> 

restart;
> 

integrand:=y^2*cot(y);

y^2*cot(y)

> 

plot(integrand,y=0..2.5);

Running integral from 0

> 

Qb:=int(integrand,y=0..b);

-(1/2)*Zeta(3)-((1/3)*I)*b^3+b^2*ln(1+exp(I*b))-(2*I)*b*polylog(2, -exp(I*b))+2*polylog(3, -exp(I*b))+b^2*ln(1-exp(I*b))-(2*I)*b*polylog(2, exp(I*b))+2*polylog(3, exp(I*b))

Looks correct

> 

plot(Qb, b=0..2.5);

The integral we want

> 

Q:=simplify(int(integrand, y=a..b));

((1/3)*I)*a^3-a^2*ln(1+exp(I*a))+I*a*polylog(2, exp((2*I)*a))-(1/2)*polylog(3, exp((2*I)*a))-a^2*ln(1-exp(I*a))-((1/3)*I)*b^3+b^2*ln(1+exp(I*b))-I*b*polylog(2, exp((2*I)*b))+(1/2)*polylog(3, exp((2*I)*b))+b^2*ln(1-exp(I*b))

> 

eval(Q,{a=1.1,b=2.5});

-2.353832100+0.1e-8*I

Complex above Pi

> 

eval(Q,{a=1.1,b=3.5});

-1.691504809-31.00627668*I

> 

plot3d(Q,b=0..Pi, a=0..b, view=[default,default,-10..2]);

> 

NULL

Download polylog.mw

tiled hexagons...

dharr 8195
April 16 2025
2 0

Not optimum. Adapted from some old code I had.

> 

restart;
> 

chexagons:=proc(n::posint,a::positive)
   # n = lattice size, a = hexagon size
   local i, j, imax, hexesA, hexesB, hexesC, hex;
   uses plottools;
   # hex returns coords of hexagon of edge length a centred at x,y.
   hex := (x,y,a) -> [[x+a,y],[x+a/2,y+a*sqrt(3)/2],[x-a/2,y+a*sqrt(3)/2],[x-a,y],[x-a/2,y-a*sqrt(3)/2],[x+a/2,y-a*sqrt(3)/2]];
   hexesA := NULL; hexesB := NULL; hexesC := NULL;
   for i to n do for j to n do
     if j-i mod 3 =0 then
       if i mod 3 = 0 then hexesA := hexesA, hex((j-1)*a+(i-1)*a/2,-(i-1)*a*sqrt(3)/2,a) end if;
       if i mod 3 = 1 then hexesB := hexesB, hex((j-1)*a+(i-1)*a/2,-(i-1)*a*sqrt(3)/2,a) end if;
       if i mod 3 = 2 then hexesC := hexesC, hex((j-1)*a+(i-1)*a/2,-(i-1)*a*sqrt(3)/2,a) end if;
     end if;
   end do end do;
   plots:-display(polygon([hexesA],color=red), polygon([hexesB],color=blue), polygon([hexesC],color=green),
           scaling=constrained,axes=none)
end proc:
> 

chexagons(10,1);

> 

 

Download hexagons.mw

naming state variable...

dharr 8195
April 14 2025
3 4

This is an answer for (1). For (2), would you please upload your worksheet (using green up-arrow) so we can understand more clearly what you want.

p.s. q_d seems to be ignored here; not sure what that means or what you wanted.

> 

restart
> 

with(DynamicSystems)

The default discretetimevar is q, so to use n

> 

SystemOptions(discretetimevar = n)

q

> 

sys2 := StateSpace([u(n) = K_virt*(q_a(n)-q_d)+D_virt*(q_a(n)-q_a(n-1))/T], inputvariable = [q_a(n)], outputvariable = [u(n)], discrete = true, sampletime = T)

Look at the matrices. Since we have c <> 0, the state variable is being used and is implicitly the delayed value - note that c is the coefficient of q_a(n-1)

> 

PrintSystem(sys2)

"[[[`State Space`],[`discrete; sampletime = T`],[`1 output(s); 1 input(s); 1 state(s)`],[inputvariable=[q_a(n)]],[outputvariable=[u(n)]],[statevariable=[x1(n)]],[a=[[[0]]]],[b=[[[1]]]],[c=[[[-D_virt/T]]]],[d=[[[(K_virt T+D_virt)/T]]]]]"

If you want to call it q_a_delayed, you have to explicitly invoke it in the difference equations

> 

sys3 := StateSpace([u(n) = K_virt*(q_a(n)-q_d)+D_virt*(q_a(n)-q_a_delayed(n))/T, q_a_delayed(n) = q_a(n-1)], inputvariable = [q_a(n)], outputvariable = [u(n)], statevariable = [q_a_delayed(n)], discrete = true, sampletime = T)

Same matrices

> 

PrintSystem(sys3)

"[[[`State Space`],[`discrete; sampletime = T`],[`1 output(s); 1 input(s); 1 state(s)`],[inputvariable=[q_a(n)]],[outputvariable=[u(n)]],[statevariable=[q_a_delayed(n)]],[a=[[[0]]]],[b=[[[1]]]],[c=[[[-D_virt/T]]]],[d=[[[(K_virt T+D_virt)/T]]]]]"

> 

NULL

Download DynamicSystems.mw

missing =...

dharr 8195
April 14 2025
2 1

In your params definition, you have varphi__8.33, which should instead be something like varphi__X=8.33.

braces...

dharr 8195
April 12 2025
0 0

Solve has an alternative way of working with braces, with slightly different syntax. 

braces.mw

2 3 4 5 6 7 8 Last Page 4 of 81