Referring to the screenshots, "J" can be converted to "N m" in MF2024.1, but not in M2024.2.
Is this some sort of bug in M2024.2?

Just around a month after the first release I am glad to announce the second public release of this project.

Changes from last release:

• added angled cuts to beam ends
• fixed bug for bolt connections with thick steel plates
• rewrite check if fasteners are placed within beams
• removed some obsolete procedures in NODEFunctions
• minor changes in XML file headers

For more information see https://github.com/Anthrazit68/NODEMaple.

With the new release of Maple Flow 2024.2 the units "Area" and "Speed" don't work.

I run a MaxBook Pro with macOS Sequoia 15.2 and uninstalled MF2024.2.

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.

This newer video also gives an explanation why internal forces must act in the plane of the rotating object but does not display them in the animation. I guess this is because the introduction of real forces would have spoiled the intuitive explanation of the video. Isolating a mass and adding an internal force now as an external force leads to an equivalent system that reproduces the effect of the rotating object. If the same force is applied in the opposite direction on the isolated mass, the isolated mass moves along the same trajectory.

4_lumped_masses_and_one_single_force_driven_mass.msim

Isolating only one mass breakes the symmetry of the model. It also gives the false impression that the introduced perturbing force acts primarily on the opposite mass. A 3-mass model does not lead to such a false interpretation. By isolating the opposite mass and introducing a second perturbing force, the discussion shifts more to the analysis of the wobble and the rotational acceleration of the orbiting masses and less to the flip.

In summary, internal forces describe how the masses interact but their orientation is counterintuitively perpendicular to direction of the flip. On the other hand, centrifugal forces that we intuitively assume acting in a 4-mass model from the perspective of an observer from an inertial reference frame do not exist. This assumption provides an intuitive explanation which is physically wrong. In the same way an accelerating radial force field does not exist. Mathematically and physically correct is a description from a rotating observer which uses fictious forces.

For me both intuitive explanations of the videos are somehow useable, but both involve centrifugal forces (in one case explicitly and in the other wrongly assumed by an observer). This is not satisfying when the goal is not to use fictious forces.

Conclusion

MapleSim visualization components can be used for more than displaying forces and moments. They are very helpful to better understand physical phenomena.

A camera view observer on a rotating reference frame would have made observation of the direction of the internal forces much easier and might have given more insights. As of now, Maple is required to reproduce the animation in the video.

There is no better intuitive visualization/explanation with a model of 3 identical masses. A 4-mass configuration provides better insight but does not explain all.

In reality every freely rotating object with more than two point masses inevitably wobbles.

Hello,

I present the result of work on a new project in the field of classical mechanics. It is a grateful and interesting topic that gives a lot of satisfaction. I am attaching the Maple worksheet.

Best regards

Rolling_Disk_3D_on_x0y_plane.mw

I am a new user of Maple Flow 2024.2.

Since I installed this version I got trouble with the following commands:

solve(x^3-2x^2+3x-2)=1.00.  Just 1 root is returned

fsolve((x-2)(x+3)(x-1))=-3.  Just 1 root is returned

ifactors(3024)=.   Maple Flow latch in and crashes without errot massage

seq(i^2,i=1..5)=. Sequence not executed

subs(...)= Substitute not executed

Optimization:-Minimize(...)=.  Latch in, error

With the help of the Maplesoft-Team I uninstalled and installed several times MF on a MacBook Pro Sequoia 15.2 and on a MacBook Pro Ventura 13.7.2 with and without Firewall and McAfee. 
No success, the problems still remain.

I'll no longer use Maple Flow in this version.
I expect a new update asap!

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.

Maple Transactions Volume 4 Issue 4 has now been published.

This issue has two Featured Contributions by people who have been plenary speakers at Maple Conferences in the past, namely Veselin Jungić and Juana Sendra. We hope you enjoy both articles.  There is an accompanying video by Professor Sendra, which we will add a link to when it becomes ready.

As usual, there is an article in the Editor's Corner, but this one is a bit different.  In this one, Michelle Hatzell (the new copyeditor for Maple Transactions, who is also a Masters' student working with me at Western) and I have written about a fun use of Maple's colour contour plots to make an image that might be used as the cover of an upcoming book, namely Perturbation Methods using backward error, which I'm just finishing now with Nic Fillion and which SIAM will publish next year.  So, while there's some math in that paper, it's more about Maple's utilities for colour plotting; so you might find it useful.  We also hope you like at least some of the images.  Some are more attractive than others!

We have several Refereed Contributions, not all of which are ready at this time of publication but which will be added as they are revised and sent in.  We have a nice paper on using continued fractions in a high school context, another on code generation, and another on using Digital Signal Processing in Engineering courses.

Finally we have a first publication in French, by Jalale Soussi.  Actually we have the paper also in English: we chose to publish both, in our Communications section, each with links to the other.  It is possible to publish in Maple Transactions solely in French, of course, but the author provided both, so why not?

Happy reading, and best wishes for 2025. 

Download dynamics_of_novel_disease_outbreak.mw

I think a new integer subtype is needed: integer greater than one, gtoint.

isgto := proc(x::anything)
  local X;
  X:=x;
  return type(X,integer) and (X>1);
end:

AddType(gtoint,z->not isgto(z));

Hello,

This is a simple dynamics example illustrating the behavior of a disc as a pendulum under ideal conditions. Source attached.

Simple_Disk_Pendulum.mw

I am very happy to announce the first public release of a project which I have been working on for the last couple of years.

NODEMaple consists of a set of Maple workbooks and a library for structural design based on the Eurocode.

Currently the main development of the workbooks is focused on "Eurocode 5: Design of timber structures" with the Norwegian Annex.

This software has been made public in the hope of that it might be useful for other structural designers, professionals as well as students. Everyone interested is very Welcome to contribute to this project. The code is published under the GPLv3 license.

For more information see https://github.com/Anthrazit68/NODEMaple.

I didn't put it in the title, but of course this is a post about Advent of Code, in particular Days 16 and 18 which feature a perenial favorite type of problem: finding shortest paths in mazes.

Your input for these is always a maze given as an ascii map.  Like so:

###############
#.......#....E#
#.#.###.#.###.#
#.....#.#...#.#
#.###.#####.#.#
#.#....#....#.#
#.#.#####.###.#
#...........#.#
###.#.#####.#.#
#...#.....#.#.#
#.#.#.###.#.#.#
#.....#...#.#.#
#.###.#.#.#.#.#
#S......#.#...#
###############

There's lots of ways to import one of these into Maple and then solve the maze, but I am to highlight how to do it with GraphTheory.  I am going to start with a GridGraph and then remove the walls in order to leave a just the vertices that represent the paths:

with(StringTools): with(GraphTheory):
maze:=
"###############
#.......#....E#
#.#.###.#.###.#
#.....#.#...#.#
#.###.#####.#.#
#.#....#....#.#
#.#.#####.###.#
#...........#.#
###.#.#####.#.#
#...#.....#.#.#
#.#.#.###.#.#.#
#.....#...#.#.#
#.###.#.#.#.#.#
#S......#.#...#
###############
":
mazelines := (Split(Trim(maze), "\n")):
sgrid := ListTools:-Reverse((map(Explode, mazelines)) ):
m,n := nops(sgrid), nops(sgrid[1]);
tgrid := table([seq(seq([i,j]=sgrid[i,j],i=1..m),j=1..n)]):
start := lhs(select(e->rhs(e)="S", [entries(tgrid,'pairs')])[]);
finish := lhs(select(e->rhs(e)="E", [entries(tgrid,'pairs')])[]);

Now the maze is stored in the table tgrid, and it is easy to find the walls and paths.  In a GridGraph the vertices are labeled with their coordinates as "x,y" and so we rewrite our list of paths in that form, so we can create the induced subgraph of the Grid that includes only those vertices.

(walls,paths) := selectremove(e->rhs(e)="#", [entries(tgrid, 'pairs')]):
paths := map(s->sprintf("%d,%d",lhs(s)[]), paths):
H := SpecialGraphs:-GridGraph(m,n);
G := InducedSubgraph(H, paths);

We can use StyleVertex to highlight the start and finish.

StyleVertex(G, sprintf("%d,%d",start[]), color="LimeGreen");
StyleVertex(G, sprintf("%d,%d",finish[]), color="Red");

plots:-display(<
DrawGraph(H, stylesheet=[vertexshape="square", vertexborder=false, vertexcolor="Black"], showlabels=false) | 
DrawGraph(G, stylesheet=[vertexshape="square", vertexborder=false, vertexcolor="Black"], showlabels=false)>);

(I omitted a step where I set the vertex locations of the maze grid, you can see that in the attached worksheet)

Now finding a path through the maze is as easy as calling GraphTheory:-ShortestPath

sp := ShortestPath(G, sprintf("%d,%d",start[]), sprintf("%d,%d",finish[]) ):

StyleVertex(G, sp[2..-2], color="Orange");
StyleEdge(G, [seq({sp[i],sp[i+1]}, i=1..nops(sp)-1)], color="Orange");
DrawGraph(G, stylesheet=[vertexshape="square", vertexpadding=10, vertexborder=false,
             vertexcolor="Black"],  showlabels=false, size=[800,800]);

Now, Advent of Code seldom gives you a completely simple maze like this, often these is a twist like having to calculate the costs of turns seperately from the cost of steps, or each direction or position has a seperate cost associated with it.

For example, Day 16 has us starting facing east, and then turns cost 1000, while moving forward costs 1. That sort of problem is no longer exactly a maze, instead of the vertices being representing an "x,y" position, instead you increase the number of vertices by a factor of 4, so that you have a vertex for every position and orientation "x,y,o" with edges of weight 1 between adjacent vertices with the same orientation and edges of wieght 1000 to connect "x,y,N" to "x,y,E" and "x,y,W" e.g.  In that sort of weighted graph, we can use GraphTheory:-DijkstrasAlgorithm to find the shortest path and it's weighted cost.

In this code, we expand our list of maze locations with directions, and the use the grid table to generate a list of weighted edges:

dtable := table([0=[0,1], 1=[1,0], 2=[0,-1], 3=[-1,0]]):
dname := table([0="N",1="E",2="S",3="W"]):
dpaths := map(s->local d;seq(cat(s,",",d), d in ["N","E","S","W"]), paths):

edges := NULL:
for i from 1 to m do for j from 1 to n do
    if tgrid[[i,j]] = "#" then next; end if;
    for d from 0 to 3 do
        dir := dtable[d];
        if tgrid[[i,j]+dir] <> "#" then
            edges := edges, [{cat("",i,",",j,",",dname[d]), cat("",i+dir[1],",",j+dir[2],",",dname[d])},1];
        end if;
        edges := edges, [{cat("",i,",",j,",",dname[d]), cat("",i,",",j,",",dname[d+1 mod 4])}, 1000],
                 [{cat("",i,",",j,",",dname[d]), cat("",i,",",j,",",dname[d-1 mod 4])}, 1000];
    end do;
end do; end do:

Gd := Graph(dpaths,weighted,{edges});

Once that is done, it's a simple matter of calling Dijkstra's Algorithm on the graph, but notice that we can reach the finsh while traveling north or east, so we need to find the sortest path to both (you can pass a list of vertices to Dijkstra, and it will efficiently calculate paths to all of them), and select the smaller of the two:

spds := DijkstrasAlgorithm(Gd, cat("",start[1],",",start[2],",E"), 
    [cat("",finish[1],",",finish[2],",N"), cat("",finish[1],",",finish[2],",E")] , 
    distance):
i := min[index](map2(op,2,spds)):
spd := spds[i];

spd := [["2,2,E", "3,2,E", "4,2,E", "4,2,N", "4,3,N", "4,4,N", "4,5,N", "4,6,N", "4,6,E", "5,6,E",
 "6,6,E", "7,6,E", "8,6,E", "8,6,N", "8,7,N", "8,8,N", "8,9,N", "8,10,N", "8,11,N", "8,12,N", 
"8,12,W", "7,12,W", "6,12,W", "5,12,W", "4,12,W", "3,12,W", "2,12,W", "2,12,N", "2,13,N", 
"2,14,N", "2,14,E", "3,14,E", "4,14,E", "5,14,E", "6,14,E", "7,14,E", "8,14,E", "9,14,E", 
"10,14,E", "11,14,E", "12,14,E", "13,14,E", "14,14,E"], 6036]

We can then plot to compare this to the unweighted shortest path:

dsp := ListTools:-MakeUnique( map(s->s[1..-3], spd[1]) );
StyleVertex(G, dsp[2..-2], color="DarkBlue");
StyleEdge(G, [seq({dsp[i],dsp[i+1]}, i=1..nops(dsp)-1)], color="DarkBlue");

DrawGraph(G, stylesheet=[vertexshape="square", vertexpadding=10,
             vertexborder=false, vertexcolor="Black"],  showlabels=false,
          size=[800,800]);

And you can see it's a path that requires more steps, but definitely uses fewer turns if we start facing east/right (6 vs. 9):

I hope this has given you a little bit of a flavor of how to use GraphTheory commands to solve path finding problems.  Like with the second part here, usually the biggest challenge is figuring out how to encode and construct a graph that represents your problem.  Then the actual commands to solve it, are easy. You can see all the code, and a couple steps I left out from above in this worksheet: Mazeblog.mw

And just for fun, here's a Maple workbook that imports a maze from an image and solves it: MazeFromImage.maple

About once a year Advent of Code give you a problem that is a gift if you are using a computer algebra system. This year that day was Day 13. The day 13 problem was one of crazy claw machines.  Each machine has two buttons that can be pressed to move the claw a given number of X and Y positions and a prize at a given position. A buttons cost 3 tokens to press, and B buttons cost 1 token to press, and we are asked to find the minimum number of tokens needed, IF it is possible to reach the prize. The data presented like so:

Button A: X+94, Y+34
Button B: X+22, Y+67
Prize: X=8400, Y=5400

Button A: X+26, Y+66
Button B: X+67, Y+21
Prize: X=12748, Y=12176

...

Some times the input is harder to parse than the problem is so solve, and this might be one of those cases.  I tend to reach for StringTools to take the input apart, but today the right tool is the old school C-style sscanf (after using StringSplit to split at every double linebreak).

machinesL := StringTools:-StringSplit(Trim(input), "\n\n");
machines := map(m->sscanf(m,"Button A: X+%d, Y+%d\nButton B: X+%d, Y+%d\nPrize: X=%d, Y=%d"), machinesL):

Now we have a list of claw machine parameters in the form [A_x, A_y, B_x, B_y, P_x, P_y] and we need to turn those into equations that we can solve. We want the number of A presses a, and B presses b to get the claw to the P_x, P_y position of the claw, it is simple to just write them down:

for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
end do;

Now because of the discrete nature of this problem, we need our variables a and b to be non-negative integers.  When solving this, I first reached for isolve like this:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol <> NULL then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Now, sometimes Advent of Code inputs contain a lot of hidden structure.  I wrote the code above, it worked on the sample input, so I tried it immediately on my real input (about 300 claw machines like the above) and IT WORKED.  But, you might notice that this code does not deal with a couple cases that could have appeared.  In particular, it doesn't check that the solutions are positive.  It also doesn't handle cases where there is more than one possible solution.  The former is easy to check

if sol <> NULL and eval(a,sol) >= 0 and eval(b,sol) >= 0 then

Unfortunately isolve does not handle inequalities, but you could try with solve, but it doesn't save us any checking, because we'd still have to check if the solutions are integers, so we might as well have just solved the equation and then checked if it were a nonnegative integer.

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := solve(eqn);
   if type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;
ans1 = tokens;

In the multiple solution case we get something like {a = 3 - 2*b, 0 <= b, b <= 3/2} which has some great information in it but might be hard to handle programmatically, so let's see what isolve does with those cases to see if it's easier to deal with

> eqn := { 17*a + 84*b = 7870, 34*a + 168*b = 15740 }:
> constr := { a >= 0, b>= 0 }:
> sol := isolve(eqn);
               sol := {a = 458 - 84 _Z1, b = 1 + 17 _Z1}

> constr := eval({ a >= 0, b>= 0 }, sol):
            const := {0 <= 1 + 17 _Z1, 0 <= 458 - 84 _Z1}

> obj := eval(3*a+b, sol);
                         obj := 1375 - 235 _Z1

You can see it's easy to tell if these show up in your input, since your "token" total will have the _Zn variables in it.  Now, since everything is simple and linear here, it seems like you could use solve to find the rational value of _Z1 that makes obj=0 and then take the closest integer but it's not so simple, we actually have to deal with the contraints that a and b be positive too. So, it really just makes sense to bring out the big hammer of Optimization:-Minimize which allows us to directly optimize over just the integers.  So a full solution looks like this:

tokens := 0:
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol = NULL then
      next;
   end if;
   constr := eval({ a >= 0, b>= 0 }, sol);
   obj := eval(3*a+b, sol);
   if not type(obj, constant) then
      tokens := tokens + Optimization:-Minimize( obj, constr, assume=integer )[1];
   elif andmap(evalb, constr) then
      tokens := tokens + obj; 
   end if;
end do;

But since we're bringing out the big hammer, why not just use Optimization in the first place.  The main reason is that Minimize doesn't simply return NULL when it doesn't work, instead it throws an exception, so we need to find all the exceptions that can occur and handle then with a try-catch, thus:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   try 
       sol := Optimization:-Minimize(3*a+b, eqn, 'assume'='nonnegint')[1];
       tokens := tokens + sol;
   catch "no feasible":
   end try; 
end do;
tokens;

(you can in fact omit the string in the catch: statement, but I can tell you from long experience that that is an excellent way to make your code much much harder to debug)

Alright, so how did people not using Maple solve this problem?  The easiest way to solve it, and the one used by all the cheaters scraping the website and using LLM-based code generators that auto-submit solutions to get into the Top 100, was to just check all possible a, b values in 0..100 and take the values than minimize 3*a+b when reaching the prize coordinates.  That's only feasible because the problem states the 100 is an upperbound for a and b, but it's also very fast (about 1/10 second in Maple):

tokens := 0:
for m in machines do;
sol := infinity;
for i from 0 to 100 do for j from 0 to 100 do
    if i*m[1]+j*m[3]=m[5] and i*m[2]+j*m[4]=m[6] and 3*i+j < sol then
        sol := 3*i+j;
    end if;
end do; end do;
tokens := tokens + ifelse(sol=infinity,0,sol);
end do;

It does not scale at all to part 2 (which modified everything to be bigger by about 10 trillion), and it seems that foiled all the LLM solvers. So, what solutions scaled in languages without integer equation solvers?  Well, the easiest solution is just to solve the general equation using paper and pencil

And you can just hard code that formula in, check that it gives integer values and compute the tokens. As long as you get unique solutions, that looks something like this

solveit := proc(m)
local asol := m[4]*m[5] - m[3]*m[6];
local bsol := m[1]*m[6] - m[2]*m[5];
local deno := m[1]*m[4] - m[2]*m[3];
if deno = 0 then return -2^63; end if; # multiple solution case - not handled
if asol mod deno = 0 and bsol mod deno = 0
   and (   ( deno>=0 and asol>=0 and bsol>=0 ) 
        or ( deno<=0 and asol<=0 and bsol<=0 ) )
then

    return 3*asol/deno+bsol/deno;
else
    return 0;
end if;
end proc:

Which if you have this in Maple, you can impress your friends by auto generating solutions in other languages. Here, for your FORTRAN friends

> CodeGeneration:-Fortran(solveit);

Warning, the following variable name replacements were made: solveit -> cg
       integer function cg (m)
        doubleprecision m(*)
        integer asol
        integer bsol
        integer deno
        asol = int(-m(3) * m(6) + m(4) * m(5))
        bsol = int(m(1) * m(6) - m(2) * m(5))
        deno = int(m(1) * m(4) - m(2) * m(3))
        if (deno .eq. 0) then
          cg = -9223372036854775808
          return
        end if
        if (mod(asol, deno) .eq. 0 .and. mod(bsol, deno) .eq. 0 .and. (0
     # .le. deno .and. 0 .le. asol .and. 0 .le. bsol .or. deno .le. 0 .a
     #nd. asol .le. 0 .and. bsol .le. 0)) then
          cg = 3 * asol / deno + bsol / deno
          return
        else
          cg = 0
          return
        end if
      end

Another way that you might solve this without solve is to use a linear algebra library to solve the linear system.  It works even if you only have a numeric solver, but you have to be careful about checking for integers:

tokens := 0:
for m in machines do
    sol := LinearAlgebra:-LinearSolve(
               Matrix(1..2,1..2,[m[[1,3]],m[[2,4]]], datatype=float), 
               Vector(m[5..6], datatype=float));
    if abs(sol[1]-round(sol[1])) < 10^(-8) and abs(sol[2]-round(sol[2])) < 10^(-8)
       and sol[1] >= 0 and sol[2] >= 0
    then
       tokens := tokens + 3*sol[1]+sol[2];
    end if;
end do;

Finally, a lot of people solved this sort of thing with the Z3 Theorem prover from Microsoft research which is also way more than you need, but it mostly just uses SMTLIB, which we also have in a library for in Maple, and it can just be used in place of solve

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := SMTLIB:-Satisfy(eqn) assuming a::nonnegint, b::nonnegint;
   if sol <> NULL and type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Notice that Satisfy handled the multiple solution case just by choosing one of the many solutions. It is possible to get SMTLib to optimize but it is slightly more involved, and this post is already too long. This time, I've put all this work in worksheet: Day13-Primes.mw

Major deficiency in Physics[Vectors]; Distinct sets of basis vectors are not recognized!

You can't define vectors in alternative bases like: {\hat{i}',\hat{j}',\hat{k}'} or {\hat{i}_{1},\hat{j}_{2},\hat{k}_{3}}.

This deficiency has been around for a while. I have found other posts regarding this problem.

The deficiency greatly reduces the allowable calculations with Physics[Vector].

Are there any plans to fix this?

Here is my example which shows this deficiency in more detail.

physics_vectors_and_multiple_unit_vectors.mw
 

Download physics_vectors_and_multiple_unit_vectors.mw