Product Tips & Techniques

Tips and Tricks on how to get the most about Maple and MapleSim

Sometimes, it’s the little things. Those little improvements that make a good tool even better. Sometimes, it’s as simple as an easy shorthand notation that allows you to create and label points on a graph with a single command. Just to pick a totally random example.

 

A screenshot of a Maple Learn document containing a math cell and a plot. The math cell reads 'A(1,2)'. The plot show a point plotted at (1,2) with the label 'A'.

Okay, maybe it’s not totally random. Maybe this new point notation is one of our newest features in Maple Learn, and maybe it’s now easy and quick for you to create labeled points to your heart’s content. Maybe you could learn more about all the ins and outs of this new feature by checking out the how-to document.

But I can’t make any guarantees, of course.

That said, if this hypothetical scenario were true, you would also be able to see it in action in our new document on the proof of the triangle inequality.

A screenshot of a Maple Learn document. The left side shows an explanation of how the triangles are constructed for Euclid's proof of the triangle inequality. The right side shows an adjustable graph of said triangles.

With this document, you can explore a detailed (and interactive!) visualization of the proof using Euclidean geometry. You can adjust the triangles to see for yourself that the sum of the lengths of any two sides must be greater than the third side, read through the explanation to see the mathematical proof, and challenge yourself with the questions it leaves you to answer. And those points on those triangles? Labeled. Smoothly and easily. I wonder how they might have done it?

We hope you enjoy the new update! Let us know what other features you want to see in Maple Learn, and we’ll do our best to turn those dreams into reality.

We all know that math is beautiful in and of itself—but sometimes students might need a little convincing. What better way to do that then sprucing up your math with a little colour? With Maple Learn, plot colours are fully customizable. We have several colour palettes to choose from—want your document to evoke the delicate tones of springtime? Looking for a palette that’s colourblind friendly? Or maybe you’re just nostalgic for the colours of Maple V? All these options and more are available for making your graphs colourful and coordinated. But maybe you’re the kind of person who wants to go against the grain, and you laugh in the face of predetermined colour coordination. Don’t worry, we’ve got you covered too! With our colour selector, you can also choose your own custom colours. The full colour spectrum is right at your fingertips. To learn more about how to customize the colours on your document, check out this How-To guide.

And of course, the potential for colour inevitably leads to the potential for art. Our Maple Learn Art Gallery has plenty of fun and colourful works you can admire and contemplate (and maybe even draw inspiration from!). One of our most recent and most colourful additions is this document showcasing the history of the rainbow pride flag, in honour of June being Pride Month. You can use the slider to move through time, letting you see how the colours on the flag evolved and read about the meanings behind them. And, thanks to the colour selector, the colours match the precise shades used for the original flags! That’s the magic of hexadecimal colours for you.


Hold on—the magic of hexadecimal colours, I hear you ask? What an enticing concept. If only we had some kind of document, perhaps one made in Maple Learn, that explained how hexadecimal colours worked and included an interactive example so that you could easily see how the red, green, and blue colour values blend together to create any given colour… Too bad we don’t!

Just kidding. Of course we do.

If all these colours have inspired you, be sure to check out our Call for Creative Works for the upcoming Maple Conference! Maybe your colourful creation could be this year’s winner.

Ever wonder how to show progress updates from your executing code without printing new lines each time?

One way to do this is to use a TextArea component and the DocumentTools package. The TextArea could be inserted from the Components Palette in Maple, or programmatically like so:

restart;

with(DocumentTools):

with(DocumentTools:-Components):

with(DocumentTools:-Layout):

s := "0": #initial text value

T := TextArea(s, identity = "TextArea0"):
xml := Worksheet(Group(Input(Textfield(T)))):

insertedname:=InsertContent(xml)[1,1]: #find the inserted component name in case changed

for i to 10 do #start the demonstration procedure
   Threads:-Sleep(1);
   SetProperty(insertedname,value,sprintf("%d",i),refresh=true);
end do:

Maplets:-Examples:-Message("Done");


Download text-area-update-progress.mw

This post discusses a solution for modeling a traveling load on Maplesim's Flexible Beam component and provides an example of a bouncing load.

The idea for the above example came from an attempt to reproduce a model of a mass sliding on a beam from MapleSim's model gallery. However, reproducing it using contact components in combination with the Flexible Beam component turned out to be not straightforward, and this will be discussed in the following.

To simulate a traveling load on the Flexible Beam component, one could apply forces at discrete locations for a certain duration. However, the fidelity of this approach is limited by the number of discrete locations, which must be defined using the Flexible Beam Frame component, as well as the way in which the forces are activated.

One potential solution to address the issue of temporal activation of forces is to attach contact elements (such as Rectangle components) at distinct locations along the beam, which are defined by Flexible Beam Frame components, and make contact using a spherical or toroidal contact element. However, this approach also introduces two new problems:

  • An additional bending moment is generated when the load is not applied at the center of the contact element's attachment point, the Flexible Beam Frame component. Depending on the length of the contact element, deformations caused by this moment can be greater than the deformation caused by the force itself when the force is applied at the ends of the contact element. Overall, this unwanted moment makes the simulation unrealistic and must be avoided.
  • When the beam bends, a gap (see below) or an overlap is created between adjacent Rectangle components. If there is a gap, the object exerting a force on the beam can fall through it. Overlaps can create differences in dynamic behavior when the radius of curvature of the beam is on the opposite side of the point of contact.

To avoid these problems, the solution presented here uses an intermediate kinematic chain (encircled in yellow below) that redistributes the contact force on the Rectangle component on two support points (ports to attach Flexible Beam Frame components) in a linear fashion.

 

 

To address gaps, the contact element (Rectangle) attached to the kinematic chain has the same width as the chain and connects to the adjacent contact elements via multibody frames. In the image below, 10 contact elements are laid on top of a single Flexible Beam component, like a belt made out of tiles. The belt has to be pinned to the flexible beam at one location (highlighted in yellow). The location of this fixed point determines how the flexible beam is loaded by tangential contact forces (friction forces) and should be selected carefully.

 

 

Some observations on the attached model:

  • Low damping and high initial potential energy of the ball can result in a failed simulation (due to constraint projection failure). Increasing the number of elastic coordinates has a similar effect. Constraint projection can be turned off in the simulation settings to continue simulation.
  • The bouncing ball excites several eigenmodes at once, causing the beam to wiggle chaotically in combination with the varying bouncing frequency of the ball. A similar looking effect can also be achieved with special initial conditions, as demonstrated with Maple in this excellent post on Euler beams and partial differential equations.
  • Repeated simulations with low damping lead to different results (an indication of chaotic behavior; see three successive simulations below (gold) compared to a saved solution(red)). The moment in the animation when the ball travels backward represents a metastable equilibrium point of the simulation. This makes predictions beyond this point difficult, as the behavior of the system is highly dependent on the model parameters. Whether the reversal is a simulation artifact or can happen in reality remains to be seen. Overall, this example could evolve into a nice experimental fun project for students.

  • Setting the gravitational constant for Mars, everything is different. I could not reproduce the fun factor on Earth. A reason more to stay ,-)

Ball_bouncing_on_a_flexible_beam.msim

 

 

# countourplot3d piggybacks on top of plot3d.
# For the "coloring=[lowColor, highColor]", the "filledregions=true" option must be present.
# If "filledregions=true" is not present, plot3d will throw an error.
# This code shows the three cases, only one of which will work.
Note to support. I cannot add a new tag. contourplot3d should be a tag.

restart;
with(plots);
with(ColorTools);
cGr4s := Color([0.50, 0.50, 0.50]);
contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, color = black, thickness = 3, coloring = [cGr4s, cGr4s], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);
contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, filledregions = false, color = black, thickness = 3, coloring = [cGr4s, cGr4s], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);
contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, filledregions = true, color = black, thickness = 3, coloring = [cGr4s, cGr4s], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);

 

The recent Maple 2023 release comes with a multitude of new features, including a new Canvas Scripting Gallery full of templates for creating interactive Maple Learn documents.

The Maple Learn Scripting Gallery can be accessed through Maple, by searching “BuildInteractiveContent Maple2023” in the search bar at the top of the application and clicking on the only result that appears. This will bring you to the help page titled “Build and Share Interactive Content”, which can also be found by searching “scripting gallery” in the search bar of a Maple help page window. The link to the Maple Learn Scripting Gallery is found under the “Canvas Scripting” section on this help page and clicking on it will open a Maple workbook full of examples and templates for you to explore.

The interactive content in the Scripting Gallery is organized into five main categories – Graphing, Visualization, Quiz, Add-ons and Options, and Applications Optimized for Maple Learn – each with its own sub-categories, templates, and examples.

One of the example scripts that I find particularly interesting is the “Normal Distribution” script, under the Visualizations category.

 

 

All of the code for each of the examples and templates in the gallery is provided, so we can see exactly how the Normal Distribution script creates a Maple Learn canvas. It displays a list of grades, a plot for the grade distribution to later appear on, math groups for the data’s mean and variance, and finally a “Calculate” button that runs a function called UpdateStats.

The initial grades loaded into the document result in the below plot, created using Maple’s DensityPlot and Histogram functions, from the Statistics package. 




 

The UpdateStats function takes the data provided in the list of grades and uses a helper function, getDist, to generate the new plot to display the data, the distribution, the mean, and the variance. Then, the function uses a Script object to update the Maple Learn canvas with the new plot and information.

The rest of the code is contained in the getDist function, which uses a variety of functions from Maple’s Statistics package. The Normal Distribution script takes advantage of Maple’s ability to easily calculate mean and variance for data sets, and to use that information to create different types of random variable distributions.

Using the “Interactive Visualization” template, provided in the gallery, many more interactive documents can be created, like this Polyhedra Visualization and this Damped Harmonic Oscillator – both from the Scripted Gallery or like my own Linear Regression: Method of Least Squares document.


 

Another new feature of Maple 2023 is the Quiz Builder, also featured in the Scripting Gallery. Quizzes created using Quiz Builder can be displayed in Maple or launched as Maple Learn quizzes, and the process for creating such a quiz is short.

The QuizBuilder template also provides access to many structured examples, available from a dropdown list:


As an example, check out this Maple Learn quiz on Expected Value: Continuous Practice. Here is what the quiz looks like when generated in Maple:


 

This quiz, in particular, is “Fill-in the blank” style, but Maple users can also choose “Multiple Choice”, “True/False”, “Multiple Select”, or “Multi-Line Feedback”. It also makes use of all of the featured code regions from the template, providing functionality for checking inputted answers, generating more questions, showing comprehensive solutions, and providing a hint at the press of a button.

Check out the Maple Learn Scripting Gallery for yourself and see what kinds of interactive content you can make for Maple and Maple Learn!

 

Physics Courseware Support: Mechanics

Hi
The attached worksheet is the final version that appears in Maple 2023 as Courseware support for Mechanics in the context of Physics courses. Everything below also works in Maple 2022.2 with the last Maplesoft Physics Updates for that release..

What follows is presented as "Topic > Problem > Solution", with typical symbolic problems and how you can solve them on a worksheet. As such, this material does not intend to compete with textbooks nor with teacher's notes but to be a helpful complement, as in "what can computer algebra really do to support the learning activity". Mainly, allow for focusing the logic and thinking while the computer takes care of the intricacies of the algebraic manipulations, that when computing with paper and pencil so frequently take mostly all of our focus.

The material, thus, has 70 solved problems covering all the sort-of-syllabus of hyperlinks below. The presentation uses notation as in textbooks and illustrates different techniques, several not present in help pages. It also shows why it is relevant to have a Vectors package that handles abstract vectors as well as projections using unit vectors, not matrix representations for them. Your feedback about everything you see in the worksheet - suggestions for new topics or problems, or anything else - can be useful and is welcome.

Due to the length of this material (~100 pages), out of the 70 problems, below I left open (visible) the Solution sections of only a few of them, illustrating different things, also new functionality e.g. the first and last ones. That is sufficient to have an idea of what this is about. At the end there is a Maple worksheet with the same contents and a PDF file of the same with all the sections open.



With the best wishes for 2023.

 

  

Explore. While learning, having success is a secondary goal: using your curiosity as a compass is what matters. Things can be done in many different ways, take full permission to make mistakes. Computer algebra can transform the algebraic computation part of physics into interesting discoveries and fun.

  

 

  

The following material assumes knowledge of how to use Maple. If you feel that is not your case, for a compact introduction on reproducing in Maple the computations you do with paper and pencil, see sections 1 to 5 of the Mini-Course: Computer Algebra for Physicists . Also, the presentation assumes an understanding of the subjects and the style is not that of a textbook. Instead, it focuses on conveniently using computer algebra to support the practice and learning process. The selection of topics follows references [1] and [2] at the end. Maple 2023.0 includes Part I. Part II is forthcoming.

NULL

Part I

1. 

Position, velocity and acceleration in Cartesian, cylindrical and spherical coordinates

a. 

The position `#mover(mi("r"),mo("→"))`(t)as a function of time

b. 

The velocity `#mover(mi("v"),mo("→"))`(t)

c. 

The acceleration   `#mover(mi("a"),mo("→"))`(t)

d. 

Deriving these formulas

e. 

Velocity and acceleration in the case of 2-dimensional motion on the x, y plane

1. 

The equations of motion

a. 

A single particle

i. 

The equations of motion - vectorial form

ii. 

The case of constant acceleration

iii. 

Motion under gravitational force close to the Earth's surface

iv. 

Motion under gravitational force not close to the Earth's surface

A. 

Circular motion

B. 

Escape velocity

i. 

Different acceleration in different regions

ii. 

The equations of motion using tensor notation

A. 

Cartesian coordinates

B. 

Curvilinear coordinates

a. 

Many-particle systems

i. 

Center of mass

ii. 

The equations of motion

iii. 

Static: reactions of planes and tensions on cables

a. 

Lagrange equations

i. 

Motion of a pendulum

1. 

Conservation laws

a. 

Work

b. 

Conservation of the total energy of a closed system or a system in a constant external field

c. 

Conservation of the total momentum of a closed system

d. 

Conservation of angular momentum

e. 

Cyclic coordinates

1. 

Integration of the equations of motion

a. 

Motion in one dimension

b. 

Reduced mass

i. 

The two-body problem

ii. 

A many-body problem

a. 

Motion in a central field

b. 

Kepler's problem

1. 

Small Oscillations

a. 

Free oscillations in one dimension

b. 

Forced oscillations

c. 

Oscillations of systems with many degrees of freedom

1. 

Rigid-body motion

a. 

Angular velocity

b. 

Inertia tensor

c. 

Angular momentum of a rigid body

d. 

The equations of motion of a rigid body

1. 

Non-inertial coordinate systems

a. 

Coriolis force and centripetal force

 

Part II (forthcoming)

1. 

The Hamiltonian and equations of motion; Poisson brackets

2. 

Canonical transformations

3. 

The Hamilton-Jacobi equation

 

Position, velocity and acceleration in Cartesian, cylindrical and spherical coordinates

 

 

Load the Physics:-Vectors  package

 

with(Physics:-Vectors)

[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface, ParametrizeVolume, Setup, Simplify, `^`, diff, int]

(1)

 

Depending on the geometry of a problem, it can be convenient to work with either Cartesian or curvilinear coordinates. In an arbitrary reference system, the position in Cartesian coordinates and the basis of unitary vectors`#mover(mi("i"),mo("∧"))`, `#mover(mi("j"),mo("∧"))`, `#mover(mi("k"),mo("∧"))`is given by

r_ = _i*x+_j*y+_k*z

r_ = _i*x+_j*y+_k*z

(2)

NULL

Problem

Rewrite the position vector `#mover(mi("r"),mo("→"))` in cylindrical and spherical coordinates

Solution

   

 

Starting from the position in the Cartesian system, now as functions of the time to allow for differentiation, first note that the Cartesian unit vectors `#mover(mi("i"),mo("∧"))`, `#mover(mi("j"),mo("∧"))`, `#mover(mi("k"),mo("∧"))` do not depend on time, they are constant vectors. So `#mover(mi("r"),mo("→"))`(t) is entered as

 

restart; with(Physics:-Vectors)

 

r_(t) = x(t)*_i+y(t)*_j+z(t)*_k

r_(t) = x(t)*_i+y(t)*_j+z(t)*_k

(20)

Before proceeding further, use a compact display to more clearly visualize the following expressions. When in doubt about the contents behind a given display, input show as shown below.

CompactDisplay((x, y, z, rho, r, theta, phi, _rho, _r, _theta, _phi)(t))

_phi(t)*`will now be displayed as`*_phi

(21)

 

For the velocity and acceleration, note the dot notation for derivatives with respect to t

v_(t) = diff(rhs(r_(t) = x(t)*_i+y(t)*_j+z(t)*_k), t)

v_(t) = (diff(x(t), t))*_i+(diff(y(t), t))*_j+(diff(z(t), t))*_k

(22)

 

show

v_(t) = (diff(x(t), t))*_i+(diff(y(t), t))*_j+(diff(z(t), t))*_k

(23)

a_(t) = diff(rhs(v_(t) = (diff(x(t), t))*_i+(diff(y(t), t))*_j+(diff(z(t), t))*_k), t)

a_(t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k

(24)

NULL

The position `#mover(mi("r"),mo("→"))`(t)as a function of time

 

 

Problem

Given the position vector as a function of the time t, rewrite it in cylindrical and spherical coordinates while making the curvilinear unit vectors' time dependency explicit.

Solution

   

 

The velocity `#mover(mi("v"),mo("→"))`(t)

 

 

Problem

Rewrite the velocity `#mover(mi("v"),mo("→"))`(t) = diff(`#mover(mi("r"),mo("→"))`(t), t) in cylindrical and spherical coordinates while making the curvilinear unit vectors' time dependency explicit .

Solution

   

 

The acceleration `#mover(mi("a"),mo("→"))`(t)

 

 

Problem

Rewrite the acceleration `#mover(mi("a"),mo("→"))`(t) = diff(`#mover(mi("r"),mo("→"))`(t), t, t)in cylindrical and spherical components while making the curvilinear unit vectors' time dependency explicit.

Solution

   

 

Deriving these formulas

 

 

All these results for the position `#mover(mi("r"),mo("→"))`, velocity `#mover(mi("v"),mo("→"))` and acceleration `#mover(mi("a"),mo("→"))` are based on the differentiation rules for cylindrical and spherical unit vectors. It is thus instructive to also be able to derive any of these formulas; for that, we need the differentiation rule for the unit vectors. For example, for the spherical ones

 

restart; with(Physics:-Vectors); CompactDisplay((x, y, z, rho, r, theta, phi, _rho, _r, _theta, _phi)(t), quiet)

 

map(%diff = diff, ([_r, _theta, _phi])(t), t)

[%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*_rho(t)]

(38)

The above result contains, in the last equation, the cylindrical radial unit vector `#mover(mi("ρ",fontstyle = "normal"),mo("∧"))`(t); rewrite it in the spherical basis

_rho(t) = ChangeBasis(_rho(t), spherical)

_rho(t) = sin(theta(t))*_r(t)+cos(theta(t))*_theta(t)

(39)

So the differentiation rules for spherical unit vectors, with the result expressed in the spherical system, are

subs(_rho(t) = sin(theta(t))*_r(t)+cos(theta(t))*_theta(t), [%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*_rho(t)])

[%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))]

(40)

NULL

Problem

With this information at hand, derive, in steps, the expressions for the velocity and acceleration in cylindrical and spherical coordinates

Solution

 

 

We want to compute

%diff(r_(t) = r(t)*_r(t), t)

%diff(r_(t) = r(t)*_r(t), t)

(41)

expand(%diff(r_(t) = r(t)*_r(t), t))

%diff(r_(t), t) = %diff(r(t), t)*_r(t)+r(t)*%diff(_r(t), t)

(42)

Introducing the differentiation rules (40) for the unit vectors

subs([%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))], %diff(r_(t), t) = %diff(r(t), t)*_r(t)+r(t)*%diff(_r(t), t))

%diff(r_(t), t) = %diff(r(t), t)*_r(t)+r(t)*((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t))

(43)

Performing the inert (grayed) derivatives

value(%diff(r_(t), t) = %diff(r(t), t)*_r(t)+r(t)*((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t)))

diff(r_(t), t) = (diff(r(t), t))*_r(t)+r(t)*((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t))

(44)

In the same way, for the acceleration

%diff(r_(t) = r(t)*_r(t), t, t)

%diff(r_(t) = r(t)*_r(t), t, t)

(45)

expand(%diff(r_(t) = r(t)*_r(t), t, t))

%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*%diff(_r(t), t)+r(t)*%diff(%diff(_r(t), t), t)

(46)

subs([%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))], %diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*%diff(_r(t), t)+r(t)*%diff(%diff(_r(t), t), t))

%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t))+r(t)*%diff((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), t)

(47)

expand(%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t))+r(t)*%diff((diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), t))

%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*(diff(theta(t), t))*_theta(t)+2*%diff(r(t), t)*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*%diff(%diff(theta(t), t), t)*_theta(t)+r(t)*(diff(theta(t), t))*%diff(_theta(t), t)+r(t)*%diff(%diff(phi(t), t), t)*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*%diff(theta(t), t)*cos(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*sin(theta(t))*%diff(_phi(t), t)

(48)

subs([%diff(_r(t), t) = (diff(theta(t), t))*_theta(t)+(diff(phi(t), t))*sin(theta(t))*_phi(t), %diff(_theta(t), t) = -(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t), %diff(_phi(t), t) = -(diff(phi(t), t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))], %diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*(diff(theta(t), t))*_theta(t)+2*%diff(r(t), t)*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*%diff(%diff(theta(t), t), t)*_theta(t)+r(t)*(diff(theta(t), t))*%diff(_theta(t), t)+r(t)*%diff(%diff(phi(t), t), t)*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*%diff(theta(t), t)*cos(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*sin(theta(t))*%diff(_phi(t), t))

%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*(diff(theta(t), t))*_theta(t)+2*%diff(r(t), t)*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*%diff(%diff(theta(t), t), t)*_theta(t)+r(t)*(diff(theta(t), t))*(-(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t))+r(t)*%diff(%diff(phi(t), t), t)*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*%diff(theta(t), t)*cos(theta(t))*_phi(t)-r(t)*(diff(phi(t), t))^2*sin(theta(t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))

(49)

value(%diff(%diff(r_(t), t), t) = %diff(%diff(r(t), t), t)*_r(t)+2*%diff(r(t), t)*(diff(theta(t), t))*_theta(t)+2*%diff(r(t), t)*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*%diff(%diff(theta(t), t), t)*_theta(t)+r(t)*(diff(theta(t), t))*(-(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t))+r(t)*%diff(%diff(phi(t), t), t)*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*%diff(theta(t), t)*cos(theta(t))*_phi(t)-r(t)*(diff(phi(t), t))^2*sin(theta(t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t)))

diff(diff(r_(t), t), t) = (diff(diff(r(t), t), t))*_r(t)+2*(diff(r(t), t))*(diff(theta(t), t))*_theta(t)+2*(diff(r(t), t))*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*(diff(diff(theta(t), t), t))*_theta(t)+r(t)*(diff(theta(t), t))*(-(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t))+r(t)*(diff(diff(phi(t), t), t))*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*(diff(theta(t), t))*cos(theta(t))*_phi(t)-r(t)*(diff(phi(t), t))^2*sin(theta(t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t))

(50)

Collect vector components

Physics:-Vectors:-Collect(diff(diff(r_(t), t), t) = (diff(diff(r(t), t), t))*_r(t)+2*(diff(r(t), t))*(diff(theta(t), t))*_theta(t)+2*(diff(r(t), t))*(diff(phi(t), t))*sin(theta(t))*_phi(t)+r(t)*(diff(diff(theta(t), t), t))*_theta(t)+r(t)*(diff(theta(t), t))*(-(diff(theta(t), t))*_r(t)+(diff(phi(t), t))*cos(theta(t))*_phi(t))+r(t)*(diff(diff(phi(t), t), t))*sin(theta(t))*_phi(t)+r(t)*(diff(phi(t), t))*(diff(theta(t), t))*cos(theta(t))*_phi(t)-r(t)*(diff(phi(t), t))^2*sin(theta(t))*(sin(theta(t))*_r(t)+cos(theta(t))*_theta(t)))

diff(diff(r_(t), t), t) = (2*r(t)*(diff(theta(t), t))*(diff(phi(t), t))*cos(theta(t))+2*(diff(r(t), t))*(diff(phi(t), t))*sin(theta(t))+r(t)*(diff(diff(phi(t), t), t))*sin(theta(t)))*_phi(t)+(-r(t)*(diff(phi(t), t))^2*sin(theta(t))^2-r(t)*(diff(theta(t), t))^2+diff(diff(r(t), t), t))*_r(t)+(-(diff(phi(t), t))^2*sin(theta(t))*r(t)*cos(theta(t))+2*(diff(r(t), t))*(diff(theta(t), t))+r(t)*(diff(diff(theta(t), t), t)))*_theta(t)

(51)

NULL

 

Summary

• 

You can express `#mover(mi("r"),mo("→"))`(t), `#mover(mi("v"),mo("→"))`(t) and `#mover(mi("a"),mo("→"))`(t)in any of the Cartesian, cylindrical or spherical systems via three different methods: 1) using the ChangeBasis command 2) differentiating 3) deriving the formulas by differentiating in steps, starting from the differentiation rules for the curvilinear unit vectors.

Velocity and acceleration in the case of 2-dimensional motion on the x, y plane

 


Problem

Derive formulas for velocity and acceleration in the case of 2-dimensional motion on the x, y plane, starting from the general 3-dimensional formulas above, e.g. (44) and (51) in spherical coordinates. Specialize the resulting formulas for the case of circular motion.

Solution

   

 

The equations of motion

 

 

A single particle

 

 

restart; with(Physics:-Vectors); CompactDisplay((r_, p_, F_, L_, N_)(t))

N_(t)*`will now be displayed as`*N_

(62)

The equation of motion of a single particle is Newton's 2^nd law

F_(t) = m*(diff(r_(t), t, t))

F_(t) = m*(diff(diff(r_(t), t), t))

(63)

where diff(`#mover(mi("r"),mo("→"))`(t), t, t) = `#mover(mi("a"),mo("→"))`(t) is the acceleration and m*(diff(`#mover(mi("r"),mo("→"))`(t), t)) = `#mover(mi("p"),mo("→"))`(t) is the linear momentum, so in terms of  `#mover(mi("p"),mo("→"))`

F_(t) = diff(p_(t), t)

F_(t) = diff(p_(t), t)

(64)

We define the angular momentum `#mover(mi("L"),mo("→"))` of a particle, and the torque `#mover(mi("N"),mo("→"))` acting upon it, as

L_(t) = `&x`(r_(t), p_(t))

L_(t) = Physics:-Vectors:-`&x`(r_(t), p_(t))

(65)

N_(t) = `&x`(r_(t), F_(t))

N_(t) = Physics:-Vectors:-`&x`(r_(t), F_(t))

(66)

Differentiating the definition of `#mover(mi("L"),mo("→"))`

diff(L_(t) = Physics:-Vectors:-`&x`(r_(t), p_(t)), t)

diff(L_(t), t) = Physics:-Vectors:-`&x`(diff(r_(t), t), p_(t))+Physics:-Vectors:-`&x`(r_(t), diff(p_(t), t))

(67)

Since diff(`#mover(mi("r"),mo("→"))`(t), t) = `#mover(mi("v"),mo("→"))` is parallel to `#mover(mi("p"),mo("→"))` = m*`#mover(mi("v"),mo("→"))`, the first term in the above cancels, and in the second term, from (64), diff(`#mover(mi("p"),mo("→"))`(t), t) = `#mover(mi("F"),mo("→"))`

eval(diff(L_(t), t) = Physics:-Vectors:-`&x`(diff(r_(t), t), p_(t))+Physics:-Vectors:-`&x`(r_(t), diff(p_(t), t)), [diff(r_(t), t) = 0, diff(p_(t), t) = F_(t)])

diff(L_(t), t) = Physics:-Vectors:-`&x`(r_(t), F_(t))

(68)

from which

subs((rhs = lhs)(diff(L_(t), t) = Physics:-Vectors:-`&x`(r_(t), F_(t))), N_(t) = Physics:-Vectors:-`&x`(r_(t), F_(t)))

N_(t) = diff(L_(t), t)

(69)

NULL

• 

As discussed below , in the case of a closed system, `#mover(mi("F"),mo("→"))` = 0 and these two equations result in

diff(`#mover(mi("p"),mo("→"))`(t), t) = 0, diff(`#mover(mi("L"),mo("→"))`(t), t) = 0

that is, the linear and angular momentum are conserved quantities. Note that diff(`#mover(mi("L"),mo("→"))`(t), t) = 0 does not require that `#mover(mi("F"),mo("→"))` = 0, only that `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("F"),mo("→"))`) = 0.

 

The equations of motion - vectorial form

 

 

Problem

Assuming that the acceleration is known as a function of t, compute:

a) The trajectory `#mover(mi("r"),mo("→"))`(t)starting from `#mover(mi("a"),mo("→"))`(t) = diff(`#mover(mi("r"),mo("→"))`(t), t, t)
b) A solution for each of the three Cartesian components

c) A solution for generic initial conditions

Solution

 

 

restart; with(Physics:-Vectors); CompactDisplay((x, y, z, rho, r, theta, phi, _rho, _r, _theta, _phi)(t), quiet)
 

 

a) Let `#mover(mi("r"),mo("→"))`(t) be the position of the particle in a reference system; then, the velocity and acceleration are given by

v_(t) = diff(r_(t), t)

v_(t) = diff(r_(t), t)

(70)

a_(t) = diff(r_(t), t, t)

a_(t) = diff(diff(r_(t), t), t)

(71)

If the acceleration is known as a function of t, the trajectory is computed by integrating (71)

dsolve(a_(t) = diff(diff(r_(t), t), t))

r_(t) = Int(Int(a_(t), t), t)+c__1_*t+c__2_

(72)

where the vectorial integration constants, c__1_ and c__2_, are specified by the initial conditions of the problem (see c) below), typically by the position and velocity at some instant, say %eval(r_(t), `=`(t, t__0)) = r__0_ and %eval(diff(r_(t), t), `=`(t, t__0)) = v__0_.

_______________________________________

 

b) The integration of vectorial equations is also frequently performed after expressing `#mover(mi("r"),mo("→"))`(t), `#mover(mi("v"),mo("→"))`(t) and `#mover(mi("a"),mo("→"))`(t) in a particular system of coordinates. For example, in the Cartesian system (71) has the form

(rhs = lhs)(a_(t) = (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k)

(diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a_(t)

(73)

Now suppose that the three components of the acceleration are known as a function of time

subs(a_(t) = a__x(t)*_i+a__y(t)*_j+a__z(t)*_k, (diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a_(t))

(diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a__x(t)*_i+a__y(t)*_j+a__z(t)*_k

(74)

Vectorial equations like this one can be integrated directly, provided that they are expressed in a particular system of coordinates and the unit vectors are constant or known expressions of the time

dsolve((diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a__x(t)*_i+a__y(t)*_j+a__z(t)*_k)

_i*x(t)+_j*y(t)+_k*z(t) = _i*(Int(Int(a__x(t), t), t)+c__3*t+c__4)+_j*(Int(Int(a__y(t), t), t)+c__5*t+c__6)+_k*(Int(Int(a__z(t), t), t)+c__7*t+c__8)

(75)

_______________________________________

 

c) The vectorial initial conditions r__0_ and v__0_, specifying the integration constants {c__3, c__4, c__5, c__6, c__7, c__8}, can also be written in components

x(t__0) = x__0, y(t__0) = y__0, z(t__0) = z__0(t), diff(x(t__0), t__0) = v__x0, diff(y(t__0), t__0) = v__y0, diff(z(t__0), t__0) = v__z0

x(t__0) = x__0, y(t__0) = y__0, z(t__0) = z__0(t), diff(x(t__0), t__0) = v__x0, diff(y(t__0), t__0) = v__y0, diff(z(t__0), t__0) = v__z0

(76)

Passing this information, the system can be solved taking these initial conditions into account -

dsolve([(diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a__x(t)*_i+a__y(t)*_j+a__z(t)*_k, x(t__0) = x__0, y(t__0) = y__0, z(t__0) = z__0(t), diff(x(t__0), t__0) = v__x0, diff(y(t__0), t__0) = v__y0, diff(z(t__0), t__0) = v__z0], ({x, y, z})(t))

{x(t) = Int(Int(a__x(tau), tau = t__0 .. tau), tau = t__0 .. t)+v__x0*t-t__0*v__x0+x__0, y(t) = Int(Int(a__y(tau), tau = t__0 .. tau), tau = t__0 .. t)+v__y0*t-t__0*v__y0+y__0, z(t) = Int(Int(a__z(tau), tau = t__0 .. tau), tau = t__0 .. t)+v__z0*t+c__4}

(77)

Note that a vectorial equation is also always equivalent to a system of equations, one for each of the components, with or without initial conditions:

convert((diff(diff(x(t), t), t))*_i+(diff(diff(y(t), t), t))*_j+(diff(diff(z(t), t), t))*_k = a__x(t)*_i+a__y(t)*_j+a__z(t)*_k, setofequations)

{diff(diff(x(t), t), t) = a__x(t), diff(diff(y(t), t), t) = a__y(t), diff(diff(z(t), t), t) = a__z(t)}

(78)

dsolve({diff(diff(x(t), t), t) = a__x(t), diff(diff(y(t), t), t) = a__y(t), diff(diff(z(t), t), t) = a__z(t)}, {x, y, z})

{x(t) = Int(Int(a__x(t), t), t)+c__7*t+c__8, y(t) = Int(Int(a__y(t), t), t)+c__5*t+c__6, z(t) = Int(Int(a__z(t), t), t)+c__3*t+c__4}

(79)

 

NULL

 

The case of constant acceleration

 

 

Problem

Starting from the vectorial equation (72) for `#mover(mi("r"),mo("→"))`(t), derive the formula for constant acceleration

Solution

   

 

Motion under gravitational force close to the Earth's surface

 

 

Problem

Derive a formula for motion under gravitational force close to the Earth's surface

Solution

   

 

Motion under gravitational force not close to the Earth's surface

 

 

The problem of two particles of masses m__1 and m__2 gravitationally attracted to each other, discarding relativistic effects, is formulated by Newton's law of gravity: the particles attract each other - so both move - with a force along the line that joins the particles and whose magnitude is proportional to 1/r^2, where r represents the distance between the particles (this problem is treated in general form  in the more advanced sections).

 

Problem

As a specific case, consider the problem of a particle of mass `≪`(m, M), where M is earth's mass, moving not close to the surface (if compared with the radius of earth).

Solution

   

 

Circular motion

 

 

Problem

Determine the angular velocity diff(phi(t), t) in the case of circular motion and show it is constant.

Solution

   

 

Escape velocity

 

 

Problem

Determine the velocity that a particle of mass m should have at Earth's surface in order to escape the planet's gravitational attraction.

Solution

   

 

Different acceleration in different regions

 

 

Problem

Suppose a particle is moving along the x axis according to

x(t) = t^3-8*t^2+18*t+3

a) Determine the regions where the motion has positive and negative acceleration. Compute the position at proc (t) options operator, arrow; infinity end proc.

b) Compute the velocity v__x(t)corresponding to x(t) = t^3-8*t^2+18*t+3, starting - not from this expression for x(t) but from the acceleration a__x(t) = diff(x(t), t, t)

Solution

   

 

The equations of motion using tensor notation

 

 

Using vector notation to formulate the equations of motion of a particle in Cartesian coordinates is relatively simple. However, for certain problems it may be advantageous to use curvilinear coordinates and / or tensor notation.

 

restart; with(Physics); with(Vectors)

NULL

Cartesian coordinates

   

Curvilinear coordinates

   

 

Many-particle systems

 

 

Center of Mass

 

Given a system of n particles of masses m__i with positions `#mover(mi("r"),mo("→"))`[i] in some frame of reference K, the center of mass of the system is defined as

 `#mover(mi("R"),mo("→"))` = (Sum(m[i]*`#mover(mi("r"),mo("→"))`[i], i = 1 .. n))/(Sum(m[i], i = 1 .. n))

The velocity of the center of mass is thus

  `#mover(mi("V"),mo("→"))` = diff(`#mover(mi("R"),mo("→"))`(t), t) and diff(`#mover(mi("R"),mo("→"))`(t), t) = (Sum(m[i]*(diff(`#mover(mi("r"),mo("→"))`[i](t), t)), i = 1 .. n))/(Sum(m[i], i = 1 .. n)) 

Problem

Consider a system of particles viewed from two systems of reference, K and K', that move with respect to each other at a constant velocity `#mover(mi("V"),mo("→"))` measured in K. Show that:

a) When `#mover(mi("V"),mo("→"))` is the velocity of the center of mass, the total momentum "(P )'" measured in K'  is equal to 0.

b) The relation between `#mover(mi("P"),mo("→"))` and the velocity `#mover(mi("V"),mo("→"))` of the center of mass, both measured in K, is the same as the relation `#mover(mi("p"),mo("→"))` = m*`#mover(mi("v"),mo("→"))` between the momentum, velocity and mass of a single particle of mass mu = Sum(m[i], i = 1 .. n).

Solution

   

 

The equations of motion

 

 

Problem
Show that, for a system of particles with total mass mu = Sum(m[i], i = 1 .. n), Newton's 2nd law for each particle `#mover(mi("F"),mo("→"))`[i] = m[i]*(diff(`#mover(mi("r"),mo("→"))`[i](t), t, t)) implies that `#mscripts(mi("F"),mi("ext"),none(),none(),mo("→"),none(),none())` = mu*(diff(`#mover(mi("R"),mo("→"))`(t), t, t)), where `#mover(mi("R"),mo("→"))` is the center of mass and `#mscripts(mi("F"),mi("ext"),none(),none(),mo("→"),none(),none())` is the external force applied to the system (it excludes the force that the particles exercise on each other).

Solution

   

 

Problem

Show that :

a) The total linear momentum `#mover(mi("P"),mo("→"))` satisfies diff(`#mover(mi("P"),mo("→"))`(t), t) = `#mscripts(mi("F"),mi("ext"),none(),none(),mo("→"),none(),none())`

b) The total torque `#mover(mi("N"),mo("→"))` = diff(`#mover(mi("L"),mo("→"))`(t), t) satisfies `#mover(mi("N"),mo("→"))` = Sum(`&x`(`#mover(mi("r"),mo("→"))`[i], `#mover(mi("f"),mo("→"))`[i, ext]), i = 1 .. n)

Solution

   

 

Static: reactions of planes and tensions on cables

 

 

Problem

A bar AB of weight w and length L has one extreme on a horizontal plane and the other on a vertical place, and is kept in that position by two cables AD and BC. The bar forms an angle alpha with the horizontal plane and its projection BC over this plane forms an angle beta with the vertical plane. The cable BC is on the same vertical plane as the bar.

 

Determine the reactions of the planes at A and B as well as the tensions on the cables.

Solution

   

 

Lagrange equations

 

 

restart; with(Physics:-Vectors); CompactDisplay((r_, v_)(t))

v_(t)*`will now be displayed as`*v_

(208)

 

In the case of a closed system, or a system in a constant external field, the equations of motions can also be derived from the knowledge of the kinetic energy T and the potential energy U . For this purpose, construct the Lagrange function L = T-U and derive the equations of motion as the Lagrange equations for L.

 

For closed systems, the potential energy U(`#mover(mi("r"),mo("→"))`[i])is related to the force acting on each particle by the equation `#mover(mi("F"),mo("→"))`[i] = -`∇__i`(U(`$`(r_[i], `=`(i, 1 .. n)))). Formally, "`∇__i`≡diff(L,r_[i]" is the gradient operator in the basis onto which `#mover(mi("r"),mo("→"))`[i] is projected, and with respect to its coordinates - say in Cartesian basis and coordinates "(∇)[i]=(i)*(`%diff`(,)+j*(`%diff`(,y))+k*(`%diff`(,z))".


The kinetic energy - say T - of a single particle is given by

T := (1/2)*m*v_(t)^2

(1/2)*m*Physics:-`^`(v_(t), 2)

(209)

Since the kinetic energy T is additive, a system of n particles has

T := %sum((1/2)*m*v_[i](t)^2, i = 1 .. n)

%sum((1/2)*m*Physics:-`^`(v_[i](t), 2), i = 1 .. n)

(210)

where `#mover(mi("v"),mo("→"))`[i] is the velocity of the ith particle. Generally speaking, the potential energy U(`$`(r_[i], `=`(i, 1 .. n))) of the system is a function of the coordinates `#mover(mi("r"),mo("→"))`[i] of the n particles, and the Lagrangian is defined asNULL

L = T-U(`$`(r_[i], `=`(i, 1 .. n)))

L = %sum((1/2)*m*Physics:-`^`(v_[i](t), 2), i = 1 .. n)-U(`$`(r_[i], i = 1 .. n))

(211)

The potential energy U(`#mover(mi("r"),mo("→"))`[i])is related to the force acting on each particle by the equation `#mover(mi("F"),mo("→"))`[i] = -`∇__i`(U(`$`(r_[i], `=`(i, 1 .. n)))). Formally, "`∇__i`≡diff(L,r_[i]" is the gradient operator in the basis onto which `#mover(mi("r"),mo("→"))`[i] is projected. Knowing the Lagrangian, we can derive the (Lagrange) equations of motion  as

%diff(%diff(L, v_[i]), t) = %diff(L, r_[i])

%diff(%diff(L, v_[i]), t) = %diff(L, r_[i])

(212)

NULL

Motion of a pendulum

 

 

Problem

Determine the Lagrangian and equation of motion of a plane pendulum with a mass m at its extremity and a suspension point which:

a) Moves uniformly over a vertical circumference with a constant frequency "omega."

 

b) Oscillates horizontally on the plane of the pendulum according to x = cos(omega*t).

c) Is fixed (omega = 0). Integrate the equation of motion for small oscillations.

Solution

   

 

Conservation laws

 

Work

 

By definition, the work performed by a force `#mover(mi("F"),mo("→"))` to move a particle an infinitesimal amount d*`#mover(mi("r"),mo("→"))` is Typesetting[delayDotProduct](`#mover(mi("F"),mo("→"))`.d, `#mover(mi("r"),mo("→"))`, true). Thus, the work to move it from `#mover(mi("A"),mo("→"))` to `#mover(mi("B"),mo("→"))` along some path C is

 

"((∫)[A]^(B) (F ) . ⅆr)[path=C]"

Problem

A particle is submitted to a force whose Cartesian components are given by" F[x]=a x^3+b x y^2+c z,  ""F[y]=a y^3+b x y^2,  "F[z] = c*z. Calculate the work done by this force when moving the particle along a straight line from the origin to a point (x[0], y[0], z[0]).

Solution

   

 

Conservation of the total energy of a closed system or system in a constant external field

 

 

Problem

Consider a closed system, or a system in a constant external field, for which the total force acting on the ith particle of the system can be derived from a potential, `#mscripts(mi("F"),mi("i"),none(),none(),mo("→"),none(),none())` = -`∇i`*U. Show that the total energy of the system is conserved.

Solution

   

 

Problem

 Consider a system of n particles in two reference systems K and K'  that move relative to each other with constant velocity `#mover(mi("V"),mo("→"))`. Show that the relation between the energies of the system, E and E', is given by

E(x) = diff(E(x), x)+(1/2)*LinearAlgebra[Norm](`#mover(mi("V"),mo("→"))`)^2*(Sum(m[a], a = 1 .. n))+`#mover(mi("V"),mo("→"))`.`#mover(mi("P'"),mo("→"))`

Solution

   

 

Conservation of the total momentum of a closed system

 


The conservation of the total momentum of a closed system of one particle is clear: if the particle does not interact with anything external, the force acting on it is zero, and from Newton's 2nd law `#mover(mi("F"),mo("→"))` = diff(`#mover(mi("p"),mo("→"))`(t), t) follows diff(`#mover(mi("p"),mo("→"))`(t), t) = 0.

 

For a closed system of many particles, while the total force acting on the system is equal to 0, there can be internal forces different from zero acting on each particle due to its interaction with the other particles. These internal forces, however, produce no acceleration of the system; in general, they cancel each other out due to Newton's 3rd law.

 

Problem

Consider a system of n particles measured in two frames of reference K and K'  that move relative to each other with velocity `#mover(mi("V"),mo("→"))`. Show that the system's momenta `#mover(mi("P"),mo("→"))` and "(P)^( ')" are related by

"P=(P)^( ')+V *(∑)m[a]".

Solution

   

 

Problem

A particle of mass m moving with velocity `#mover(mi("v"),mo("→"))`[1] leaves a half-space in which its potential energy is a constant U[1] and enters another in which its potential energy is a different constant U[2]. Determine the change in direction of motion of the particle; that is, sin(`θ__1`)/sin(`θ__2`) where `θ__1` and `θ__2` are the angles between an axis perpendicular to the separating plane and the momentum `#mover(mi("p"),mo("→"))` in the regions 1 and 2.

Solution

   

 

Conservation of angular momentum

 


Like the conservation of linear momentum, the conservation of the total angular momentum of a closed system of one particle is natural: if the particle does not interact with anything external, the force acting on it is zero and therefore its torque `#mover(mi("N"),mo("→"))` = `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("F"),mo("→"))`) and `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("F"),mo("→"))`) = 0. From diff(`#mover(mi("L"),mo("→"))`(t), t) = `#mover(mi("N"),mo("→"))` it follows that diff(`#mover(mi("L"),mo("→"))`(t), t) = 0, that is, the angular momentum `#mover(mi("L"),mo("→"))` is conserved.

 

Problem

a) Express the Cartesian components of the angular momentum `#mover(mi("L"),mo("→"))`, as well as its norm, in cylindrical and spherical coordinates.

b) Rewrite `#mover(mi("L"),mo("→"))` in cylindrical coordinates and using the cylindrical orthonormal basis of unit vectors, then do the same using spherical coordinates and the spherical basis.

Solution

   

 

Problem

Consider a system of n particles measured in two frames of reference K and K' whose origins have distance `#mover(mi("A"),mo("→"))` from each other. Show that the momenta `#mover(mi("L"),mo("→"))` and "((L))^( ')" of the system are related by

 "(L)=((L))^( ')+(A) *P"

where `#mover(mi("P"),mo("→"))` = Sum(`#mover(mi("p"),mo("→"))`[a], a = 1 .. n) is the total momentum of the system as seen from K.

Solution

   

 

Problem

a) Consider a closed system of n particles, and two frames of reference K and K'  that move relative to each other with a constant velocity `#mover(mi("V"),mo("→"))`. Show that the momenta `#mover(mi("L"),mo("→"))` and "((L))^( ')"respectively measured in K and K' are related by

`#mover(mi("L"),mo("→"))` = (Sum(m[a], a = 1 .. n))*`&x`(`#mover(mi("R"),mo("→"))`, `#mover(mi("V"),mo("→"))`)+`&x`(`#mover(mi("A"),mo("→"))`, `#mover(mi("P'"),mo("→"))`)+`#mover(mi("L'"),mo("→"))`

where `#mover(mi("A"),mo("→"))` is the distance from the origin of K to the origin of K',  "R=(∑)m[a] (r)[a]/(∑)m[a]" is the position of the center of mass as seen from K, and `#mover(mi("P'"),mo("→"))` and `#mover(mi("L'"),mo("→"))` are the total linear and angular momenta measured in K'.
b) Show that when the origin of K' is the center of mass `#mover(mi("R"),mo("→"))`, this formula reduces to

`#mover(mi("L"),mo("→"))` = `#mover(mi("L'"),mo("→"))`+`&x`(`#mover(mi("R"),mo("→"))`, `#mover(mi("P"),mo("→"))`),

where `#mover(mi("P"),mo("→"))` = Sum(m[a]*`#mover(mi("v"),mo("→"))`[a], a = 1 .. n) is the total linear momentum in K.

Solution

 

 

restart; with(Physics:-Vectors)

 

a) The momentum `#mover(mi("L"),mo("→"))` of the system measured in K is given by

L_ = Sum(`&x`(r_[a], p_[a]), a = 1 .. n)

L_ = Sum(Physics:-Vectors:-`&x`(r_[a], p_[a]), a = 1 .. n)

(322)

The right-hand side of the above can be expressed in terms of `#mover(mi("P'"),mo("→"))` and `#mover(mi("L'"),mo("→"))` using  `#mover(mi("p"),mo("→"))`[a] = m[a]*`#mover(mi("v"),mo("→"))`[a], where `#mover(mi("v"),mo("→"))`[a] = `#mover(mi("v'"),mo("→"))`[a]+`#mover(mi("V"),mo("→"))` :

p_[a] = (diff(p(x), x))*_[a]+m[a]*V_

p_[a] = V_*m[a]+`p'_`[a]

(323)

subs(p_[a] = V_*m[a]+`p'_`[a], L_ = Sum(Physics[Vectors][`&x`](r_[a], p_[a]), a = 1 .. n))

L_ = Sum(Physics:-Vectors:-`&x`(r_[a], V_*m[a]+`p'_`[a]), a = 1 .. n)

(324)

ee := expand(L_ = Sum(Physics[Vectors][`&x`](r_[a], V_*m[a]+`p'_`[a]), a = 1 .. n))

L_ = -Physics:-Vectors:-`&x`(V_, Sum(m[a]*r_[a], a = 1 .. n))+Sum(Physics:-Vectors:-`&x`(r_[a], `p'_`[a]), a = 1 .. n)

(325)

The term with Sum(m[a]*`#mover(mi("r"),mo("→"))`[a], a = 1 .. n) can be expressed in terms of the position vector of the center of mass `#mover(mi("R"),mo("→"))` (copy the subexpression from (325) , paste, then edit)

subs(Sum(m[a]*`#mover(mi("r"),mo("→"))`[a], a = 1 .. n) = (Sum(m[a], a = 1 .. n))*R_, L_ = -Physics[Vectors][`&x`](V_, Sum(m[a]*r_[a], a = 1 .. n))+Sum(Physics[Vectors][`&x`](r_[a], `p'_`[a]), a = 1 .. n))

L_ = -Physics:-Vectors:-`&x`(V_, (Sum(m[a], a = 1 .. n))*R_)+Sum(Physics:-Vectors:-`&x`(r_[a], `p'_`[a]), a = 1 .. n)

(326)

To express Sum(`&x`(`#mover(mi("r"),mo("→"))`[a], `#mover(mi("p'"),mo("→"))`[a]), a = 1 .. n) in terms of `#mover(mi("L'"),mo("→"))` and `#mover(mi("A"),mo("→"))`, introduce the relation between `#mover(mi("r"),mo("→"))`[a] and `#mover(mi("r'"),mo("→"))`[a]

r_[a] = A_+(diff(r(x), x))*_[a]

r_[a] = A_+`r'_`[a]

(327)

subs(r_[a] = A_+`r'_`[a], L_ = -Physics[Vectors][`&x`](V_, (Sum(m[a], a = 1 .. n))*R_)+Sum(Physics[Vectors][`&x`](r_[a], `p'_`[a]), a = 1 .. n))

L_ = -(Sum(m[a], a = 1 .. n))*Physics:-Vectors:-`&x`(V_, R_)+Sum(Physics:-Vectors:-`&x`(A_+`r'_`[a], `p'_`[a]), a = 1 .. n)

(328)

expand(L_ = -(Sum(m[a], a = 1 .. n))*Physics[Vectors][`&x`](V_, R_)+Sum(Physics[Vectors][`&x`](A_+`r'_`[a], `p'_`[a]), a = 1 .. n))

L_ = -(Sum(m[a], a = 1 .. n))*Physics:-Vectors:-`&x`(V_, R_)+Physics:-Vectors:-`&x`(A_, Sum(`p'_`[a], a = 1 .. n))+Sum(Physics:-Vectors:-`&x`(`r'_`[a], `p'_`[a]), a = 1 .. n)

(329)

On the right-hand side, two of the sums represent `#mover(mi("P'"),mo("→"))` and `#mover(mi("L'"),mo("→"))` (copy the sum subexpressions, paste into the next line, then edit)

subs(Sum(`&x`(`#mover(mi("r'"),mo("→"))`[a], `#mover(mi("p'"),mo("→"))`[a]), a = 1 .. n) = (diff(L(x), x))*_, Sum(`#mover(mi("p'"),mo("→"))`[a], a = 1 .. n) = (diff(P(x), x))*_, L_ = -(Sum(m[a], a = 1 .. n))*Physics[Vectors][`&x`](V_, R_)+Physics[Vectors][`&x`](A_, Sum(`p'_`[a], a = 1 .. n))+Sum(Physics[Vectors][`&x`](`r'_`[a], `p'_`[a]), a = 1 .. n))

L_ = -(Sum(m[a], a = 1 .. n))*Physics:-Vectors:-`&x`(V_, R_)+Physics:-Vectors:-`&x`(A_, `P'_`)+`L'_`

(330)

This is already the desired result.

_______________________________________

 

b) If the origin of K' is the center of mass `#mover(mi("R"),mo("→"))`, then `#mover(mi("A"),mo("→"))` = `#mover(mi("R"),mo("→"))`

subs(A_ = R_, L_ = -(Sum(m[a], a = 1 .. n))*Physics[Vectors][`&x`](V_, R_)+Physics[Vectors][`&x`](A_, `P'_`)+`L'_`)

L_ = -(Sum(m[a], a = 1 .. n))*Physics:-Vectors:-`&x`(V_, R_)+Physics:-Vectors:-`&x`(R_, `P'_`)+`L'_`

(331)

`#mover(mi("P"),mo("→"))` and `#mover(mi("P'"),mo("→"))` were related in "P_ = V_*Sum(m[a],a = 1 .. n)+`P'_`≡ "`#mover(mi("P"),mo("→"))` = `#mover(mi("V"),mo("→"))`*(Sum(m[a], a = 1 .. n))+`#mover(mi("P'"),mo("→"))`.  This relation can be used to express `#mover(mi("L"),mo("→"))` in terms of  `#mover(mi("P"),mo("→"))` instead of `#mover(mi("P'"),mo("→"))`

simplify(L_ = -(Sum(m[a], a = 1 .. n))*Physics[Vectors][`&x`](V_, R_)+Physics[Vectors][`&x`](R_, `P'_`)+`L'_`, {P_ = V_*(Sum(m[a], a = 1 .. n))+`P'_`}, {(diff(P(x), x))*_})

L_ = -(Sum(m[a], a = 1 .. n))*Physics:-Vectors:-`&x`(V_, R_)+Physics:-Vectors:-`&x`(R_, -V_*(Sum(m[a], a = 1 .. n))+P_)+`L'_`

(332)

expand(L_ = -(Sum(m[a], a = 1 .. n))*Physics[Vectors][`&x`](V_, R_)+Physics[Vectors][`&x`](R_, -V_*(Sum(m[a], a = 1 .. n))+P_)+`L'_`)

L_ = Physics:-Vectors:-`&x`(R_, P_)+`L'_`

(333)

NULL

which is the result we were looking for.

 

 

Cyclic coordinates

 

 

Any generalized coordinate q__i which does not appear explicitly in the Lagrangian is called cyclic. To any cyclic coordinate corresponds a conserved quantity. From

%diff(%diff(L, diff(q[i](t), t)), t) = %diff(L, q[i])

when q__i is cyclic, the right-hand side is 0 and so the quantity %diff(L, diff(q[i](t), t)) is conserved.

 

Problem

The Lagrangian describing the movement of a particle in a central field has phi as a cyclic coordinate. Using cylindrical coordinates, show that the corresponding conserved quantity is the z component of the angular momentum `#mover(mi("L"),mo("→"))`

Solution

   

 

Integration of the equations of motion

 

 

Motion in one dimension

 

Problem
For a closed system, or any system where the total energy E = T+U is conserved, show the following:

a) The trajectory in implicit form can always be computed directly from E.

b) The turning points, if any, can be computed directly from U.

Solution

   

 

Problem

Determine the period of oscillations of a pendulum of mass m and length l in a gravitational field as a function of the amplitude of the oscillations

Solution

   

 

Problem

Integrate the equations of motion for a particle of mass m moving in a field whose potential energy is U = A*abs(x)^n.

Solution

   

Reduced mass

 

The two-body problem

 

 

Problem
Show that by placing the origin of the reference system at the center of mass  `#mover(mi("R"),mo("→"))` = (Sum(m[i]*`#mover(mi("r"),mo("→"))`[i], i = 1 .. n))/(Sum(m[i], i = 1 .. n)), the problem of two particles that interact with each other can be reduced to the problem of a single particle of mass mu = m__1*m__2/(m__1+m__2), herein called reduced mass, in an external field U(LinearAlgebra[Norm](r__2_-r__1_)).

Solution

   

 

A many-body problem

 


Problem

A system consists of one particle of mass M and n particles of equal masses m.

a) Show, in steps, that eliminating the motion of the center of mass reduces the problem to one involving only n particles.

b) Show that when n = 1the result of a) becomes the result obtained for the previous two-body problem, equation `≡`(L = (1/2)*Physics[Vectors][Norm](diff(`ℝ_`(t), t))^2*mu-U(Physics[Vectors][Norm](`ℝ_`(t))), L) = (1/2)*mu*LinearAlgebra[Norm](diff(`#mover(mi("ℝ"),mo("→"))`(t), t))^2-U.

Solution

 

 

restart; with(Physics:-Vectors)

 

a) Let R__M_ represent the position vector of the particle of mass M and `#mover(mi("r"),mo("→"))`[a] those of the particles of mass m. The Lagrangian L = T-U is given by

CompactDisplay((r__M_, r_, `ℝ_`)(t))

`ℝ_`(t)*`will now be displayed as`*`ℝ_`

(386)

Physics:-Assume(m > 0, M > 0)

{m::(RealRange(Open(0), infinity))}, {M::(RealRange(Open(0), infinity))}

(387)

L = (1/2)*M*(diff(r__M_(t), t))^2+(1/2)*(Sum(m*(diff(r_[a](t), t))^2, a = 1 .. n))-U

L = (1/2)*M*Physics:-`^`(diff(r__M_(t), t), 2)+(1/2)*(Sum(m*Physics:-`^`(diff(r_[a](t), t), 2), a = 1 .. n))-U

(388)

Expand the formal powers of vectors to express them in terms of their Norm

expand(L = (1/2)*M*Physics[`^`](diff(r__M_(t), t), 2)+(1/2)*(Sum(m*Physics[`^`](diff(r_[a](t), t), 2), a = 1 .. n))-U)

L = (1/2)*M*Physics:-Vectors:-Norm(diff(r__M_(t), t))^2+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(r_[a](t), t))^2, a = 1 .. n))-U

(389)

As done in the two-body problem, introduce the relative vector positions of the n particles with respect to the particle of mass M

`ℝ_`[a](t) = r_[a](t)-r__M_(t)

`ℝ_`[a](t) = r_[a](t)-r__M_(t)

(390)

and place the origin of the reference system at the center of mass,

(M*r__M_(t)+Sum(m*r_[a](t), a = 1 .. n))/(M+Sum(m, a = 1 .. n)) = 0

(M*r__M_(t)+Sum(m*r_[a](t), a = 1 .. n))/(M+Sum(m, a = 1 .. n)) = 0

(391)

expand((M*r__M_(t)+Sum(m*r_[a](t), a = 1 .. n))/(M+Sum(m, a = 1 .. n)) = 0)

M*r__M_(t)/(m*n+M)+m*(Sum(r_[a](t), a = 1 .. n))/(m*n+M) = 0

(392)

Using this equation and (390)`#mover(mi("r"),mo("→"))`[a] = `#mover(mi("ℝ"),mo("→"))`[a]-R__M_, fontstyle = "normal", Typesetting:-msemantics = "function" is sufficient to eliminate R__M_ and `#mover(mi("r"),mo("→"))`[a] from the Lagrangian (389). In steps, eliminating `#mover(mi("r"),mo("→"))`[a] from the equation for the center of mass

simplify(M*r__M_(t)/(m*n+M)+m*(Sum(r_[a](t), a = 1 .. n))/(m*n+M) = 0, {`ℝ_`[a](t) = r_[a](t)-r__M_(t)}, {r_[a](t)})

(M*r__M_(t)+m*(Sum(`ℝ_`[a](t)+r__M_(t), a = 1 .. n)))/(m*n+M) = 0

(393)

expand((M*r__M_(t)+m*(Sum(`ℝ_`[a](t)+r__M_(t), a = 1 .. n)))/(m*n+M) = 0)

M*r__M_(t)/(m*n+M)+m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M)+m*r__M_(t)*n/(m*n+M) = 0

(394)

To eliminate `#msub(mi("R"),mi("M_"))` , use

simplify(isolate(M*r__M_(t)/(m*n+M)+m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M)+m*r__M_(t)*n/(m*n+M) = 0, r__M_(t)))

r__M_(t) = -m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M)

(395)

Likewise, to eliminate R__M_ use

isolate(`ℝ_`[a](t) = r_[a](t)-r__M_(t), r_[a](t))

r_[a](t) = `ℝ_`[a](t)+r__M_(t)

(396)

Substitute, sequentially, these two equations into the Lagrangian (389)

subs(r_[a](t) = `ℝ_`[a](t)+r__M_(t), r__M_(t) = -m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M), L = (1/2)*M*Physics[Vectors][Norm](diff(r__M_(t), t))^2+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(r_[a](t), t))^2, a = 1 .. n))-U)

L = (1/2)*M*Physics:-Vectors:-Norm(diff(-m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M), t))^2+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t)-m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M), t))^2, a = 1 .. n))-U

(397)

To verify by eye each step for correctness, one can manipulate this expression surgically using subsindets . First expand only the Norms 

subsindets(L = (1/2)*M*Physics[Vectors][Norm](diff(-m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M), t))^2+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t)-m*(Sum(`ℝ_`[a](t), a = 1 .. n))/(m*n+M), t))^2, a = 1 .. n))-U, specfunc(Norm), expand)

L = (1/2)*M*m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2-2*m*Physics:-Vectors:-`.`(diff(`ℝ_`[a](t), t), Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))/(m*n+M)+m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2, a = 1 .. n))-U

(398)

This result is correct. Next expand only the Sums 

subsindets(L = (1/2)*M*m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2-2*m*Physics[Vectors][`.`](diff(`ℝ_`[a](t), t), Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))/(m*n+M)+m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2, a = 1 .. n))-U, specfunc(Sum), expand)

L = (1/2)*M*m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2+(1/2)*m*(m^2*n^2*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2+2*M*m*n*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2+M^2*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2-m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2*n/(m*n+M)^2-2*M*m*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2)-U

(399)

This result also verifies visually. Collect the terms polynomial in Norm then Sum and normalize the coefficients

collect(L = (1/2)*M*m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2+(1/2)*m*(m^2*n^2*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2+2*M*m*n*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2+M^2*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. n))/(m*n+M)^2-m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2*n/(m*n+M)^2-2*M*m*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2)-U, [Norm, Sum], normal)

L = -(1/2)*m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2, a = 1 .. n))-U

(400)

This Lagrangian involves only the n relative position vectors `#mover(mi("ℝ"),mo("→"))`[a], achieving the reduction of the original problem of n+1 particles to a problem of only n particles.

_______________________________________

 

b) To obtain the result `≡`(L = (1/2)*Physics[Vectors][Norm](diff(`ℝ_`(t), t))^2*mu-U(Physics[Vectors][Norm](`ℝ_`(t))), L) = (1/2)*LinearAlgebra[Norm](diff(`#mover(mi("ℝ"),mo("→"))`(t), t))^2*mu-U for the two-body problem, substitute n = 1 into `≡`(L = -(1/2)*m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. n))-U, L) = -(1/2)*m^2*LinearAlgebra[Norm](Sum(diff(`#mover(mi("ℝ"),mo("→"))`[a](t), t), a = 1 .. n))^2/(m*n+M)+(1/2)*m*(Sum(LinearAlgebra[Norm](diff(`#mover(mi("ℝ",fontstyle = "normal"),mo("→"))`[a](t), t))^2, a = 1 .. n))-U

subs(n = 1, L = -(1/2)*m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. n))^2/(m*n+M)+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. n))-U)

L = -(1/2)*m^2*Physics:-Vectors:-Norm(Sum(diff(`ℝ_`[a](t), t), a = 1 .. 1))^2/(M+m)+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`ℝ_`[a](t), t))^2, a = 1 .. 1))-U

(401)

Expand only the Sums and collect the Norm

subsindets(L = -(1/2)*m^2*Physics[Vectors][Norm](Sum(diff(`ℝ_`[a](t), t), a = 1 .. 1))^2/(M+m)+(1/2)*m*(Sum(Physics[Vectors][Norm](diff(`ℝ_`[a](t), t))^2, a = 1 .. 1))-U, specfunc(Sum), expand)

L = -(1/2)*m^2*Physics:-Vectors:-Norm(diff(`ℝ_`[1](t), t))^2/(M+m)+(1/2)*m*Physics:-Vectors:-Norm(diff(`ℝ_`[1](t), t))^2-U

(402)

collect(L = -(1/2)*m^2*Physics[Vectors][Norm](diff(`ℝ_`[1](t), t))^2/(M+m)+(1/2)*m*Physics[Vectors][Norm](diff(`ℝ_`[1](t), t))^2-U, Norm, normal)

L = (1/2)*M*m*Physics:-Vectors:-Norm(diff(`ℝ_`[1](t), t))^2/(M+m)-U

(403)

Comparing with the definition of reduced mass (384)

m[1]*m[2]/(m[1]+m[2]) = mu

m[1]*m[2]/(m[1]+m[2]) = mu

(404)

we see the substitutions to transform (403) into (384) are

M = m[1], m = m[2], r_[1] = R_

M = m[1], m = m[2], r_[1] = R_

(405)

Substitute them, simultaneously (enclose the sequence of equations into a list), then simplify taking m[1]*m[2]/(m[1]+m[2]) = mu into account

simplify(subs([M = m[1], m = m[2], r_[1] = R_], L = (1/2)*M*m*Physics[Vectors][Norm](diff(`ℝ_`[1](t), t))^2/(M+m)-U), {m[1]*m[2]/(m[1]+m[2]) = mu})

L = (1/2)*Physics:-Vectors:-Norm(diff(`ℝ_`[1](t), t))^2*mu-U

(406)

NULL

This is the same as (385), the reduced Lagrangian for the two-body problem.

 

Motion in a central field

 

A one-body problem in a central field, is about the motion of a single particle of mass m in a field where the potential energy, and so the magnitude of the force, depends only on the distance between the particle and a fixed point, frequently chosen as the origin of the reference system. As seen above, the two-body problem , is reducible to a one-body problem in a central field.

 

Problem

The angular momentum `#mover(mi("L"),mo("→"))` = `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("p"),mo("→"))`) of a particle that moves in a central field is conserved, so `#mover(mi("r"),mo("→"))` evolves in time on a fixed plane perpendicular to the constant `#mover(mi("L"),mo("→"))`. Show that the surface swept per second by the position vector `#mover(mi("r"),mo("→"))` is constant (Kepler's second law).

Solution

   

 

Problem
Starting from the constancy of the energy E and the angular momentum `#mover(mi("L"),mo("→"))`, compute the equations of movement and integrate them according to:

a) using differential elimination techniques to uncouple the system of equations of movement that involve both of rho(t) and phi(t) 

b) interactively, one step at a time, uncouple the equations of movement eliminating phi from the problem, resulting in an implicit solution `≡`(t, t(rho)).

c) eliminate t from the problem to obtain an equation for diff(phi(rho), rho), whose solution is the trajectory as phi(rho)

Solution

   

 

Kepler's problem

 

 

Problem

Show that, when U(rho) = -alpha/rho, with "alpha>0,"the solution (442) for the motion in a central field, part c) of the previous problem,

phi(rho) = Int(L__z/(rho^2*(2*m*(E-U(rho))-L__z^2/rho^2)^(1/2)), rho)+c__1

becomes the equation of a conic section

p/rho = 1+xi*cos(phi) 

where p = L__z^2/(m*alpha), xi = sqrt(1+2*E*L__z^2/(m*alpha^2)).

Solution

   

 

Small Oscillations

 

Free oscillations in one dimension

 

 

Problem

Consider the case of 1-dimensional motion where the acting force opposes the motion as a function of the position, `#mover(mi("F"),mo("→"))` = -k*x*`#mover(mi("i"),mo("∧"))`. This is the case, for example, of a spring, the more you stretch it (the bigger x), the more it opposes the stretching in the opposite direction -`#mover(mi("i"),mo("∧"))`. Write the equation of motion as Newton's 2nd law, then the Lagrangian and Lagrange equation for it, and integrate the equation for generic initial conditions

Solution

   

 

Forced oscillations

 

 

Problem

Consider oscillations in 1 dimension of a system on which an external force `#mover(mi("F"),mo("→"))`[ext](t) acts. For the oscillations to be small, `#mover(mi("F"),mo("→"))`[ext](t) must produce only small displacements. The total force is  `#mover(mi("F"),mo("→"))` = -k*x*`#mover(mi("i"),mo("∧"))`+`#mover(mi("F"),mo("→"))`[ext](t).

a) Write the equation of motion as Newton's 2nd law, then write the Lagrangian and Lagrange equations.

b) Integrate the equation of motion for generic initial conditions.

c) Specialize the solution computed in b) for LinearAlgebra[Norm](`#mover(mi("F"),mo("→"))`[ext](t)) = f__0*cos(lambda*t+beta) to obtain

x(t) = (a*cos(omega*t+alpha)-f__0*omega*cos(lambda*t+beta)/(m*(lambda^2-omega^2))-f__0*cos(-omega*t+beta)/(2*m*(lambda+omega))+f__0*cos(omega*t+beta)/((2*(lambda-omega))*m))/omega

for some constants a, b, alpha, beta.

d) Show that a solution for the case considered in c), that is, LinearAlgebra[Norm](`#mover(mi("F"),mo("→"))`[ext](t)) = f__0*cos(lambda*t+beta), can be computed manually to get

x(t) = a*cos(omega*t+alpha)+f__0*cos(lambda*t+beta)/(m*(-lambda^2+omega^2))

for some other constants a, b, alpha, beta.

e) Specialize the solution of item d) in the case of resonance, when lambda = omega, by taking limits, thus obtaining

x(t) = a*cos(omega*t+alpha)+f__0*t*sin(omega*t+beta)/(2*omega*m)

f) Show that the solution computed taking limits in e) can be compute directly by using dsolve and specializing the integration constants c__1 and c__2 that appear when solving the underlying differential equation.

Solution

   

 

Oscillations of systems with many degrees of freedom

 

 

Problem
Formulate the equations of motion for the free oscillations of a system with n degrees of freedom as

m[a, c]*(diff(x[a](t), t, t))+`k__a,c`*x[a] = 0

where x[a] represents the displacement of the a^th generalized coordinate q__a, the index a runs from 1 to n and there is an implicit sum over repeated indices (Einstein's convention).

Solution

 

 

restart; with(Physics); with(Vectors); CompactDisplay(x(t))

x(t)*`will now be displayed as`*x

(524)

 

Denoting the generalized coordinates by q__a, the potential energy U(q__a) can be expanded in series around the minimums q__a0. Since the movement consist only in small displacements x__a = q__a-q__a0 around q__a0, it is sufficient to keep terms in the expansion up to order 2, resulting in an expression analogous to U = (1/2)*k*x^2 of the 1-dimensional case:

U = (1/2)*(Sum(Sum(k[a, b]*x[a](t)*x[b](t), i = 1 .. n), j = 1 .. n))

U = (1/2)*(Sum(Sum(k[a, b]*x[a](t)*x[b](t), i = 1 .. n), j = 1 .. n))

(525)

Likewise, for the kinetic energy,

T = (1/2)*(Sum(Sum(A[a, b](q__0)*(diff(x[a](t), t))*(diff(x[b](t), t)), a = 1 .. n), b = 1 .. n))

T = (1/2)*(Sum(Sum(A[a, b](q__0)*(diff(x[a](t), t))*(diff(x[b](t), t)), a = 1 .. n), b = 1 .. n))

(526)

where we take the A[a, b] at the minimums `#msub(mi("q"),mn("0"))` and denote them as `m__i,j`

subs(A[a, b](q__0) = m[a, b], T = (1/2)*(Sum(Sum(A[a, b](q__0)*(diff(x[a](t), t))*(diff(x[b](t), t)), a = 1 .. n), b = 1 .. n)))

T = (1/2)*(Sum(Sum(m[a, b]*(diff(x[a](t), t))*(diff(x[b](t), t)), a = 1 .. n), b = 1 .. n))

(527)

Both `k__a,b` and `m__a,b` can be split into symmetric and antisymmetric parts, with the antisymmetric parts canceling out in view of the symmetric character of x[a]*x[b] in U and (diff(x[a](t), t))*(diff(x[b](t), t)) in T. Therefore, we can take `k__a,b` and `m__a,b` as symmetric without any loss of generality.

 

Before proceeding, note the similarity in notation between the three formulas (525) to (527) for T and U and tensor notation. In T and U the x__a describe independent displacements, so one can think of `#msub(mi("x"),mi("a"))` as a tensor in an Euclidean space of displacements of generic abstract dimension n, with KroneckerDelta  as the metric. It is then simpler to write the Lagrangian using tensor notation with a generic type of index (that admits and abstract n-dimension). For this purpose, introduce lowercaselatin indices from a to h to represent generic indices, and when necessary use KroneckerDelta  as the metric.

Setup(genericindices = lowercaselatin_ah)

[genericindices = lowercaselatin_ah]

(528)

Now introduce the tensors while making sure to indicate `m__a,b` and `k__a,b` are symmetric (passing x__a together is not a problem, it has only one index)

Define(x[a], m[a, b], k[a, b], symmetric)

{Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], k[a, b], m[a, b], x[a], Physics:-LeviCivita[alpha, beta, mu, nu]}

(529)

The Lagrangian L = T-U can now be written as

L := 1/2*(m[a, b]*(diff(x[a](t), t))*(diff(x[b](t), t))-k[a, b]*x[a](t)*x[b](t))

(1/2)*m[a, b]*(diff(x[a](t), t))*(diff(x[b](t), t))-(1/2)*k[a, b]*x[a](t)*x[b](t)

(530)

where Einstein's summation rule for repeated indices is used. Einstein's rule is taken into account by the system when differentiating, computing products and simplifying tensor indices. The simplest way to compute the Lagrange equations is to use the LagrangeEquations  command

LagrangeEquations(L, x)

k[a, c]*x[a](t)+m[a, c]*(diff(diff(x[a](t), t), t)) = 0

(531)

The same result can be computed via %diff(%diff(L, diff(x[c](t), t)), t) = %diff(L, x[c](t)). Although not necessary, enclose the operation with forward quotes to delay its evaluation in order to see what is being computed

'%diff(%diff(L, diff(x[c](t), t)), t) = %diff(L, x[c](t))'

%diff(%diff(L, Physics:-Vectors:-diff(x[c](t), t)), t) = %diff(L, x[c](t))

(532)

value(%diff(%diff(L, Physics[Vectors][diff](x[c](t), t)), t) = %diff(L, x[c](t)))

(1/2)*m[a, b]*(diff(diff(x[a](t), t), t))*Physics:-KroneckerDelta[b, c]+(1/2)*m[a, b]*(diff(diff(x[b](t), t), t))*Physics:-KroneckerDelta[a, c] = -(1/2)*k[a, b]*x[a](t)*Physics:-KroneckerDelta[b, c]-(1/2)*k[a, b]*x[b](t)*Physics:-KroneckerDelta[a, c]

(533)

Simplifying tensor indices,

Simplify((1/2)*m[a, b]*(diff(diff(x[a](t), t), t))*Physics[KroneckerDelta][b, c]+(1/2)*m[a, b]*(diff(diff(x[b](t), t), t))*Physics[KroneckerDelta][a, c] = -(1/2)*k[a, b]*x[a](t)*Physics[KroneckerDelta][b, c]-(1/2)*k[a, b]*x[b](t)*Physics[KroneckerDelta][a, c])

m[a, c]*(diff(diff(x[a](t), t), t)) = -k[a, c]*x[a](t)

(534)

NULL

which is the same as (531).

``

 

Rigid-body motion

 

 

A rigid body is one where (in approximation) the distances between the body's parts remain unchanged. In what follows, for simplicity, the body is considered as discrete set of particles; the formulas for a continuous body can be obtained from those by replacing the masses m__i of each particle by rho(`#mover(mi("r"),mo("→"))`)*dV, where rho(`#mover(mi("r"),mo("→"))`) is the mass density as a function of the position and dV is the volume element, whose integration represent the body's volume.

 

This problem is systematically treated by using two reference systems: an inertial one K, where the observer is, and another one K', rigidly attached to the body, that moves with it and thus it is typically non-inertial. It is customary (not necessary) to place the origin `#mover(mi("R"),mo("→"))`(t) of K'  at the body's center of mass .

 

A rigid body is thus a system with six degrees of freedom: three indicating the position `#mover(mi("R"),mo("→"))`(t) of the center of mass plus three angles specifying the orientation of the axes of K' with respect to those of K.

 

Angular velocity

 

 

Problem
a) Show, using graphs, that the velocity `#mover(mi("v"),mo("→"))` of a point P of a body, measured in an inertial reference system K, can be written as

`#mover(mi("v"),mo("→"))` = `#mover(mi("V"),mo("→"))`+`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `#mover(mi("r'"),mo("→"))`)

where `#mover(mi("V"),mo("→"))` is the velocity of the origin of K', a frame of reference attached to the body's center of mass, `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` = (diff(Phi(t), t))*`#mover(mi("Φ",fontstyle = "normal"),mo("∧"))` is the body's angular velocity (its instantaneous counter-clockwise rotation speed around some axis in the direction of a unit vector`#mover(mi("Φ",fontstyle = "normal"),mo("∧"))`) and `#mover(mi("r'"),mo("→"))`  is the distance from the center of mass (origin of K') to the point P.

b) Derive algebraically the same result of a) , using the fact that vectors are defined up to parallel translation and so `#mover(mi("r"),mo("→"))` and `#mover(mi("r'"),mo("→"))` are related by a rotation matrix omega[a, b] which, as all rotation matrices, is orthogonal.

Solution

   

 

Inertia tensor

 

 

Problem

a) Show that using `#mover(mi("v"),mo("→"))` = `#mover(mi("V"),mo("→"))`+`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `#mover(mi("r'"),mo("→"))`), derived in Angular velocity  for the velocity `#mover(mi("v"),mo("→"))` of a point P of a rigid body in terms of `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` and the position `#mover(mi("r'"),mo("→"))` of P viewed from the center of mass `#mover(mi("R"),mo("→"))`, the kinetic energy of the rigid body can be written in terms of the positions `#mover(mi("r'"),mo("→"))` of the n particles (not their velocities), the velocity of the center of mass `#mover(mi("V"),mo("→"))` and angular velocity `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` as

T = (1/2)*LinearAlgebra[Norm](`#mover(mi("V"),mo("→"))`)^2*mu+Sum((1/2)*m[i]*LinearAlgebra[Norm](`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`)^2*LinearAlgebra[Norm](`#mover(mi("r'"),mo("→"))`[i])^2-(1/2)*m[i]*(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`.`#mover(mi("r'"),mo("→"))`[i])^2, i = 1 .. n)

b) Use tensor notation to show that this result can be rewritten as

T = (1/2)*LinearAlgebra[Norm](`#mover(mi("V"),mo("→"))`(t))^2*mu+(1/2)*`𝕀`[a, b]*Omega[b]*Omega[a]

where

"`𝕀`[a,b]=(∑)m[i] (r[i,c]^2 delta[a,b]-`r'`[i,a] `r'`[i,b])"

is the Inertia tensor, `Ω__a` represents the components of the vector `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` and `r'`[i, a] represents the components of the position vector `#mover(mi("r'"),mo("→"))`[i] of the i^th particle.

Solution

   

 

NULL

Problem

Determine the Inertia tensor corresponding to a triatomic molecule that has the form of an isosceles triangle with two masses m[1] in the extremes of the base and a mass m[2] in the third vertex. The distance between the two masses m[1] is equal to a, and the height of the triangle is equal to h.

Solution

   

 

Angular momentum of a rigid body

 

 

In the section related to the conservation of angular momentum , the solution to the second Problem  shows that the value of the angular momentum `#mover(mi("L"),mo("→"))` of a system of particles depends on the origin of the frame of reference. In this section, it is assumed that the origin is at the center of mass, so Sum(m[i]*`#mover(mi("r"),mo("→"))`[i], i = 1 .. n) = 0.

 

Problem
Show, using tensor notation, that in the K' system whose origin is at the center of mass, the components (diff(L(x), x))[a] of the angular momentum of a rigid body can be expressed in terms of the inertia tensor `𝕀`[a, b] and the components of the angular velocity Omega[b] as

(diff(L(x), x))[a] = `𝕀`[a, b]*Omega[b]

Solution

   

 

The equations of motion of a rigid body

 

 

A rigid body is a system with six degrees of freedom: three indicating the position `#mover(mi("R"),mo("→"))`(t) plus three angles specifying the orientation of the axes of K' with respect to those of K. As discussed in the equations of motion  for many-particle systems , the two vectorial equations of motion are

value(%diff(P_(t), t) = F_[ext])

diff(P_(t), t) = F_[ext]

(615)

diff(L_(t), t) = N_(t)

diff(L_(t), t) = N_(t)

(616)

N_ = Sum(Physics[Vectors][`&x`](r_[i](t), f_[i, ext]), i = 1 .. n)

N_ = Sum(Physics:-Vectors:-`&x`(r_[i](t), f_[i, ext]), i = 1 .. n)

(617)

NULL

where `#mover(mi("P"),mo("→"))` is the total momentum, `#mover(mi("F"),mo("→"))`[ext] is the total external force acting upon the body, `#mover(mi("f"),mo("→"))`[i, ext] is the external force acting upon the i^th particle, `#mover(mi("L"),mo("→"))`the total angular momentum and `#mover(mi("N"),mo("→"))` is the total torque.

 

Problem
Show that the equations of movement of a rigid body can be computed as the Lagrange equations for `#mover(mi("R"),mo("→"))` and `#mover(mi("Φ",fontstyle = "normal"),mo("→"))` from the Lagrangian

L = (1/2)*mu*LinearAlgebra[Norm](`#mover(mi("V"),mo("→"))`)^2+(1/2)*`𝕀`[a, b]*Omega[b]*Omega[a]+U(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`)

where mu = Sum(m[i], i = 1 .. n) is the total mass, `𝕀`[a, b] is the inertia tensor, U(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`) is the potential energy for the external force `#mover(mi("F"),mo("→"))`[ext], and  `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` = diff(`#mover(mi("Φ",fontstyle = "normal"),mo("→"))`(t), t).

Solution

 

 

restart; with(Physics:-Vectors); CompactDisplay((R_, V_, Phi_, Omega)(t))

Omega(t)*`will now be displayed as`*Omega

(618)

The required derivation is easy by expressing the Lagrangian in tensor notation from the beginning. In what follows, however, with the purpose of illustrating different techniques, Lagrange's equations are computed using vectorial notation, only switching to tensor notation at the time of expressing the time derivative of the angular momentum.

 

The kinetic energy T in vectorial form is derived in this problem  for the inertia tensor

subs((diff(r(x), x))*_[i](t) = (diff(r(x), x))*_[i], T = (1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i](t))^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i](t))^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu)

T = (1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu

(619)

Adding the potential energy as a function of the coordinates `#mover(mi("R"),mo("→"))` (location of the center of mass) and `#mover(mi("Φ",fontstyle = "normal"),mo("→"))` (the rotation axis, so that `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` = diff(`#mover(mi("Φ",fontstyle = "normal"),mo("→"))`(t), t)), the Lagrangian in vectorial form is given by

L = rhs(T = (1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i])^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu)+U(R_(t), Phi_(t))

L = (1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu+U(R_(t), Phi_(t))

(620)

The first equation of movement, for the total momentum `#mover(mi("P"),mo("→"))`(t) is derived as the Lagrange equation for this Lagrangian

%diff(%diff(L, V_(t)), t) = %diff(L, R_(t))

%diff(%diff(L, V_(t)), t) = %diff(L, R_(t))

(621)

subs(L = (1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i])^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu+U(R_(t), Phi_(t)), %diff(%diff(L, V_(t)), t) = %diff(L, R_(t)))

%diff(%diff((1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu+U(R_(t), Phi_(t)), V_(t)), t) = %diff((1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu+U(R_(t), Phi_(t)), R_(t))

(622)

value(%)

mu*(diff(V_(t), t)) = (D[1](U))(R_(t), Phi_(t))

(623)

On the right-hand side, (D[1](U))(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`)is*the*derivative*ofU(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`) with respect to its first argument `#mover(mi("R"),mo("→"))`, equivalent to the gradient taking `#mover(mi("R"),mo("→"))` as the coordinates, and mu*(diff(`#mover(mi("V"),mo("→"))`(t), t)) = diff(`#mover(mi("P"),mo("→"))`(t), t) is the total momentum

subs(rhs(mu*(diff(V_(t), t)) = (D[1](U))(R_(t), Phi_(t))) = F_[ext](t), diff(V_(t), t) = (diff(P_(t), t))/mu, mu*(diff(V_(t), t)) = (D[1](U))(R_(t), Phi_(t)))

diff(P_(t), t) = F_[ext](t)

(624)

Note however that in Maple the Vectors:-Gradient  command computes the gradient with respect to Cartesian, cylindrical or spherical coordinates, not an arbitrary vector `#mover(mi("R"),mo("→"))`. If more precision is required, the dependency of the potential energy could be expressed in terms of the norm of `#mover(mi("R"),mo("→"))`, as in

U(Norm(R_))

U(Physics:-Vectors:-Norm(R_))

(625)

in which case differentiating with respect to `#mover(mi("R"),mo("→"))` is equivalent to the definition of directional derivative

(%diff = diff)(U(Physics[Vectors][Norm](R_)), `#mover(mi("R"),mo("→"))`)

%diff(U(Physics:-Vectors:-Norm(R_)), R_) = (D(U))(Physics:-Vectors:-Norm(R_))*R_/Physics:-Vectors:-Norm(R_)

(626)

The same computation, this time with respect to the coordinates `#mover(mi("Φ",fontstyle = "normal"),mo("→"))` where `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` = diff(`#mover(mi("Φ",fontstyle = "normal"),mo("→"))`(t), t) is the corresponding velocity,

%diff(%diff(L, Omega_(t)), t) = %diff(L, Phi_(t))

%diff(%diff(L, Omega_(t)), t) = %diff(L, Phi_(t))

(627)

subs(L = (1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i])^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu+U(R_(t), Phi_(t)), %diff(%diff(L, Omega_(t)), t) = %diff(L, Phi_(t)))

%diff(%diff((1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu+U(R_(t), Phi_(t)), Omega_(t)), t) = %diff((1/2)*(Sum(m[i]*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(`r'_`[i])^2-m[i]*Physics:-Vectors:-`.`(Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics:-Vectors:-Norm(V_(t))^2*mu+U(R_(t), Phi_(t)), Phi_(t))

(628)

subs(sum = Sum, value(%diff(%diff((1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i])^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu+U(R_(t), Phi_(t)), Omega_(t)), t) = %diff((1/2)*(Sum(m[i]*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](`r'_`[i])^2-m[i]*Physics[Vectors][`.`](Omega_(t), `r'_`[i])^2, i = 1 .. n))+(1/2)*Physics[Vectors][Norm](V_(t))^2*mu+U(R_(t), Phi_(t)), Phi_(t))))

(1/2)*(Sum(2*m[i]*Physics:-Vectors:-Norm(`r'_`[i])^2*(diff(Omega_(t), t))-2*m[i]*Physics:-Vectors:-`.`(diff(Omega_(t), t), `r'_`[i])*`r'_`[i], i = 1 .. n)) = (D[2](U))(R_(t), Phi_(t))

(629)

Switching to tensor notation as done in the other problems, the left-hand side can be rewritten as the time derivative of the components (diff(L(x), x))[a] of the angular momentum .

with(Physics)

Setup(spaceindices = lowercase_ah, generic = lowercase_is, tensors = {Omega[a], (diff(r(x), x))[i, a]})

[genericindices = lowercaselatin_is, spaceindices = lowercaselatin_ah, tensors = {Physics:-Dgamma[mu], Omega[a], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[a, b], `r'`[i, a], Physics:-LeviCivita[alpha, beta, mu, nu]}]

(630)

Introducing tensor components for the vectors of (629)

Typesetting[delayDotProduct]((diff(Omega_(t), t)).(diff(r(x), x)), _[i], true) = (diff(Omega[b](t), t))*(diff(r(x), x))[i, b], diff(Omega_(t), t) = KroneckerDelta[a, b]*(diff(Omega[b](t), t)), Norm((diff(r(x), x))*_[i])^2 = (diff(r(x), x))[i, c]^2, (diff(r(x), x))*_[i] = (diff(r(x), x))[i, a], (D[2](U))(R_(t), Phi_(t)) = Component((D[2](U))(R_(t), Phi_(t)), a)

Physics:-Vectors:-`.`(diff(Omega_(t), t), `r'_`[i]) = (diff(Omega[b](t), t))*`r'`[i, b], diff(Omega_(t), t) = Physics:-KroneckerDelta[a, b]*(diff(Omega[b](t), t)), Physics:-Vectors:-Norm(`r'_`[i])^2 = `r'`[i, c]^2, `r'_`[i] = `r'`[i, a], (D[2](U))(R_(t), Phi_(t)) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(631)

eval((1/2)*(Sum(2*m[i]*Physics[Vectors][Norm](`r'_`[i])^2*(diff(Omega_(t), t))-2*m[i]*Physics[Vectors][`.`](diff(Omega_(t), t), `r'_`[i])*`r'_`[i], i = 1 .. n)) = (D[2](U))(R_(t), Phi_(t)), [Physics[Vectors][`.`](diff(Omega_(t), t), `r'_`[i]) = (diff(Omega[b](t), t))*`r'`[i, b], diff(Omega_(t), t) = Physics[KroneckerDelta][a, b]*(diff(Omega[b](t), t)), Physics[Vectors][Norm](`r'_`[i])^2 = `r'`[i, c]^2, `r'_`[i] = `r'`[i, a], (D[2](U))(R_(t), Phi_(t)) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a)])

(1/2)*(Sum(2*m[i]*`r'`[i, c]^2*Physics:-KroneckerDelta[a, b]*(diff(Omega[b](t), t))-2*m[i]*(diff(Omega[b](t), t))*`r'`[i, b]*`r'`[i, a], i = 1 .. n)) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(632)

expand((1/2)*(Sum(2*m[i]*`r'`[i, c]^2*Physics[KroneckerDelta][a, b]*(diff(Omega[b](t), t))-2*m[i]*(diff(Omega[b](t), t))*`r'`[i, b]*`r'`[i, a], i = 1 .. n)) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a))

Physics:-KroneckerDelta[a, b]*(diff(Omega[b](t), t))*(Sum(m[i]*`r'`[i, c]^2, i = 1 .. n))-(diff(Omega[b](t), t))*(Sum(m[i]*`r'`[i, a]*`r'`[i, b], i = 1 .. n)) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(633)

collect(Physics[KroneckerDelta][a, b]*(diff(Omega[b](t), t))*(Sum(m[i]*`r'`[i, c]^2, i = 1 .. n))-(diff(Omega[b](t), t))*(Sum(m[i]*`r'`[i, a]*`r'`[i, b], i = 1 .. n)) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a), diff(Omega[b](t), t), `@`(factor, combine))

(Sum(-m[i]*(-`r'`[i, c]^2*Physics:-KroneckerDelta[a, b]+`r'`[i, a]*`r'`[i, b]), i = 1 .. n))*(diff(Omega[b](t), t)) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(634)

Introducing the expression for the inertia tensor,

`𝕀`[a, b] = Sum(-m[i]*(-`r'`[i, c]^2*Physics[KroneckerDelta][a, b]+`r'`[i, a]*`r'`[i, b]), i = 1 .. n)

`𝕀`[a, b] = Sum(-m[i]*(-`r'`[i, c]^2*Physics:-KroneckerDelta[a, b]+`r'`[i, a]*`r'`[i, b]), i = 1 .. n)

(635)

subs((rhs = lhs)(`𝕀`[a, b] = Sum(-m[i]*(-`r'`[i, c]^2*Physics[KroneckerDelta][a, b]+`r'`[i, a]*`r'`[i, b]), i = 1 .. n)), (Sum(-m[i]*(-`r'`[i, c]^2*Physics[KroneckerDelta][a, b]+`r'`[i, a]*`r'`[i, b]), i = 1 .. n))*(diff(Omega[b](t), t)) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a))

`𝕀`[a, b]*(diff(Omega[b](t), t)) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(636)

From the result (614) for the problem of representing the angular momentum of a rigid body in terms of the inertia tensor  the left-hand side is the derivative of the components of the angular momentum

`L'`[a] = `𝕀`[a, b]*Omega[b]

`L'`[a] = `𝕀`[a, b]*Omega[b]

(637)

subs(Omega[b] = Omega[b](t), (diff(L(x), x))[a] = (diff(L(x), x))[a](t), `L'`[a] = `𝕀`[a, b]*Omega[b])

`L'`[a](t) = `𝕀`[a, b]*Omega[b](t)

(638)

diff(`L'`[a](t) = `𝕀`[a, b]*Omega[b](t), t)

diff(`L'`[a](t), t) = `𝕀`[a, b]*(diff(Omega[b](t), t))

(639)

subs((rhs = lhs)(diff(`L'`[a](t), t) = `𝕀`[a, b]*(diff(Omega[b](t), t))), `𝕀`[a, b]*(diff(Omega[b](t), t)) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a))

diff(`L'`[a](t), t) = Physics:-Vectors:-Component((D[2](U))(R_(t), Phi_(t)), a)

(640)

The right-hand side is the a^th component of the variation of the potential energy U(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`) with respect to `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`. A graphics analysis of (D[2](U))(`#mover(mi("R"),mo("→"))`, `#mover(mi("Φ",fontstyle = "normal"),mo("→"))`) and algebraic derivation as done in the problem of the section Angular velocity  results in

(D[2](U))(R_(t), Phi_(t)) = `&x`(R_(t), F_[ext])

(D[2](U))(R_(t), Phi_(t)) = Physics:-Vectors:-`&x`(R_(t), F_[ext])

(641)

subs((D[2](U))(R_(t), Phi_(t)) = Physics[Vectors][`&x`](R_(t), F_[ext]), diff(`L'`[a](t), t) = Physics[Vectors][Component]((D[2](U))(R_(t), Phi_(t)), a))

diff(`L'`[a](t), t) = Physics:-Vectors:-Component(Physics:-Vectors:-`&x`(R_(t), F_[ext]), a)

(642)

As shown in this problem  of the section on the equations of motion for many-particle systems , the right-hand side of this result of that is the total torque `#mover(mi("N"),mo("→"))`(t), resulting in the second equation of movement of a rigid body, for the time derivative of the angular momentum

subs(rhs(diff(`L'`[a](t), t) = Physics[Vectors][Component](Physics[Vectors][`&x`](R_(t), F_[ext]), a)) = N_[a](t), diff(`L'`[a](t), t) = Physics[Vectors][Component](Physics[Vectors][`&x`](R_(t), F_[ext]), a))

diff(`L'`[a](t), t) = N_[a](t)

(643)

NULL

Non-inertial coordinate systems

 

 

When describing the motion of a particle as seen from a non-inertial reference system (e.g. a rotating planet, like the Earth), we also see "acceleration" that is not due to any force but instead to the fact that the reference system itself is accelerated.

 

Problem
Consider a non-inertial reference system K' which moves with non-constant translational velocity `#mover(mi("V"),mo("→"))`(t) with regards to an inertial reference system K__0. 

a) Show that the Lagrangian L' in K' is given by

diff(L(x), x) = (1/2)*m*`#mover(mi("v'"),mo("→"))`-m*`#mover(mi("W"),mo("→"))`.`#mover(mi("r'"),mo("→"))`-U

where `#mover(mi("W"),mo("→"))` = diff(`#mover(mi("V"),mo("→"))`(t), t) is the translational acceleration of the frame K'  as seen from K__0.

 

b) Show that the Lagrange equation derived from this Lagrangian in the frame K'  is

m*(diff(`#mover(mi("v'"),mo("→"))`(t), t)) = -Nabla(U)-m*`#mover(mi("W"),mo("→"))`

Solution

   

 

Coriolis force and centripetal force

 

 

Problem

Consider a second non-inertial frame of reference J whose origin coincides with that of K', but which rotates relative to K' with variable angular velocity `#mover(mi("Ω",fontstyle = "normal"),mo("→"))`(t). Denote the position vector and velocity in the non-inertial frame J as `#mover(mi("r"),mo("→"))` and `#mover(mi("v"),mo("→"))`,

a) Show that the Lagrangian L in the non-inertial frame J is given by

L = (1/2)*m*`#mover(mi("v"),mo("→"))`^2+m*`#mover(mi("v"),mo("→"))`.`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `#mover(mi("r"),mo("→"))`)+(1/2)*m*`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `#mover(mi("r"),mo("→"))`)^2-m*(`#mover(mi("W"),mo("→"))`.`#mover(mi("r"),mo("→"))`)-U

b) Show that the Lagrange equation derived from this Lagrangian in the frame J  is

m*(diff(`#mover(mi("v"),mo("→"))`(t), t)) = -Nabla(U)-m*`#mover(mi("W"),mo("→"))`+m*`&x`(`#mover(mi("r"),mo("→"))`, diff(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`(t), t))+2*m*`&x`(`#mover(mi("v"),mo("→"))`, `#mover(mi("Ω",fontstyle = "normal"),mo("→"))`)+m*`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("Ω",fontstyle = "normal"),mo("→"))`))

where 2*m*`&x`(`#mover(mi("v"),mo("→"))`, `#mover(mi("Ω",fontstyle = "normal"),mo("→"))`) is the Coriolis force and m*`&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `&x`(`#mover(mi("r"),mo("→"))`, `#mover(mi("Ω",fontstyle = "normal"),mo("→"))`)) is the centrifugal force.

 

Solution

 

 

restart; with(Physics:-Vectors)

 

a) The starting point is the Lagrangian in the frame K__0. Denoting vectors and the Lagrangian in K__0 with the suffix 0, L__0 is given by

CompactDisplay((r__0_, v__0_, (diff(r(x), x))*_, (diff(v(x), x))*_, r_, v_, V_, W_, Omega_)(t))

Omega_(t)*`will now be displayed as`*Omega_

(658)

L__0 = (1/2)*m*v__0_(t)^2-U(r__0_(t))

L__0 = (1/2)*m*Physics:-`^`(v__0_(t), 2)-U(r__0_(t))

(659)

The velocities of the particle in the frames K__0 and K' are related by

v__0_(t) = (diff(v(x), x))*_(t)+V_(t)

v__0_(t) = `v'_`(t)+V_(t)

(660)

where `#mover(mi("V"),mo("→"))` is the translational velocity of K' viewed from K__0. Inserting this relation (660) into L__0 gives diff(L(x), x), the result of the previous problem

`L'` = (1/2)*m*Physics[Vectors][Norm](`v'_`(t))^2-m*Physics[Vectors][`.`](`r'_`(t), W_(t))-U(`r'_`(t))

`L'` = (1/2)*m*Physics:-Vectors:-Norm(`v'_`(t))^2-m*Physics:-Vectors:-`.`(`r'_`(t), W_(t))-U(`r'_`(t))

(661)

In turn, the velocities `#mover(mi("v'"),mo("→"))` and `#mover(mi("v"),mo("→"))` in the frames K' and J are related by

(diff(v(x), x))*_(t) = v_(t)+`&x`(Omega_(t), r_(t))

`v'_`(t) = v_(t)+Physics:-Vectors:-`&x`(Omega_(t), r_(t))

(662)

where `#mover(mi("Ω",fontstyle = "normal"),mo("→"))` is the angular velocity of the frame J viewed from K'. Also, since K' and J have the same origin, `#mover(mi("r'"),mo("→"))` = `#mover(mi("r"),mo("→"))` and the Lagrangian L in J  is

subs(`v'_`(t) = v_(t)+Physics[Vectors][`&x`](Omega_(t), r_(t)), (diff(r(x), x))*_ = r_, diff(L(x), x) = L(x), `L'` = (1/2)*m*Physics[Vectors][Norm](`v'_`(t))^2-m*Physics[Vectors][`.`](`r'_`(t), W_(t))-U(`r'_`(t)))

L = (1/2)*m*Physics:-Vectors:-Norm(v_(t)+Physics:-Vectors:-`&x`(Omega_(t), r_(t)))^2-m*Physics:-Vectors:-`.`(r_(t), W_(t))-U(r_(t))

(663)

expand(L = (1/2)*m*Physics[Vectors][Norm](v_(t)+Physics[Vectors][`&x`](Omega_(t), r_(t)))^2-m*Physics[Vectors][`.`](r_(t), W_(t))-U(r_(t)))

L = (1/2)*m*Physics:-Vectors:-Norm(v_(t))^2+m*Physics:-Vectors:-`.`(v_(t), Physics:-Vectors:-`&x`(Omega_(t), r_(t)))+(1/2)*m*Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(r_(t))^2-(1/2)*m*Physics:-Vectors:-`.`(Omega_(t), r_(t))^2-m*Physics:-Vectors:-`.`(r_(t), W_(t))-U(r_(t))

(664)

Two of these terms can be regrouped as `&x`(`#mover(mi("Ω",fontstyle = "normal"),mo("→"))`, `#mover(mi("r"),mo("→"))`)^2

`&x`(Omega_(t), r_(t))^2

Physics:-`^`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), 2)

(665)

expand(Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)) = Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)

Physics:-Vectors:-Norm(Omega_(t))^2*Physics:-Vectors:-Norm(r_(t))^2-Physics:-Vectors:-`.`(Omega_(t), r_(t))^2 = Physics:-`^`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), 2)

(666)

simplify(L = (1/2)*m*Physics[Vectors][Norm](v_(t))^2+m*Physics[Vectors][`.`](v_(t), Physics[Vectors][`&x`](Omega_(t), r_(t)))+(1/2)*m*Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](r_(t))^2-(1/2)*m*Physics[Vectors][`.`](Omega_(t), r_(t))^2-m*Physics[Vectors][`.`](r_(t), W_(t))-U(r_(t)), {Physics[Vectors][Norm](Omega_(t))^2*Physics[Vectors][Norm](r_(t))^2-Physics[Vectors][`.`](Omega_(t), r_(t))^2 = Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)})

L = (1/2)*m*Physics:-Vectors:-Norm(v_(t))^2+m*Physics:-Vectors:-`.`(v_(t), Physics:-Vectors:-`&x`(Omega_(t), r_(t)))-m*Physics:-Vectors:-`.`(r_(t), W_(t))+(1/2)*m*Physics:-`^`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), 2)-U(r_(t))

(667)

which is the expected result.

_______________________________________

 

b) To compute Lagrange's equation in vectorial form, one can use the standard formula as in the previous problem

 

%diff(%diff(L = (1/2)*m*Physics[Vectors][Norm](v_(t))^2+m*Physics[Vectors][`.`](v_(t), Physics[Vectors][`&x`](Omega_(t), r_(t)))-m*Physics[Vectors][`.`](r_(t), W_(t))+(1/2)*m*Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)-U(r_(t)), v_(t)), t)-%diff(L = (1/2)*m*Physics[Vectors][Norm](v_(t))^2+m*Physics[Vectors][`.`](v_(t), Physics[Vectors][`&x`](Omega_(t), r_(t)))-m*Physics[Vectors][`.`](r_(t), W_(t))+(1/2)*m*Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)-U(r_(t)), r_(t))

%diff(%diff(L = (1/2)*m*Physics:-Vectors:-Norm(v_(t))^2+m*Physics:-Vectors:-`.`(v_(t), Physics:-Vectors:-`&x`(Omega_(t), r_(t)))-m*Physics:-Vectors:-`.`(r_(t), W_(t))+(1/2)*m*Physics:-`^`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), 2)-U(r_(t)), v_(t)), t)-%diff(L = (1/2)*m*Physics:-Vectors:-Norm(v_(t))^2+m*Physics:-Vectors:-`.`(v_(t), Physics:-Vectors:-`&x`(Omega_(t), r_(t)))-m*Physics:-Vectors:-`.`(r_(t), W_(t))+(1/2)*m*Physics:-`^`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), 2)-U(r_(t)), r_(t))

(668)

 

Evaluate these derivatives and replace diff(`#mover(mi("r"),mo("→"))`(t), t) = `#mover(mi("v"),mo("→"))` 

subs(diff(r_(t), t) = v_(t), value(%diff(%diff(L = (1/2)*m*Physics[Vectors][Norm](v_(t))^2+m*Physics[Vectors][`.`](v_(t), Physics[Vectors][`&x`](Omega_(t), r_(t)))-m*Physics[Vectors][`.`](r_(t), W_(t))+(1/2)*m*Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)-U(r_(t)), v_(t)), t)-%diff(L = (1/2)*m*Physics[Vectors][Norm](v_(t))^2+m*Physics[Vectors][`.`](v_(t), Physics[Vectors][`&x`](Omega_(t), r_(t)))-m*Physics[Vectors][`.`](r_(t), W_(t))+(1/2)*m*Physics[`^`](Physics[Vectors][`&x`](Omega_(t), r_(t)), 2)-U(r_(t)), r_(t))))

0 = m*(diff(v_(t), t))+m*(Physics:-Vectors:-`&x`(diff(Omega_(t), t), r_(t))+Physics:-Vectors:-`&x`(Omega_(t), v_(t)))-m*Physics:-Vectors:-`&x`(v_(t), Omega_(t))+m*W_(t)-m*Physics:-Vectors:-`&x`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), Omega_(t))+(D(U))(r_(t))

(669)

Isolating m*(diff(`#mover(mi("v"),mo("→"))`(t), t))  and collecting vector products we get the expected result

isolate(0 = m*(diff(v_(t), t))+m*(Physics[Vectors][`&x`](diff(Omega_(t), t), r_(t))+Physics[Vectors][`&x`](Omega_(t), v_(t)))-m*Physics[Vectors][`&x`](v_(t), Omega_(t))+m*W_(t)-m*Physics[Vectors][`&x`](Physics[Vectors][`&x`](Omega_(t), r_(t)), Omega_(t))+(D(U))(r_(t)), m*(diff(v_(t), t)))

m*(diff(v_(t), t)) = -m*(Physics:-Vectors:-`&x`(diff(Omega_(t), t), r_(t))-Physics:-Vectors:-`&x`(v_(t), Omega_(t)))+m*Physics:-Vectors:-`&x`(v_(t), Omega_(t))-m*W_(t)+m*Physics:-Vectors:-`&x`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), Omega_(t))-(D(U))(r_(t))

(670)

collect(m*(diff(v_(t), t)) = -m*(Physics[Vectors][`&x`](diff(Omega_(t), t), r_(t))-Physics[Vectors][`&x`](v_(t), Omega_(t)))+m*Physics[Vectors][`&x`](v_(t), Omega_(t))-m*W_(t)+m*Physics[Vectors][`&x`](Physics[Vectors][`&x`](Omega_(t), r_(t)), Omega_(t))-(D(U))(r_(t)), `&x`)

m*(diff(v_(t), t)) = m*Physics:-Vectors:-`&x`(Physics:-Vectors:-`&x`(Omega_(t), r_(t)), Omega_(t))-m*Physics:-Vectors:-`&x`(diff(Omega_(t), t), r_(t))+2*m*Physics:-Vectors:-`&x`(v_(t), Omega_(t))-m*W_(t)-(D(U))(r_(t))

(671)

NULL

 

 

Part II (forthcoming)

The Hamiltonian and equations of motion; Poisson brackets

   

Canonical transformations

   

The Hamilton-Jacobi equation

   

References

   

 

NULL

Download Mechanics_(2022).mw

Download Mechanics_(2022).pdf
 

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

This is about functionality introduced in Maple 2022, which however is still not well known: Integral Vector Calculus and parametrization using symbolic (algebraic) vector notation. Four new commands were added to the Physics:-Vectors package, implementing the parametrization of curves, surfaces and volumes, as well as the computation of path, surface and volume vector integrals. Those are integrals where the integrand is a scalar or vector function. The computation is done from any description (algebraic, parametric, vectorial) of the region of integration - a path, surface or volume.
 
There are three kinds of line or path integrals:

NOTE Jan 1: Updated the worksheet linked below; it runs in Maple 2022.
Download Integral_Vector_Calculus_and_Parametrization.mw

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

This command should have been in Physics on day one. Being more familiar with functional differentiation, and Physics:-Fundiff was the first Physics command that ever existed, I postponed writing LagrangeEquations year after year. In general, however, functional differentiation is seen as a more advanced topic. So there is now a new command, Physics:-LagrangeEquations, taking advantage of functional differentiation on background, and distributed for everybody using Maple 2022.2 within the Maplesoft Physics Updates. This is the first version of its help page.


Download LagrangeEquations.mw

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

Welcome back to another Maple Learn blog post! Today we’re going to talk about the gift-wrapping algorithm, used to find the convex hull of a set of points. If you’re not sure what that means yet, don’t worry! We’re going to go through it with four Maple Learn documents; two which are background information on the topic, one that is a visualization for the gift-wrapping algorithm, and another that goes through the steps. Each will be under their own heading, so feel free to skip ahead to your skill level!

Before we can get into the gift-wrapping algorithm we need to define a few terms. Let’s start by defining polygons and simple polygons.

A picture containing text, person

Description automatically generated

Polygon: A closed shape created by joining a series of line segments.

Simple polygon: A polygon without holes and that does not intersect itself.

Shape, polygon

Description automatically generated

So, what are convex and concave polygons? Well, there are three criteria that define a convex polygon. A polygon that is not convex is called concave. The criteria are…

  1. Any line segment connecting any two points within the polygon stays within the polygon.
  2. Any line intersects a polygon’s boundary at most twice.
  3. All interior angles are less than 180 degrees or pi radians.

A picture containing chart

Description automatically generated

Because the criteria are equivalent, if any one is missing, the shape is concave. AKA, all three criteria must be present for a shape to be convex. Most “regular shapes”, such as trapezoids, are convex polygons!

A shape that satisfies convex criteria but not the criteria for being a polygon is called a convex set.

As mentioned at the start of this post, the gift-wrapping algorithm is used to find the convex hull of a set of points. Now that we know what convex polygons and convex sets are, we can define the convex hull!

Convex hull: The convex set of a shape or several shapes that fully contains the object and has the smallest possible area.

A picture containing text, stationary, envelope, businesscard

Description automatically generated

Why was the convex polygon important? Well, the convex hull of a set of points is always a convex polygon. Some of the points in the set are the vertices of said polygon, and are called extreme points. You can find the convex hull of either concave or convex polygons.

This document amazed me when I tried it for the first time. Here, you can generate a set of points with the “Generate Another” button, and then press the “Visualize” button. The document then calculates the perimeter of the convex hull of the set of points! The set can be further customized below the buttons, by changing the number of points. The other option below it allows you to slow down or speed up the visualization. Pretty cool, huh? It’s like it’s thinking!

Try the document out a few times, or watch the gif below to get a quick idea of it.

This final document walks you through the steps of how to use the gift wrapping algorithm. It is a simple loop of 4 steps, with one set-up step. Unlike the other documents in this post, I won’t be delving too far into the math behind the steps. I want to encourage you to check this one out yourself, as it’s really quite a fun problem to solve once you have some time!

Chart, line chart

Description automatically generated

I hope you check out the documents in this post. Please let us know below if there’s any other documents you’d like to see featured!

A failing slinky is another intriguing physics phenome that can be easily reproduced with MapleSim.

The bottom of a vertically suspended slinky does not move when the top is released until the slinky is fully collapsed.

 

 

To model this realistically in MapleSim, it is necessary to

  • Establish a stretched equilibrium state at the start of the fall
  • Avoid penetration of windings when windings collapse (i.e. get into contact)

The equilibrium state is achieved with the snapshot option. Penetration is avoided with the Elasto Gap component. Details can be found in the attached model.

A good overview of “Slinky research” is given here. The paper provides a continuous description of the collapse process (using an inhomogenous wave equation combined with contact modeling!!!) and introduces a finite time for the collapse of all windings. Results for a slinky are presented that collapses after 0.27s. The attached model has sufficient fidelity to collapse at the same time.

Real Slinkies also feature a torsional wave that precedes the compression wave and disturbs an ideal collapse. This can be seen on slowmo footage and advanced computer models. With a torsion spring constant at hand (are there formulas for coil springs?), it could also be modeled with MapleSim.

Falling_slinky.msim

Have you ever heard of the Maurer Rose?

The Maurer Rose was demonstrated in 1987 by Peter Maurer and is created by connecting certain points on a rose curve. This creates petal-like patterns, caused by the oscillation of a sine curve.

 

Chart, radar chart

Description automatically generated

So, how are these created? A "rose curve" is created in polar coordinates with the equation sin(nt) for a (positive integer) value of n.  To create the Maurer Rose, straight line segments are drawn connecting points on the curve at incrementing angle values.  The size of this increment (called d in our examples) leads to different patterns of lines across the curve.

This can be done in Maple Learn! One example of the Maurer Rose already exists, complete with a full interactive visualization and a more detailed overview of the Maurer Rose.

Play around with it and look below at some of the different shapes that can be created using this document! The first is created with an n value of 31 and a d value of 65, with blue and red. The second uses an n value of 4 and a d value of 133, and purple and green.

Chart, diagram, radar chart

Description automatically generated

Are there any other concepts you’d like to see represented in Maple Learn’s document gallery? Please let us know in the comments below!

Hello MaplePrimes community,

We just created a Frequently Asked Question article that may address some Primes questions about updates to Physics in Maple 2022.2: 

Why does Maple 2022.2 throw an error executing Physics:-Version(latest)?

For searchability, the specific error in question is

Error, (in Physics:-Version) unable to determine the Physics Updates version, could you please report the problem to support [at] maplesoft [dot] com

 

  • Maplesoft will work to improve package updating in future versions of Maple.
  • In Maple 2022.2, the workaround is to install and/or update the Maplesoft Physics Updates using the MapleCloud toolbar.

 

Have you heard of Maple Scripting before? Do you want to extend your Maple Learn documents with your Maple knowledge? Scripting is the process of using Maple to create Maple Learn documents. If you’re already used to Maple, this may be a piece of cake for you, but we wanted to start from the basics for anyone who wants to extend their Maple Learn and Maple knowledge. This process can be used for many different types of documents, from quizzes to intensive 3D visualizations.

So, let’s get started! All Maple Learn document scripting needs the DocumentTools:-Canvas package. The canvas, as you know, is that white space in a Maple Learn document. Therefore, this package is the core content of a scripted document! Always put:

with(DocumentTools:-Canvas):

At the top of your code, or put

uses DocumentTools:-Canvas:

At the start of your procedures.

Now that we’ve told Maple to use the DocumentTools:-Canvas, we need to create a canvas.

Canvases are created as variables, using the command NewCanvas. Inside NewCanvas, you will add a square-bracket list of all the content you want to see inside. For now, just know that you can add text cells with Text(“YOUR TEXT”) and a math cell with Math(YOUR MATH). On the next line, make sure to put either ShareCanvas(YOUR CANVAS VARIABLE) or ShowCanvas(YOUR CANVAS VARIABLE).  ShareCanvas creates a Maple Learn sharelink, while ShowCanvas shows the canvas directly in Maple. Note that ShowCanvas does not have every Maple Learn feature, but makes quick work of fast error checking.

canvas := NewCanvas([Text(“My first canvas”), Math(3*x+2*y)]):

ShareCanvas(canvas);

There are two more things I want to show you in this post: How to make a group have multiple cells (instead of just the one), and how to position your items on the canvas. Let’s start with group making.

To create a group with multiple cells, use the Group() command within the NewCanvas command, and separate the cells with commas, in a list. You don’t need to specify Text() or Math() when using Group().

canvas := NewCanvas([Group([“This is the first cell…”, “The second….”, “and the third.”])]):

At the end of any command/canvas element, within the brackets, you can define position=[x,y] to specify where on the canvas the object should go. You can adjust the precision pixel by pixel until you are happy with the layout.

When we put all these together, we get code that looks like this:

with(DocumentTools:-Canvas):

canvas := NewCanvas([

Group(["This is the first cell…", "The second…", "and the third."], position=[200,200]),

Math(3*x+2, position=[100,100]),

Text("This is text!", position=[400,400])]):

ShareCanvas(canvas);

And in the end, your scripted document looks like this.

We hope this helps you get started with Maple Scripting. There will be another post on even more of what we can do with Maple Scripting, and how we can make these documents even more interactive. Let us know if there’s anything specific you want to see in that post!

 

Greetings, fellow educators, researchers, engineers, students, and folx who love mathematics! 

 

I believe in the importance of mathematics as a structure to our society, as a gateway to better financial decision making, and as a crucial subject to teach problem solving. I also believe in the success of all students, through self-discovery and creativity, while working with others to create their own knowledge. Consequently, I’ve designed my examples in the Maple Learn gallery to suit these needs. Many of my documents are meant to be “stand-alone” investigations, summary pages, or real-world applications of mathematical concepts meant to captivate the interest of students in using mathematics beyond the basic textbook work most curricula entail. Thus, I believe in the reciprocal teaching and learning relationship, through the independence and creativity that technology has afforded us. The following is an example of roller coaster track creation using functions. Split into a five part investigation, students are tasked to design the next roller coaster in a theme park, while keeping in mind the elements of safety, feasibility, and of course fun!

Common elements we take for granted such as having a starting and ending platform that is the same height (since most coasters begin and end at the same location), boarding the coaster on a flat surface, and smooth connections between curves translate into modeling with functions. 

 

Aside from interning with Maplesoft, I am an educator, researcher, student, financial educator, and above all, someone who just loves mathematics and wishes to share that joy with the whole world. As a practicing secondary mathematics and science teacher in Ontario, Canada, I have the privilege of taking what I learned in my doctorate studies and applying it to my classrooms on a daily basis. I gave this assignment to my students and they really enjoyed creating their coasters as it finally gave them a reason to learn why transformations of quadratics, amongst other functions, were important to learn, and where a “real life” application of a piecewise function could be used. 

Graphical user interface, text, application, Word

Description automatically generated

 

Having worked with the Ontario and International Baccalaureate mathematics curricula for over a decade, I have seen its evolution over time and in particular, what concepts students struggled to understand, and apply them to the “real world.” Concurrently, working with international mathematics curricula as part of my collaboration with Maplesoft, I have also seen trends and emergent patterns as many countries’ curricula have evolved to incorporate more mathematical literacy along with competencies and skills. In my future posts, you will see Maple Learn examples on financial literacy since working as a financial educator has allowed me to see just how ill prepared families are towards their retirement and how we can get lost amongst a plethora of options provided by mass media. Hence, I have 2 main goals I dedicate to a lifelong learning experience; financial literacy and greater comprehension of mathematics topics in the classroom. 

 

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