Advanced Search

dharr

Dr. David Harrington

9152 Reputation

22 Badges

21 years, 261 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 replies submitted by dharr

Back to Fig 1 - corrected...

dharr 9152
9 hours ago
0 0

Here's further back to psi and fig 1. I had v=u^3/3 instead of v=u^2/3 earlier; now everything looks correct.

Subcases (1) and (3)

> 

restart
> 

with(PDEtools)
> 

undeclare(prime, quiet)
> 

declare(u(xi), quiet); declare(w(`ξ__1`), quiet)

Eq 3.5

> 

ode := (diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+g[0]

(diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+g[0]

Eq 3.6 with epsilon = 1

> 

itr := {`ξ__1` = delta^(1/4)*xi, w(`ξ__1`) = delta^(1/4)*u(xi)}

{xi__1 = delta^(1/4)*xi, w(xi__1) = delta^(1/4)*u(xi)}

> 

tr := solve(itr, {xi, u(xi)})

{xi = xi__1/delta^(1/4), u(xi) = w(xi__1)/delta^(1/4)}

Eq. 3.7

> 

ode2 := dchange(tr, ode, [`ξ__1`, w(`ξ__1`)], params = {delta})

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+(beta*p^3+p*q-r)*w(xi__1)^2/(delta^(1/2)*(3*beta*p+1))+g[0]

Eq 3.7 in terms of e2 and e0

> 

ode3 := subs(1/sqrt(delta) = e__2*(3*beta*p+1)/(beta*p^3+p*q-r), g[0] = e__0, ode2)

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+e__2*w(xi__1)^2+e__0

> 

F := subs(w(`ξ__1`) = W, rhs(ode3))

W^4+W^2*e__2+e__0

> 

F2 := collect(((W-w__1)^2+w__2^2)^2, W)

W^4-4*w__1*W^3+(6*w__1^2+2*w__2^2)*W^2-4*(w__1^2+w__2^2)*w__1*W+(w__1^2+w__2^2)^2

 

For the coefficients of w^3 to be the same we need w1 = 0. Then for the coeffs of w^2 and 1 we need to solve

> 

F21 := eval(F2, w__1 = 0); eqs := {coeff(F21, W^2) = coeff(F, W^2), coeff(F21, W, 0) = coeff(F, W, 0)}

W^4+2*W^2*w__2^2+w__2^4

{w__2^4 = e__0, 2*w__2^2 = e__2}

There is no solution for w2 unless we have some special relationship between e[2] and e[0]

> 

eliminate(eqs, w__2); eq := %[1][2][]; eq0 := e__0 = solve(eq, e__0)

[{w__2 = -(1/2)*2^(1/2)*e__2^(1/2)}, {-e__2^2+4*e__0}], [{w__2 = (1/2)*2^(1/2)*e__2^(1/2)}, {-e__2^2+4*e__0}]

-e__2^2+4*e__0

e__0 = (1/4)*e__2^2

So with this substitution we can work out the integral. The result has a slightly different form since w1 =0

> 

factor(eval(F, eq0)); `assuming`([Int(1/sqrt(%), W)], [2*W^2+e__2 > 0]); eqxi1 := `ξ__1` = value(%)+`ξ__10`

(1/4)*(2*W^2+e__2)^2

Int(1/(W^2+(1/2)*e__2), W)

xi__1 = 2^(1/2)*arctan(W*2^(1/2)/e__2^(1/2))/e__2^(1/2)+xi__10

So the solution to ode3 for this is

> 

sol3 := w(`ξ__1`) = solve(eqxi1, W)

w(xi__1) = (1/2)*tan((1/2)*e__2^(1/2)*(xi__1-xi__10)*2^(1/2))*e__2^(1/2)*2^(1/2)

> 

odetest(sol3, eval(ode3, eq0))

0

Convert back to earlier constants. We have e0 = g[0] =e2^2/4 so change

> 

eqe2 := e__2 = (beta*p^3+p*q-r)/(sqrt(delta)*(3*beta*p+1)); eqg0 := g[0] = rhs(eval(eq0, eqe2)); ode2b := eval(ode2, eqg0); sol2 := simplify(eval(sol3, eqe2))

e__2 = (beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1))

g[0] = (1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+(beta*p^3+p*q-r)*w(xi__1)^2/(delta^(1/2)*(3*beta*p+1))+(1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

w(xi__1) = (1/2)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(xi__1-xi__10)*2^(1/2))*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*2^(1/2)

> 

odetest(sol2, ode2b)

0

So now transform back to u(xi)

> 

itr2 := {`ξ__10` = delta^(1/4)*`ξ__0`}
> 

solu := dchange(`union`(itr, itr2), sol2/delta^(1/4), [u(xi), xi, `ξ__0`], params = {delta}); odeb := eval(ode, eqg0)

u(xi) = (1/2)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*xi-delta^(1/4)*xi__0)*2^(1/2))*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*2^(1/2)/delta^(1/4)

(diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+(1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

> 

odetest(solu, odeb)

0

So v(xi) is, from Eq. 3.2

> 

solv := v(xi) = eval(-(1/3)*u(xi)^2, solu)

v(xi) = -(1/6)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*xi-delta^(1/4)*xi__0)*2^(1/2))^2*(beta*p^3+p*q-r)/(delta*(3*beta*p+1))

And so psi is perhaps? (There is something about taking the real part I didn't understand.)

> 

psi := MakeFunction(eval(rhs(solv), xi = c*t+x+y), x, y, t)

proc (x, y, t) options operator, arrow; -(1/6)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*(c*t+x+y)-delta^(1/4)*xi__0)*2^(1/2))^2*(beta*p^3+p*q-r)/(delta*(3*beta*p+1)) end proc

Figure 1 (psi may be scaled differently). Not sure which is x and which is t

> 

psi1 := eval(psi(x, 0, t), {`ξ__0` = 0, beta = 1, c = 1/4, delta = 1/6, p = -1/2, q = 2*sqrt(6)*(1/3), r = -1/8})

-(2/3)*tan((1/12)*4^(1/2)*6^(3/4)*((1/4)*t+x)*2^(1/2))^2*6^(1/2)

> 

plot3d(psi1, t = -10 .. 10, x = -10 .. 10, view = [default, default, -40 .. 2], orientation = [55, 55])

> 

plot(eval(psi1, t = 0), x = -10 .. 10, -40 .. 2)

Subcase (3) is the same function, just written in a different way. Fig. 3:

> 

psi4 := eval(psi(x, 0, t), {`ξ__0` = 0, beta = 1, c = 1/4, delta = 1/6, p = -1/2, q = -(1/2)*sqrt(6), r = -1/8})

(1/2)*tan((1/12)*(-3)^(1/2)*6^(3/4)*((1/4)*t+x)*2^(1/2))^2*6^(1/2)

> 

plot3d(psi4, t = -5 .. 5, x = -5 .. 5, view = [default, default, -2 .. 0], orientation = [55, 55])

NULL

Download fp3.mw

transform back to Eq 3.5...

dharr 9152
10 hours ago
0 0

@salim-barzani Here's how to transform back to Eq 3.5, which was the original starting point of this question

fp3.mw

subcases...

dharr 9152
14 hours ago
0 0

I think the subcases in the paper are not all different when e1=0, so it is simpler to go through them without the classification scheme. Knowing when the quartic is a perfect square is useful to simplify the integral, but it is not clear to me that knowing when the various roots are real or not has much to do with evaluating the integral. Probably the solutions can be interconverted. In the end, Maple can just solve the ode directly, so I'm not sure I can see the point.

roots.mw

interesting...

dharr 9152
17 hours ago
0 0

@Alfred_F That's a nice observation. e0 is lost, but there are now two integration constants. But the square root is gone, which I think is more important.

simpler...

dharr 9152
Yesterday at 11:50 PM
0 0

Here's a simpler way to find two of the solutions, and of course the general solution is easy to find.

roots2.mw

special case...

dharr 9152
Yesterday at 7:32 PM
0 0

@dharr I found a slightly different tan solution which only works for a special relationship between e0 and e2, so probably that is what is meant. 

fp3.mw

no tan solution...

dharr 9152
Yesterday at 3:16 PM
0 0

@salim-barzani OK, lets take epsilon=1 then I agree with Eq 3.7. I needed W because you can't integrate wrt a function in Maple, i.e. int(..., w(xi__1)) doesn't work.

But as I showed there is no way to transform F to F2, so the integral that gives the arctan result is invalid, and so Eq 3.13 is incorrect, so it will not satisfy the ode. You can perfom the integral with e2 and e0 in (just change the Int to int). The problem is then that you get xi__1 = f(w) but here and in general f(w) has multiple w in, so you cannot easily invert it to get w as a function of xi__1.

The solution to this is to solve the ode in Eq 3.7. dsolve will do this directly, or you could take the square root and integrate. Since you get a general solution it will work for arbitrary e2 and e0. If you want to then find cases that are real for various choices of e2 and e0 you can make assumptions that make the things under the square roots positive.

closer...

dharr 9152
May 11 2026
0 0

@salim-barzani Almost there.

try...

dharr 9152
May 11 2026
0 0

@salim-barzani Well if it were changed I might look at it. I have several times pointed out that it makes it harder to help you, and I have ignored several of your recent questions for this reason. In the present case using w instead of r will be especially difficult since the paper uses w later for a different purpose.

notation...

dharr 9152
May 11 2026
0 0

I notice that as usual you changed the notation in your worksheet to make it harder to follow the paper. Perhaps you wish to correct that.

Works...

dharr 9152
May 10 2026
0 0

@sand15 Sorry; I must have google searched it instead of clicking the link

link problem...

dharr 9152
May 10 2026
0 0

@sand15 Your second link does not work (not a web address).

nice numbers...

dharr 9152
May 10 2026
0 0

@Alfred_F It is a pity that the nice numbers aren't really proven because the maximization works in floating point numbers, so r = 0.749999999968849584 we guess is 3/4, and identify guesses at the exact values of the angles. Possibly by doing the differentiations involved in the gradients involved in the maximization and solving, you could do it all symbolically.

geom3d is a nice package and easily works here for the duals since the radius of the circumsphere and in-sphere are built-in. It has a duality command that I didn't make use of. One limitation, unlike its 2D equivalent geometry, is that the sizes have to be fixed at the time of creation, so you can't have a cube with sides a, where a is a name.

improved...

dharr 9152
May 09 2026
0 0

@dharr Yes it works with 3 Euler angles; see updated answer. I should have realized that right away. As usual you have to play around with the initial point.

agreed...

dharr 9152
May 09 2026
0 0

@Alfred_F If I do the further tilt you describe, I can get 9/16.

cubeoctOptimization2.mw

I suspect I need 3 Euler rotations, so I'll try that and update if it works.

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