dharr

Dr. David Harrington

9207 Reputation

22 Badges

21 years, 300 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 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

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

@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.

@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.

@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.

Can't be sure it is a global maximum. [Edit: updated to use Euler angles, which gives a larger maximum than my first attempt.]

restart

with(geom3d); _local(O, gamma)

Let C and O be the cube and octahedron. The cube has sides of length 1. Take the circumsphere of the octahedron to be 1 for now.

point(o, 0, 0, 0); cube(C, o, side = 1); octahedron(O, o, 1)

Rotate the octahedron by Euler angles alpha, beta, gamma.

line(zax, [0, 0, t], t); line(xax, [t, 0, 0], t); rotation(OR, rotation(Ozx, rotation(Oz, O, gamma, zax), beta, xax), alpha, zax)

Scale by r and consider the vertices

verts := map(proc (q) options operator, arrow; `~`[`*`](r, q) end proc, vertices(OR))

[[r*(cos(alpha)*cos(gamma)-sin(alpha)*cos(beta)*sin(gamma)), r*(cos(alpha)*cos(beta)*sin(gamma)+sin(alpha)*cos(gamma)), r*sin(beta)*sin(gamma)], [r*(-cos(alpha)*cos(gamma)+sin(alpha)*cos(beta)*sin(gamma)), r*(-cos(alpha)*cos(beta)*sin(gamma)-sin(alpha)*cos(gamma)), -r*sin(beta)*sin(gamma)], [r*(-cos(alpha)*sin(gamma)-sin(alpha)*cos(beta)*cos(gamma)), r*(cos(alpha)*cos(beta)*cos(gamma)-sin(alpha)*sin(gamma)), r*sin(beta)*cos(gamma)], [r*(cos(alpha)*sin(gamma)+sin(alpha)*cos(beta)*cos(gamma)), r*(-cos(alpha)*cos(beta)*cos(gamma)+sin(alpha)*sin(gamma)), -r*sin(beta)*cos(gamma)], [r*sin(alpha)*sin(beta), -r*cos(alpha)*sin(beta), r*cos(beta)], [-r*sin(alpha)*sin(beta), r*cos(alpha)*sin(beta), -r*cos(beta)]]

The constraints are that each coordinate lies between -1/2 and +1/2, so no vertex "sticks out" of the cube

c1 := `~`[`>=`](`~`[op](verts), -1/2); c2 := `~`[`<=`](`~`[op](verts), 1/2); constr := remove(is, {c1[], c2[]}, true)

{-1/2 <= r*(-cos(alpha)*cos(gamma)+sin(alpha)*cos(beta)*sin(gamma)), -1/2 <= r*(cos(alpha)*cos(gamma)-sin(alpha)*cos(beta)*sin(gamma)), -1/2 <= r*(-cos(alpha)*sin(gamma)-sin(alpha)*cos(beta)*cos(gamma)), -1/2 <= r*(cos(alpha)*sin(gamma)+sin(alpha)*cos(beta)*cos(gamma)), -1/2 <= r*(-cos(alpha)*cos(beta)*cos(gamma)+sin(alpha)*sin(gamma)), -1/2 <= r*(cos(alpha)*cos(beta)*cos(gamma)-sin(alpha)*sin(gamma)), -1/2 <= r*(-cos(alpha)*cos(beta)*sin(gamma)-sin(alpha)*cos(gamma)), -1/2 <= r*(cos(alpha)*cos(beta)*sin(gamma)+sin(alpha)*cos(gamma)), -1/2 <= r*cos(beta), -1/2 <= r*cos(alpha)*sin(beta), -1/2 <= r*sin(alpha)*sin(beta), -1/2 <= r*sin(beta)*cos(gamma), -1/2 <= r*sin(beta)*sin(gamma), -1/2 <= -r*cos(beta), -1/2 <= -r*cos(alpha)*sin(beta), -1/2 <= -r*sin(alpha)*sin(beta), -1/2 <= -r*sin(beta)*cos(gamma), -1/2 <= -r*sin(beta)*sin(gamma), r*(-cos(alpha)*cos(gamma)+sin(alpha)*cos(beta)*sin(gamma)) <= 1/2, r*(cos(alpha)*cos(gamma)-sin(alpha)*cos(beta)*sin(gamma)) <= 1/2, r*(-cos(alpha)*sin(gamma)-sin(alpha)*cos(beta)*cos(gamma)) <= 1/2, r*(cos(alpha)*sin(gamma)+sin(alpha)*cos(beta)*cos(gamma)) <= 1/2, r*(-cos(alpha)*cos(beta)*cos(gamma)+sin(alpha)*sin(gamma)) <= 1/2, r*(cos(alpha)*cos(beta)*cos(gamma)-sin(alpha)*sin(gamma)) <= 1/2, r*(-cos(alpha)*cos(beta)*sin(gamma)-sin(alpha)*cos(gamma)) <= 1/2, r*(cos(alpha)*cos(beta)*sin(gamma)+sin(alpha)*cos(gamma)) <= 1/2, r*cos(beta) <= 1/2, r*cos(alpha)*sin(beta) <= 1/2, r*sin(alpha)*sin(beta) <= 1/2, r*sin(beta)*cos(gamma) <= 1/2, r*sin(beta)*sin(gamma) <= 1/2, -r*cos(beta) <= 1/2, -r*cos(alpha)*sin(beta) <= 1/2, -r*sin(alpha)*sin(beta) <= 1/2, -r*sin(beta)*cos(gamma) <= 1/2, -r*sin(beta)*sin(gamma) <= 1/2}

So we want to maximize r subject to these constraints

ans := Optimization:-Maximize(r, constr, alpha = -Pi .. Pi, gamma = -Pi .. Pi, beta = 0 .. Pi, initialpoint = {alpha = (1/4)*Pi, beta = (1/4)*Pi, gamma = (1/4)*Pi, r = 1})

[.749999999968849584, [alpha = HFloat(0.7853981633974485), beta = HFloat(1.230959417307807), gamma = HFloat(0.7853981633974483), r = HFloat(0.7499999999688496)]]

Evidently r is 3/4, alpha=gamma = Pi/4

identify(eval(alpha, ans[2])); identify(eval(beta, ans[2])); identify(eval(gamma, ans[2]))

(1/4)*Pi

HFloat(1.230959417307807)

(1/4)*Pi

Some further exploration at higher Digits suggests beta = arccos(1/3)

arccos(1/3); evalf(%)

arccos(1/3)

1.230959417

So go back and scale and rotate a new octahedron accordingly

octahedron(O2, o, 3/4); rotation(OR2, rotation(Ozx2, rotation(Oz2, O2, (1/4)*Pi, zax), arccos(1/3), xax), (1/4)*Pi, zax)

draw([C(cutout = .85), OR2])

And the maximum volume is

volume(OR2)

9/16

``

Download cubeoctOptimization2.mw

Even though kernelopts(version) is incorrect, interface(version) correctly reports 2026.1.

@C_R Converting the lprint(%) to 1-D resolves the issue. But I was not able to generate it with a fresh 2026.1 worksheet (I didn't try document mode). But I have a recollection that recently 2D math with % in has not worked consistently. So it might not be to do wth lprint but with % in 2D. Unfortunately I can't get it to not work when I try a few examples. If I come up with one I'll post it here.

You added a second question with the uploaded file instead of editing your own question and uploading the file to it. So I did this for you, moving your file to where it should have been. I then deleted the now redundant question, with the reason "duplicate question" and not "spam", though perhaps the system shows you that.

I am happy to accept donations for my brain upgrade. On the other hand you may wish to delete this post.

Maple's standard diff can differentiate with respect to theta, but not to theta(s). However diff is modified to accept theta(s) when the Physics package is loaded. So with(Physics) solves the error and returns some output, but I don't know if it is what you wanted.

@Mth2026 That download page for the assistant implies that it is not included with Maple Flow but has to be downloaded separately.

There is a link to download it at https://www.maplesoft.com/products/mapleflow/migration-assistant/ 

If you are having diffculties with the download and installation, you should contact technical support.

@Ronan Well one documented way to stay in the superscipt mode after x^2 is to enter a space; then + stays on the same line as 2. The same works here, i.e.,  1/2 ! works as expected. But one doesn't want to have to think about this non-intuitive solution.

A version or two ago, 2D math was modified so that the keystrokes x^2+ automatically brought the + sign back to the baseline. At that time, the behaviour you observed was introduced. I assume this was an unintended side effect of the changes, and I consider it a bug. There was some discussion of this before and I am assuming it is on a list of things to be fixed; but if not please add it to the list, Maplesoft.

My recollection is that there are other cases where the cursor jumps out of the denominator, but now that I try to find them I'm out of luck, so maybe I'm wrong about this.

@Harry Garst So if I have some Matrices, I might want to manipulate them with their entries. In your Zehfuss example, I can see why you might want to use that rule to evaluate Determinant(Kronecker(A,B) since using that rule gives a much simpler expression than the direct evaluation. But for that purpose, you could write a procedure DetKron:=(A,B)->Determinant(A)^etc.

I can see a different application where you want to simplify some expression with abstract matrices A,B that are just the symbols A,B but have the properties of matrices (in particular non-cummutativity) - such as implementing these Knonecker properties.

You seem to want to make your life difficult by wanting to have A,B real matrices and then not use that fact - requiring uneval and uneval quotes. So everything has to be done inside a call to ZehfussReplace or similar. Or you can use inert operators and then use value on your result. But it is awkward and not the way I tend to use Maple.

For the abstract case, you can do some things with define or pattern matching. But I think the physics package is a better option. QuantumOperators act like matrices, and I suspect that Knonecker products are also there somewhere since in an abstract sense they are just tensor products. So these rules may already exist. I don't use the package much so I can't be more specific.

In general (for symbols or real matrices), if I want to know if two things are equal, I simplify their difference and see if it is zero (or the zero matrix). And in general I manipulate things in abstract form and when I get to the end of those manipulations I put in the numbers with eval. That is a bit harder with matrices, but try

A . B;
eval(%, {A = Matrix(2, 2, symbol = a), B = Matrix(2, 2, symbol = b)});

or for an expanded polynomial

2*x^3 + x^2 + 2;
eval(%, x = Matrix(2, 2, [1, 0, 3, 5]));

@acer That's much better. I struggled with getting the two square matrices into the same type specification. For the inerts, I was thinking that the OP might be wanting an inert output, e.g.,

ZehfussReplace := proc(expr::uneval)
  uses LinearAlgebra;
  subsindets(expr,
    specfunc(specfunc('KroneckerProduct'),'Determinant'),
    proc(q)
      local AB;
      AB:=op(op(q));
      if not eval([AB])::['Matrix'('square'),'Matrix'('square')] then return eval(q) end if;
      %Determinant(AB[1])^(upperbound(AB[2])[1])*%Determinant(AB[2])^(upperbound(AB[1])[1])
    end proc
   );
end proc:

(or using @acers version),and then value could be used on the output.

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