Approximation_of_ODEs_with_Cubic_Splines.mw

Download New_in_Physics_2025.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions

Download Einstein_Equations_From_First_Principles.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions

I’m absolutely delighted to announce the launch of Maple 2025!

Although you see a new release every year, new features take anything from a few fast-paced weeks to develop, to months of careful cultivation.

Working on so many features in parallel, each with varying time scales, isn't easy! We have to fastidiously manage and track our work.

So it's easy to lose ourselves in the daily minutiae of software development. To help us maintain perspective, we constantly ask ourselves questions like:

• What user problem are we solving and how often does this problem occur?
• Can we validate our proposed solution with preliminary user feedback?
• Is this a solution to a problem that doesn't exist and will never exist, or are we pre-empting a future need?
• Are we offering value to our users?

Given the answers, we course-correct to make sure we stay on track for our central mission - to make you happy, and to keep you coming back year-after-year.

With Maple 2025, I think we've smashed that goal. We have many new features that'll appeal to many different types of users - from students, educators and mathematicians, to engineers, scientists and technical professionals

Let me walk you through some of my personal highlights.

It’ll be difficult for anyone to miss this - Maple 2025 has a new interface! It’s a ribbon-based UI that look clean and contemporary, and helps you find and discover tools more quickly than before.

You have large, meaningful icons.

Items are logically grouped.

The ribbons is contextual. If you click on a plot, you get new tabs for interacting with and drawing on the plot.

A new Education tab collects pedagogical resources that were scattered around the interface in prior releases.

This is the biggest visual overhaul to Maple in many years. We hope you like it! 

We also appreciate that changes in look and feel can be divisive. Please rest assured that we will refine and finesse the interface with each successive release; your comments and suggestions are most welcome.

The new interface is available on Windows and Linux, and as a technology preview on Mac.

The right arrow key on my keyboard is wearing out…and it’s all because of Maple. I’m knee deep in Maple nearly every day entering equations, and I’m always using right-arrow to move the cursor. It gets kind of tedious!

This anecdote reflects some investigative work we did. We comprehensively examined our internal library of thousands of Maple worksheets and discovered that these three input patterns are extremely common.

Previously, you’d use the right-arrow key to move the cursor out of the exponential, division or subscript.

Now, in Maple 2025, when you

• type ^, /, or enter a literal subscript with a double-underscore,
• followed by a number or symbol
• and then input another operator (such as +)

the operator is automatically inserted on the baseline (except when y = 1).

Of course, you can also make the cursor return to or stay in the exponent or denominator with a simple keystroke, when that is what is needed.

This is one of those little quality of life refinements that I’m very fond of - it’s a little visual and usability dopamine hit.

The sum command (and its typeset form) now indexes into vectors without you needing to spam unevaluation quotes all over your expression.

We’ve been integrating units deeper into the Maple system, release after release. Much of this is driven by our engineering users.

A few releases ago, we made int(numeric) compatible with units. With Maple 2025, you can now numerically differentiate  expressions and procedures that have units.

I’m a grizzled thermodynamics hack, so here’s an example in which I calculate the specific heat capacity of water by differentiating enthalpy with respect to temperature (and then confirm the result with the built-in value):

This is in addition to many other improvements to the units experience.

Although this is a part of Maple that I don’t touch often (my colleague Karishma takes point on the education side), I REALLY wish I’d had this when I was struggling with math.

You can now automatically generate unlimited variants of the same problem for students to solve with the Try Another feature, which has been added to Maple’s Check My Work tools (another feature I really could have used!). This is available for many common math principles, including factorization, simplification, integration and more.

This is just one of the improvements in Maple 2025 for teaching and learning.

 If you’ve ever found yourself going back and forth (and back and forth) between two large, almost identical-looking Maple expressions, trying to figure out how they are different, you’re going to love this one.  ExpressionTools is a new package that lets you compare the differences between two expressions.

I really like the use of color to highlight differences. Less squinting at the screen!

You can now run Maple Flow worksheets from Maple (you don’t need Maple Flow installed to do this). You can send parameters into the Flow worksheet and extract your desired results.

This means you can use the entire flexibility of Maple to analysis and manipulate your Flow worksheet. You could, for example:

• Attach a Flow worksheet to a Maple workbook and create an interactive application
• Carry out parameter studies of a Flow worksheet by evaluating it over many parameter sets in Maple
• Create an Excel interface for a Flow worksheet using the Maple add-in for Excel

Simplify is one of those functions that literally tens of thousands of people use each day. Every time we make an incremental improvement, the cumulative benefits across our entire user base are significant.

We’ve refined simplify in a number of critical ways. For example, simplify now recognizes when exponentials can be profitably converted to hyperbolic trig functions:

The analysis of many scientific phenomena result in Laplace transforms that do not have a symbolic inverse which can be expressed in terms of elementary functions. This includes applications in heat transfer, fluid mechanics, fractional diffusion processes, control systems and electrical transmission.

For example, this monster Laplace transform results from an analysis of voltage on a transmission line:

You can now numerically invert this transform courtesy of an enhancement to inttrans:-invlaplace - a fast quadrature method.

I’ve saved what I think has the most future potential for last.

I’m sure nearly all of you have experimented with the various AI tools. They’re an inevitable part of our present and future, whether we're comfortable with it or not.

This is something we've been mulling over for some time.

• In Maple 2019, the DeepLearning package made its debut. This package provides tools for machine learning, supporting operations such as classification and regression using neural networks.
• In Maple 2024, we introduced an AI-powered formula lookup feature.

In Maple 2025, we’re giving you an early-stage technology preview of AI-powered document generation.

You can automatically generate worksheet content by prompting an AI, and then gradually refine the content

If you’re an educator, you might want some content that describes applications of calculus. So you might ask the AI “How do I derive the formula for the area of a circle” by entering your prompt into this text box:

This is the worksheet content that may be returned:

If you’re structural engineer who wants to know how to calculate the hardness of concrete, you might ask the AI: “How do I calculate the compressive strength of slow hardening concrete as a function of time? Use the CEB-FIP Model Code 90. Include a worked example with Maple code”.

This worksheet content that could be generated (note the live Maple code):

We’re labelling AI-generated worksheet content as a technology preview. You might see

• text that might be misleading (but sounds plausible)
• code that doesn’t work (but looks plausible)
• or different results each time you click “Generate Document”

For the moment, I would not rely on AI-generated worksheet content without realistic expectations, a healthy dose of scepticism and a modicum of detached analysis. But AI models are rapidly growing in robustness, and we want to position ourselves to best exploit their future potential. The next few years will be VERY exciting.

We can never cover everything in a short blog post like this. So if you want to know more, head on over to the What’s New pages for Maple 2025!

Hi MaplePrimes,
I made a quick code to calculate Aliquot Sequences.

The code works.

I was inspired by a Numberphile YouTube video.

see

aliquot_sequence_cut_1.mw

aliquot_sequence_cut_1.pdf

This should be helpful.

Regards,

Matt

Recently, @zenterix asked a question related to solving the quantum mechanics of the finite-wall-height particle-in-a-box problem. In the last two decades or so I've been teaching a quantum mechanics course every second year, in which I sketch the method of solution of this problem for the class, I but do not get into the mathematical details. A few times I've started to work on it in Maple, but abandoned it because of lack of time. So I decided to persist this time.

This post is more about some of the challenges and tradeoffs in making it work, hence the title. The title is also a nod to the fact that this problem appears in Engel's textbook [1] in a chapter entitled "Applying the particle in a box model to real-world topics". You don't need to know much about quantum mechanics to understand it. From the mathematical point of view you are solving three ordinary differential equations (in three regions) and stitching the solutions together so that the solution, a wavefunction, is continuous with continuous derivative and tends to zero at plus infinity and minus infinity.

The three regions and the wavefunctions for three eigenvalues (energies) are shown here in the figure, which is the final output of the worksheet. Next is the worksheet and I'll comment further on it below.

(It's a quirk of quantum mechanics that we plot two things with different units (wavefunctions and energy) on the same axis, and worse, we offset one by the other.)

Download finite_box_non-dim2.mw

Comments
1. The first general comment is about the issue of how much Maple can do. In principle one should be able to give dsolve the three ODEs and the boundary conditions and it should output the result. We are still far away from that. Even for one region, the boundary conditions at infinity are not handled by dsolve. The other extreme is to solve the ODEs for their general solutions, and follow the algebraic manipulations for implementing the boundary conditions as one might do it on paper. Engel for example has a comment "At this point we notice that by dividing the equations in each pair, the coefficients can be eliminated to give [...]". Maple will not easily do that trick or the manipulations that led to the special form on which the trick is applied. Levine's textbook [2] gives a different non-intuitive set of steps.

But Maple can solve complicated equations, so I wanted a more general strategy that avoids as much as possible the requirement to manipulate into special forms. On the other hand, it is tempting (as @zenterix tried) to just give the three general solutions with 6 arbitrary constants and the boundary conditions to Maple's solve, and solve for the six constants. One finds that all six are zero. This us a consequence of the fact that the ODEs are linear, and misses the crucial point that it is an eigenvalue problem. So there must be some sort of step-by-step guidance for Maple and there is a general but non-obvious strategy to solving such problems.

2. The second general comment is that there is a trade-off between presenting the solution steps without any arcane Maple code, and getting to the solution efficiently and programmatically, without too many manual manipulations. One of my suggestions to @zenterix involved manually dividing both sides of an equations by a fairly complicated expression, to which @acer expressed a preference for programmatic solutions. At the time, I mainly did that to keep close the to existing solution, and was thinking that I usually set things up so that is not necessary. But I found here that I do it a lot, and it is a natural consequence of the interactive way I use Maple. Here's two examples:

In the second line of label (4), I originally got a more complicated form of de_outside, then manually figured out the factor to put inside simplify to get the form you see. I do this frequently with Maple, especially when I use dchange.

In the last line of label (8), I manually divided by B to remove it. Had B been in the denominator, I could have removed it programmatically with numer (an operation I use frequently, though it somewhat recklessly ignores denominator zeroes).

Programmatically doing things can be relatively unreadable and disrupt the readability for the non-Maple reader. For example in label (6) there is "extraneous" code to find which coefficient to set to zero so that the solution goes to zero at +/- infinity. This was forced on me since after generating the nice figure I re-ran the worksheet only to have multiple errors. I originally used eval(outside, {c__1 = B__L, c__2 = A__L}); to change Maple's dsolve constants to read the way I wanted. But I realized that dsolve does not reproducibly assign c__1 and c__2 to the same term each time, raising the general issue of: 

3. Session to session reproducibility. If I am only going to use a worksheet myself, I don't usually worry too much about this, since I can usually debug this fairly easily. On the other hand if the worksheet is for others, it needs to be foolproof. I have taken to always sorting the output of solve, e.g., the assignment to bval in the execution group after label (15). The code added for the dsolve constants works, but is non-trivial Maple (suffixed type testing). (Solving this issue for Eigenvectors is yet more difficult.)

4. Some efforts to make the worksheet more readable were hard to do.

- For some reason, Engel chose to call the two constants in the left region B and B'. Entering B*exp(-k*x) + B'*exp(k*x) in 2D leads to B(x)*exp(-k*x) + diff(B(x), x)*exp(k*x). No doubt this can be fixed, but I didn't pursue this.

- I would normally use upper case Psi for the non-dimensionalized psi, but Psi is the name of a special function in Maple. But wait, now we can use "local Phi;" to solve this problem. However diff(Psi(x),x,x) displays as Psi(2,x), so I had to settle for Phi. (SCR submitted.)

- I build up the left wavefunctions incrementally as expressions assigned to left1, left2, etc, which is rather ugly. Why not use arrow procedures so we can use psi(x), D(psi)(x) etc. Aside from some awkwardness in updating such procedures by redefining them as needed, I ran into the following problem:

limit(A*exp(-k*x) + B*exp(k*x), x = infinity) assuming k > 0; gives as expected
signum(B)*infinity;

but

psi := x -> A*exp(-k*x) + B*exp(k*x);
limit(psi(x), x = infinity) assuming k > 0;
returns unevaluated. What happened to full evaluation?

- I'm used to entering hbar as hbar in 1D; it displays as h with the bar through. In 2-D entering hbar^2 is displayed nicely, but internally it is `ℏ`. We can have the following puzzle:

- I originally wanted a vertical axis label Psi, E/V__0, but labels = [X, Psi, E/V__0] obviously won't work. It took some time to figure out the answer is to make "Psi, E/V__0" an atomic variable. I finally settled on the more honest, if cryptic, label shown.

5. solve gives an "invalid" solution. Originally I worked with solutions left2 etc (see label (6)) after simplify(eval(left2, eqns)). Every second wavefunction was about 10^9 times higher than the others, which was difficult to debug. It turns out the denominator of those wavefunctions was exactly zero but numerically 10^(-10), meaning they plotted as large-scaled versions of the right shape, without any division-by-zero error. The fix is to remove the denominators by scaling (see label (14)).

This is just a consequence of Maple giving generic solutions, and the only way around it is to be aware of it.

Summary

To solve such problems with Maple requires one to play around in an interactive way. After hacking to a solution, it can take some effort to make it presentable and usable to others. But the final result is pleasing, and is certainly much easier than working through the problem on paper. The ability to manipulate complicated expressions means circumventing tricks that might be needed in manual solutions. 

References
[1] T. Engel, Quantum Chemistry and spectroscopy. Pearson 2019, 4th ed. Sec. 5.1, Prob. 5.3. 
[2] I. N. Levine, Quantum Chemistry. Pearson 2009, 6th ed. Sec. 2.4. 

Hi Maple lovers, and others,

I wrote up some simple Maple code that 
checks if an expression (a^3+2) is a prime number.

The output fits on one page of paper.
n_cubed_plus_2.mw

n_cubed_plus_2.pdf

Also, I have an account at mersenneforum.org

Kindest Regards,

Matt

PS Does this help?  The code works.

I'm happy to announce the publication of Volume 5, Issue #1 of Maple Transactions.  You can find it at

mapletransactions.org
 

We have a survey paper by Veselin Jungic and Naomi Borwein on teaching Experimental Mathematics courses as our Featured Contribution.  Many of you will find it interesting and useful.

In the refereed paper section we have a paper on Metaprogramming with Maple and C by Ilias Kotsireas; a paper on fast transposed Vandermonde solving by Hyukho Kwon & Michael Monagan; a paper by David Ulgenes (an honours student in Oslo) on Gamma, Pseudogamma, and Inverse Gamma functions; and a paper by John Campbell on applications of Gosper's nonlocal derangement identity (which, if you don't know that the word "derangement" has a technical mathematical meaning, may give you the wrong impression!).

As usual I've also written something, and I hope you like it: it's about Chladni figures and standing modes in an elliptical drum, and visualizing such in Maple.  It uses Mathieu functions in Maple and noodles a bit about zerofinding (but winds up using fsolve because that's so convenient).
 

Keep the papers coming.  This is the 12th issue of Maple Transactions, and I remind you that it has a "Diamond Class" designation, which means there are no page charges to authors, and the articles are free to read for everyone.  This means that there's some volunteer labour needed, of course: you have to write the articles, and what we want is that you write articles that people in the Maple community actually want to read.

I'd also like to thank the copyeditor, Michelle Hatzel, for her very hard work on this issue.  She's really made a difference, and I think you will be able to see it.   

In the book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the another remarkable theorem was proved: any two flat bounded regions can be cut by a single straight line so that each of these regions is divided into two regions of equal area (the second  pancake problem). This is an existence theorem which does not provide any way to find this cut. In this post I made an attempt to find such cut for any 2 convex regions on the plane bounded by a piecewise smooth self-non-intersecting curves.
The Each_Into_2_Equal_Areas procedure returns a list of coordinates of 4 endpoints of the cutting segments. This procedure significantly uses my old procedures  Area  and  Picture , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal arguments of the Each_Into_2_Equal_Areas procedure are the lists  L1  and  L2 specifying the boundaries of the regions to be cut. When specifying  L1  and  L2 , 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 regions and cutting segments. For ease of use, the code for the  Area  and  Picture   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:

Picture := proc(L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border;
uses plottools; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else 
P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) fi else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) fi; fi; od; 
Q := [seq(P[i], i = 1 .. nops(L))]; Border := curve([op(Q), Q[1]], op(Boundary)); [polygon(Q, C), Border] 
end proc:

Each_Into_2_Equal_Areas:=proc(L1::list, L2::list)
local D, n, m, L10, L20, S1,S2, f, L11, L21, i, j, k, s, A, B, C , sol;

f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L10:=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), L1);
L20:=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), L2);
S1:=Area(L1); S2:=Area(L2);  
n:=nops(L1); m:=nops(L2);

for i from 1 to n do
for j from i to n do

for k from 1 to m do
for s from k to m do

if not ((nops({i,j})=1 and type(L1[i],listlist)) or (nops({k,s})=1 and type(L2[k],listlist))) then

A:=eval(L10[i,1],t=t1): 
B:=eval(L10[j,1],t=t2):
C:=eval(L20[k,1],t=t3): 
D:=eval(L20[s,1],t=t4):

L11:=`if`(j=i,[subsop([2,2]=t1..t2,L10[i]),[B,A]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]], [subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),op(L10[i+1..j-1]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]])):

L21:=`if`(s=k,[subsop([2,2]=t3..t4,L20[k]),[D,C]],`if`(s=k+1,[subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]], [subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),op(L20[k+1..s-1]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]])):

sol:=fsolve(simplify({Area(L11)-S1/2,Area(L21)-S2/2,eval(f(A,B),[x=C[1],y=C[2]]),eval(f(A,B),[x=D[1],y=D[2]])}),{t1=op([2,2,1],L10[i])..op([2,2,2],L10[i]),t2=op([2,2,1],L10[j])..op([2,2,2],L10[j]),t3=op([2,2,1],L20[k])..op([2,2,2],L20[k]),t4=op([2,2,1],L20[s])..op([2,2,2],L20[s])});

if type(sol,set(`=`)) then  return eval([A,B,C,D],sol) fi;

fi;
od: od: od: od:
end proc:

Pic := proc(L1,L2,col1,col2,Size:=[800,400])
local P1, P2, P3, T, P;
uses plots, plottools;
P1, P2 := Picture(L1, color=col1), Picture(L2, color=col2):
P3 := line(Sol[1..2][],color=red,thickness=3), line(Sol[3..4][],color=red,thickness=3), line(Sol[1],Sol[4],linestyle=2,thickness=3,color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
P:=pointplot(Sol, symbol=solidcircle, color=red, symbolsize=10);
display(P1,P2,P3,T,P, scaling=constrained, size=Size, axes=none);
end proc: 

Examples of use.

local D:
L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+16,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+16,5*sin(t)/2],t=Pi..2*Pi],[[21,0],[16,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue",[900,400]);

The specified regions may overlap:

L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+9,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+9,5*sin(t)/2],t=Pi..2*Pi],[[14,0],[9,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2):  Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

If there is a solution for which the cutting segments intersect the boundary of each of the regions at 2 points, then the procedure also works for such non-convex regions:

L1:=[[[cos(t),sin(t)],t=Pi/3..2*Pi-Pi/3],[[cos(-t)+1,sin(-t)],t=-Pi-Pi/3..-Pi+Pi/3]]:
L2:=[[[cos(t)+2,sin(t)],t=-Pi/6..Pi+Pi/6],[[cos(-t)+2,sin(-t)-1],t=-5*Pi/6..-Pi/6]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

A number of other interesting examples can be found in the attached file.

Each_Into_2_Equal_Areas1.mw

After a long chat with ChatGPT I finally received a fully working code for a proc for a generalized Woodbury Identity for the inversion of the sum of two or more positive definite matrices.
I was inspired by a family member who is a trained professional programmer, who told me that in his professional work he uses ChatGTP for an initial draft of his program.  
I did find out that ChatGTP makes errors: from simple ones like writing 'Simplify' instead of 'simplify' to serious conceptual errors, for example in recursive loops. However, ChatGTP seems to 'understand' the error after given specific feedback. Although, this does not mean that the next proposal does not contain the same logical error. But after a long chat I received a nice proc that seems to work. 
My second surprise was that Gemini suggested a formula for the generalized Woodbury lemma that was unknown to me, and I was unable to find on Scholar Google or https://math.stackexchange.com. Based on a special case of that formula, I was able to write the second proc myself. 
In conclusion, to start working it can be helpful to collaborate with AI friend, a little patience may help, AI may not be astute as someone on Mapleprimes wrote, but neither am I. I am now retired, and it is fun to play with Maple and AI. 
By the way, the search term Woodbury did not give a single hit on Mapleprimes.With_a_little_help_from_my_friends.mw
kind regards,Harry 

 

Download blogpost-algebra.mw

Download lim-int-vv.mw

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.