In the excellent book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the theorem is proved which states that any bounded planar region can be cut into 4 regions of equal area by 2 perpendicular cuts (the pancake problem). This is an existence theorem which does not provide any way to find these cuts. In this post I made an attempt to find such cuts for any convex region on the plane bounded by a piecewise smooth self-non-intersecting curve.
The Into_4_Equal_Areas procedure returns a list of coordinates of 5 points: the first 4 points are the endpoints of the cutting segments, the fifth point is the intersection point of these segments. This procedure significantly uses my old procedure Area , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal argument of the Into_4_Equal_Areas procedure is a list  L specifying the boundary of the region to be cut. When specifying  L, the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source region and cutting segments. For ease of use, the code for the  Area  procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Into_4_Equal_Areas:=proc(L::list,N::symbol:='OneSolution', eps::numeric:=0)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, S, sol, Sol;
f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
S:=Area(L); c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
if not ((nops({i,j,k})=1 and type(L[i],listlist)) or (nops({j,k,m})=1 and type(L[j],listlist)))then
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
P:=eval([x,y], solve({f(A,C),f(B,D)},{x,y})):
L1:=`if`(j=i,[subsop([2,2]=t1..t2,L0[i]),[convert(B,list),P],[P,convert(A,list)]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]], [subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),L0[i+1..j-1][],subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]])):
L2:=`if`(k=j,[subsop([2,2]=t2..t3,L0[j]),[convert(C,list),P],[P,convert(B,list)]],`if`(k=j+1,[subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]], [subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),L0[j+1..k-1][],subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]])):
L3:=`if`(m=k,[subsop([2,2]=t3..t4,L0[k]),[convert(D,list),P],[P,convert(C,list)]],`if`(m=k+1,[subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]], [subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),L0[k+1..m-1][],subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]])):
sol:=fsolve({Area(L1)-S/4,Area(L2)-S/4,Area(L3)-S/4,LinearAlgebra:-DotProduct(D-B,C-A, conjugate=false)},{t1=op([2,2,1],L0[i])-eps..op([2,2,2],L0[i])+eps,t2=op([2,2,1],L0[j])-eps..op([2,2,2],L0[j])+eps,t3=op([2,2,1],L0[k])-eps..op([2,2,2],L0[k])+eps,t4=op([2,2,1],L0[m])-eps..op([2,2,2],L0[m])+eps}) assuming real:
if type(sol,set(`=`)) then if N='OneSolution' then return convert~(eval([A,B,C,D,P],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D,P],sol),list) fi;
 fi; fi;
od: od: od: od:
convert(Sol,list);
end proc:

Pic:=proc(L,Sol)
local P1, P2, T;
uses plots, plottools;
P1:=seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L):
P2:=line(Sol[1],Sol[3],color=red, thickness=2), line(Sol[2],Sol[4],color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"],[Sol[5][],"P"]], font=[times,18], align=[left,above]);
display(P1,P2,T, scaling=constrained, size=[800,500], axes=none);
end proc: 

Examples of use:

L:=[[[0,0],[1,4]],[[1,4],[6,7]],[[6,7],[12,0]],[[12,0],[0,0]]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L, Sol);

# Check (areas of all 4 regions)
Area([[L[1,1],Sol[4],Sol[5],Sol[1],L[1,1]]]),
Area([[Sol[4],Sol[5],Sol[3],L[4,1],Sol[4]]]),
Area([[Sol[5],Sol[2],L[3,1],Sol[3],Sol[5]]]),
Area([[Sol[5],Sol[2],L[1,2],Sol[1],Sol[5]]]);

        

 

L:=[[[1+cos(-t),1+sin(-t)],t=-3*Pi/2..-Pi],[[0,1],[-1,0]],[[cos(t),sin(t)],t=Pi..2*Pi]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L,Sol);

    

# The boundary is the Archimedes spiral and the arc of a circle

L:=[[[t*cos(t),t*sin(t)],t=0..2*Pi],[[Pi+5*cos(-t),sqrt(25-Pi^2)+5*sin(-t)],t=arccos(Pi/5)..Pi-arccos(Pi/5)]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

     

 

L:=[[[0,0],[2,0]],[[2,0],[1,sqrt(3)]],[[1,sqrt(3)],[0,0]]]:
Sol:=evalf[5](Into_4_Equal_Areas(L, AllSolutions, 0.1)); # All 3 solutions
plots:-display(<Pic(L, Sol[1]) | Pic(L, Sol[2])  | Pic(L, Sol[3])>, size=[300,300]);  

 

L:=[[[-t,-sin(-t)],t=-5*Pi/4..0],[[cos(t),sin(t)-1],t=Pi/2..3*Pi/2],[[t,cos(t)-3],t=0..3*Pi/2],[[3*Pi/2,-3],[5*Pi/4,sqrt(2)/2]]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

More examples can be found in the attached file.

4_Equal_Area1.mw

[Edit]. The post has been edited. One inaccuracy in the code has been corrected, which could sometimes lead to errors. Two options have been added to the code of Into_4_Equal_Areas procedure. The first option is the formal argument  N . If N=OneSolution  (by default), the procedure returns one solution. If  N=AllSolutions , the procedure returns all solutions that it can find. The  eps  option has also been added (by default, eps=0). It is advisable to use it when we are looking for all solutions, and the ends of the cutting segments fall on the boundaries of intervals (this option slightly expands the boundaries of intervals, otherwise the  fsolve  command sometimes misses solutions). Two new examples have also been added.

 

Keywords: Intermediate axis theorem, Tennis racket theorem, Dzhanibekov effect, Coriolis force, Euler equations

In 1988 I witnesses the instability of the rotation about the intermediate axis of a foam brick.

Since then I have been fascinated by this effect. It was one of the many experiments which enriched a lecture series on kinetics and on that day Euler equations were on the agenda. Colored surfaces of the brick made it possible to observe the effect without micro gravity and slow-motion equipment.

This post is about reproducing an “intuitive” visualization of an explanation of the effect by Terry Tao from 2011 using 4 rigidly connected point masses. 8 years later the explanation was animated in a YouTube video (The Bizarre Behavior of Rotating Bodies) and considered to be the “best intuitive” explanation.

Motivated by the video, I wondered whether a similar animation with acting forces is possible with MapleSoft products and whether there might be a better intuitive explanation without the use of centrifugal forces. Initially I saw this more as a good test of MapleSim’s visualization capabilities. Finally, it took over 3 years and numerous attempts (mostly during vacation, kind of a substitute for drawing circles in the sand...) to come to a conclusion on the effect.

Intermediate_axis_theoreme_with_3_point_masses.msim

About the model:

Unlike the YouTube video, I decided to simulate 3 identical point masses because a 3-mass model fits better to a T-handle (overlayed in the animation above), video footage from space experiments and discussions in this forum (221298, 225760, 228066).

The movement of the model generates acceleration forces on each mass. The clip displays the corresponding opposing forces that act in the model (i.e. act on the massless T-structure). The blue mass, which is not perfectly centered on the axis of rotation at the start of the simulation perturbs the orbits of the red and the green masses. That was my initial intuitive attempt to explain the effect.

The 3 masses form an isosceles triangle. Here it is helpful to think of a rotating arrowhead where the shape determines stability of the rotation. The aspect ratio (the ratio of the height to the base length) of the triangle determines the stability of rotation about the mirror symmetry axis of the triangle (i.e. the symmetry axis of the T-structure). An obtuse triangle (“blunt”, aspect ratio < sqrt(3)/2) is unstable when rotating about an axis that is slightly inclined with respect to this axis of symmetry. The inclination can be in the plane of the triangle or out of plane. An acute (“pointy”) triangle only wobbles.

About the MapleSim model:

A supplementary rigid body component without mass and rotational inertia is used at the center of mass of the three masses to impose initial conditions. Rotating the triangle at the start of the simulation about the center of mass of the 3 masses prevents the triangle from drifting laterally away from its initial position. This effect of lateral drift is visible in video footage from space with the T-handle.

The rotational inertia of the other rigid body components is set to zero. Without rotational inertia it could be assumed that only Newtonian mechanics are used in the simulation (i.e. no Euler equations are integrated). This is however wrong. MapleSim generates automatically from a system with 3x6=18 coordinates a system with 3 Newtonian equations for translation and 3 Euler equations for rotation.

Forces and moments are measured with sensor components. Visualization is done with force and moment visualization components. These components are “abused” to display the following other physical quantities:

The angular momentum of the masses

The vectors of the angular velocity and the angular acceleration

Moments of the forces with respect to the center of mass

Moments of the forces with respect to the center of the base of the triangle

For a clean model, sensor components and mathematical components to calculate physical quantities are grouped in three subsystems (one per mass, indicated with a colored dot in the image below).

The model contains parameter sets for in plane and out of plane inclination of the axis of the T with respect to the initial axis of rotation (the x-axis).

 

Visualization of physical components can be turned on by enabling the corresponding subsystems which are labeled accordingly (in the image above the display of the angular momentum is enabled). The subsystem “Verification” computes quantities that should either be conserved or should be equal to zero.  Calculation of quantities is done with MapleSim’s mathematical components (i.e. no embedded code or custom components are used).

 

Some observations

Kinetic energies are exchanged between the masses.  During a flip of the T (see animation above), the red and green masses “exchange” their energy. The blue mass mediates this exchange.  Depending on the initial conditions (in plane or out of plane), the energy of the red mass decreases first during the flip and the energy of the green mass increases (and vice versa, as seen below for the out of plane case which exhibits symmetric energy distributions).

Energy peaks are a good measure for the flip frequency. The frequency increases with the initial misalignment of the rotation axis to the symmetry axis of the T.

 

Tracing the blue perturbing mass reveals that the mass never gets closer to the (initial) rotation axis than its initial off-axis position.

The angular momenta of the masses vary, but the total angular momentum is, as expected, conserved. In the image below the angular momenta of the three masses are visualized to the left. The change of kinetic energy can be appreciated from the change in magnitude of the angular momenta.

The vector of the angular velocity (violet, at the origin) wobbles during the flip but does not flip direction. The vector of the angular acceleration (orange) rotates in the yz-plane

Forces act in the plane of the triangle. There is no component normal to the plane, as in the YouTube video, that could cause a flip. Thus, the displayed forces measured in the inertial reference frame do not provide an intuitive explanation why the flip occurs.

The same applies for the moments of the forces at the center of mass: They are perfectly balanced. There is no net component that could be attributed to an in-plane rotation.

 

Why are the animations different: Apparent vs. internal reactive forces.

The MapleSim animation shows internal reactive forces that illustrate the interplay of the moving masses which are bound to each other. They act in the model and obey actio = reactio, which means that the same vectors of opposite sign pull on the masses when the masses are isolated (they follow Newtons second law and equate to mass times the vector of acceleration; the last image in this post displays an isolated mass and the opposing force). 

On the contrary, the YouTube animation shows apparent forces (centrifugal forces) that appear when accelerations are described in a reference frame that moves (accelerates or rotates) with respect to the inertial reference frame. They look like external forces acting on the model, but they are not real. Since apparent forces are fictitious (not real), not everyone is satisfied with using them for an intuitive explanation.

 

Can the MapleSim animation be improved?

Calculation of apparent forces is possible but less straight forward for the simple reason that the Mathematical components library does not provide operators for coordinate transform and matrix multiplication. Those operators are normally not required for simulation purposes. (It would be interesting to see how calcualtion of apparent forces can be done in MapleSim. Verification of code implementation might not be as easy as in the inertial reference frame.)

What ultimately prevents a reproduction of the video is the observer/camera view that rotates with the model. This feature does not exist in the current version of MapleSim 2024. To reproduce the video, Maple has to be used. This would also make the implementation of the calculation of apparent forces much easier as compared to, for example, Modelica code implementation (at least for me).

 

Is the 3-mass model equally intuitive as a 4-mass model?

The initial idea was to have two orbiting masses that are perturbed by a third mass. The third mass flips like a pointer back and forth while the two masses still follow their orbit. This is in case of 3 identical masses only possible with a short-legged T as shown here:

Only a reduced mass would allow for a longer leg. Since the T has only one axis of symmetry, the two orbiting masses do not orbit in a plane. They perform a wobbling motion and shift laterally in position during a flip since the rotation is performed about the common center of mass. Only when 4 masses are used in a symmetrical cross configuration, two masses can orbit closer to a plane that contains the common center of mass while the two perturbing masses flip sides of the plane (the wobble is less pronounced but still visible by the enlarging blue trace in the animation below).

With a mass ratio of 1:100 in the animation below the two orbiting masses create kind of a centrifugal potential field in which the two perturbing masses swing like a pendulum. In this configuration the two perturbing masses can no longer be regarded as strongly disturbing, but rather as oscillating satellites. The sudden flip is created by the increasing accelerating field strength which increases with the distance from the axis of rotation. This lets a pendulum swing with a stronger than expected acceleration and is perhaps a new insight.

Both models represent the simplest possible implementation to generate the effect in terms of number of parameters. The 4-mass configuration has more objects but is simpler to understand because of the higher degree of symmetry.  Either identical masses at varying distances or identical distances at varying masses can be used in both models. No more reduction of parameters is possible to generate the effect. A two mass object cannot even wobble.

Out of plane initial inclination makes the acceptance of an explanation easier since the orbiting masses do not generate a momentum as in the case of an in plane inclination. For the latter case an intuitve explanation is more difficult and perhaps there is none.

Although the pendulum swing of the out of plane case might provide an intuitive explanation of the effect it is not fully satisfying. It does not explain why larger masses than the orbiting masses do not lead to a swing but smaller masses do. Another well-made video provides an explanation for that.