Posts - MaplePrimes

Physics Courseware Support: Mechanics

by: ecterrab 14540 Product: Maple

February 09 2023

12 1

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.

Part I

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

a.	The position as a function of time

b.	The velocity

c.	The acceleration

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

[`&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 vectorsis given by

(2)

Problem

Rewrite the position vector 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 do not depend on time, they are constant vectors. So is entered as

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

(21)

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

(22)

(23)

(24)

The position 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

Problem

Rewrite the velocity in cylindrical and spherical coordinates while making the curvilinear unit vectors' time dependency explicit .

	Solution

The acceleration

Problem

Rewrite the acceleration in cylindrical and spherical components while making the curvilinear unit vectors' time dependency explicit.

	Solution

Deriving these formulas

All these results for the position , velocity and acceleration 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)

(38)

The above result contains, in the last equation, the cylindrical radial unit vector ; rewrite it in the spherical basis

(39)

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

(40)

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

(41)

(42)

Introducing the differentiation rules (40) for the unit vectors

(43)

Performing the inert (grayed) derivatives

(44)

In the same way, for the acceleration

(45)

%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)

%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)

(48)

(49)

(50)

Collect vector components

(51)

Summary

•	You can express and 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 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))

(62)

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

(63)

where is the acceleration and is the linear momentum, so in terms of

(64)

We define the angular momentum of a particle, and the torque acting upon it, as

(65)

(66)

Differentiating the definition of

(67)

Since is parallel to , the first term in the above cancels, and in the second term, from (64),

(68)

from which

(69)

•	As discussed below , in the case of a closed system, and these two equations result in

that is, the linear and angular momentum are conserved quantities. Note that does not require that , only that .

The equations of motion - vectorial form

Problem

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

a) The trajectory starting from
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 be the position of the particle in a reference system; then, the velocity and acceleration are given by

(70)

(71)

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

(72)

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

_______________________________________

b) The integration of vectorial equations is also frequently performed after expressing , and in a particular system of coordinates. For example, in the Cartesian system (71) has the form

(73)

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

(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

_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 and , specifying the integration constants ${c__3, c__4, c__5, c__6, c__7, c__8}$ , can also be written in components

(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:

${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)

The case of constant acceleration

Problem

Starting from the vectorial equation (72) for , 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 and 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 , where 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 , 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 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

a) Determine the regions where the motion has positive and negative acceleration. Compute the position at .

b) Compute the velocity corresponding to , starting - not from this expression for but from the acceleration

	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.

	Cartesian coordinates

	Curvilinear coordinates

Many-particle systems

Center of Mass

Given a system of n particles of masses with positions 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 measured in K. Show that:

a) When is the velocity of the center of mass, the total momentum measured in K' is equal to 0.

b) The relation between and the velocity of the center of mass, both measured in K, is the same as the relation 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 2^nd law for each particle implies that , where is the center of mass and 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 satisfies

b) The total torque 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 with the horizontal plane and its projection BC over this plane forms an angle 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))

(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 and derive the equations of motion as the Lagrange equations for L.

For closed systems, the potential energy is related to the force acting on each particle by the equation . Formally, is the gradient operator in the basis onto which is projected, and with respect to its coordinates - say in Cartesian basis and coordinates .

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

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

(209)

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

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

(210)

where is the velocity of the i^th particle. Generally speaking, the potential energy of the system is a function of the coordinates of the n particles, and the Lagrangian is defined as

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

(211)

The potential energy is related to the force acting on each particle by the equation . Formally, is the gradient operator in the basis onto which is projected. Knowing the Lagrangian, we can derive the (Lagrange) equations of motion as

(212)

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

b) Oscillates horizontally on the plane of the pendulum according to .

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

	Solution

Conservation laws

Work

By definition, the work performed by a force to move a particle an infinitesimal amount is . Thus, the work to move it from to along some path C is

"((∫)[A]^(B) (F ) . &DifferentialD;r)[path=C]"

Problem

A particle is submitted to a force whose Cartesian components are given by. Calculate the work done by this force when moving the particle along a straight line from the origin to a point ().

	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 i^th particle of the system can be derived from a potential, . 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 . 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 2^nd law follows .

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 3^rdlaw.

Problem

Consider a system of n particles measured in two frames of reference K and K' that move relative to each other with velocity . Show that the system's momenta and are related by

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

	Solution

Problem

A particle of mass m moving with velocity leaves a half-space in which its potential energy is a constant and enters another in which its potential energy is a different constant . Determine the change in direction of motion of the particle; that is, sin(`θ__1`)/sin(`θ__2`) where and are the angles between an axis perpendicular to the separating plane and the momentum 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 . From it follows that , that is, the angular momentum is conserved.

Problem

a) Express the Cartesian components of the angular momentum , as well as its norm, in cylindrical and spherical coordinates.

b) Rewrite 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 from each other. Show that the momenta and of the system are related by

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 . Show that the momenta and 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 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 and are the total linear and angular momenta measured in K'.
b) Show that when the origin of K' is the center of mass , this formula reduces to

,

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

Solution

a) The momentum of the system measured in K is given by

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 and using , where :

(323)

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

(324)

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 (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 and , introduce the relation between and

(327)

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)

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 and (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 , then

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

(331)

and were related in `#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 in terms of instead of

$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)

(333)

which is the result we were looking for.

Cyclic coordinates

Any generalized coordinate 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 is cyclic, the right-hand side is 0 and so the quantity is conserved.

Problem

The Lagrangian describing the movement of a particle in a central field has as a cyclic coordinate. Using cylindrical coordinates, show that the corresponding conserved quantity is the z component of the angular momentum

	Solution

Integration of the equations of motion

Motion in one dimension

Problem
For a closed system, or any system where the total energy 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 .

	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 .

	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 the result of a) becomes the result obtained for the previous two-body problem, equation .

Solution

a) Let represent the position vector of the particle of mass M and those of the particles of mass m. The Lagrangian is given by

(386)

(387)

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

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

(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

(391)

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)≡ is sufficient to eliminate and from the Lagrangian (389). In steps, eliminating 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, {`&Ropf;_`[a](t) = r_[a](t)-r__M_(t)}, {r_[a](t)})$

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

(393)

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

(394)

To eliminate , use

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

(395)

Likewise, to eliminate use

(396)

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

L = (1/2)*M*Physics:-Vectors:-Norm(diff(-m*(Sum(`&Ropf;_`[a](t), a = 1 .. n))/(m*n+M), t))^2+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`&Ropf;_`[a](t)-m*(Sum(`&Ropf;_`[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

L = (1/2)*M*m^2*Physics:-Vectors:-Norm(Sum(diff(`&Ropf;_`[a](t), t), a = 1 .. n))^2/(m*n+M)^2+(1/2)*m*(Sum(Physics:-Vectors:-Norm(diff(`&Ropf;_`[a](t), t))^2-2*m*Physics:-Vectors:-`.`(diff(`&Ropf;_`[a](t), t), Sum(diff(`&Ropf;_`[a](t), t), a = 1 .. n))/(m*n+M)+m^2*Physics:-Vectors:-Norm(Sum(diff(`&Ropf;_`[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

(399)

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

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

(400)

This Lagrangian involves only the n relative position vectors , achieving the reduction of the original problem of particles to a problem of only n particles.

_______________________________________

b) To obtain the result for the two-body problem, substitute into

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

(401)

Expand only the Sums and collect the Norm

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

(402)

L = (1/2)*M*m*Physics:-Vectors:-Norm(diff(`&Ropf;_`[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

(404)

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

(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(`&Ropf;_`[1](t), t))^2/(M+m)-U), {m[1]*m[2]/(m[1]+m[2]) = mu})$

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

(406)

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 of a particle that moves in a central field is conserved, so evolves in time on a fixed plane perpendicular to the constant . Show that the surface swept per second by the position vector is constant (Kepler's second law).

	Solution

Problem
Starting from the constancy of the energy E and the angular momentum , 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 and

b) interactively, one step at a time, uncouple the equations of movement eliminating from the problem, resulting in an implicit solution .

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

	Solution

Kepler's problem

Problem

Show that, when , with 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

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, . 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 . Write the equation of motion as Newton's 2^nd 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 acts. For the oscillations to be small, must produce only small displacements. The total force is .

a) Write the equation of motion as Newton's 2^nd 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 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 .

d) Show that a solution for the case considered in c), that is, , 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 .

e) Specialize the solution of item d) in the case of resonance, when , 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 and 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

where represents the displacement of the generalized coordinate , 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))

(524)

Denoting the generalized coordinates by , the potential energy can be expanded in series around the minimums . Since the movement consist only in small displacements around , it is sufficient to keep terms in the expansion up to order 2, resulting in an expression analogous to 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))

(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))

(526)

where we take the at the minimums and denote them as

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 and can be split into symmetric and antisymmetric parts, with the antisymmetric parts canceling out in view of the symmetric character of in U and in T. Therefore, we can take and 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 describe independent displacements, so one can think of 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.

(528)

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

${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 can now be written as

(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

(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, Physics:-Vectors:-diff(x[c](t), t)), t) = %diff(L, x[c](t))

(532)

(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,

(534)

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 of each particle by , where is the mass density as a function of the position and 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 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 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 of a point P of a body, measured in an inertial reference system K, can be written as

where is the velocity of the origin of K', a frame of reference attached to the body's center of mass, is the body's angular velocity (its instantaneous counter-clockwise rotation speed around some axis in the direction of a unit vector) and 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 and are related by a rotation matrix which, as all rotation matrices, is orthogonal.

	Solution

Inertia tensor

Problem

a) Show that using , derived in Angular velocity for the velocity of a point P of a rigid body in terms of and the position of P viewed from the center of mass , the kinetic energy of the rigid body can be written in terms of the positions of the n particles (not their velocities), the velocity of the center of mass and angular velocity 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

where

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

is the Inertia tensor, represents the components of the vector and represents the components of the position vector of the particle.

	Solution

Problem

Determine the Inertia tensor corresponding to a triatomic molecule that has the form of an isosceles triangle with two masses in the extremes of the base and a mass in the third vertex. The distance between the two masses 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 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 of the angular momentum of a rigid body can be expressed in terms of the inertia tensor and the components of the angular velocity as

	Solution

The equations of motion of a rigid body

A rigid body is a system with six degrees of freedom: three indicating the position 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

(615)

(616)

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

(617)

where is the total momentum, is the total external force acting upon the body, is the external force acting upon the particle, the total angular momentum and is the total torque.

Problem
Show that the equations of movement of a rigid body can be computed as the Lagrange equations for and from the Lagrangian

where mu = Sum(m[i], i = 1 .. n) is the total mass, is the inertia tensor, is the potential energy for the external force , and .

Solution

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

(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

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 (location of the center of mass) and (the rotation axis, so that ), the Lagrangian in vectorial form is given by

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 is derived as the Lagrange equation for this Lagrangian

(621)

(622)

(623)

On the right-hand side, with respect to its first argument , equivalent to the gradient taking as the coordinates, and is the total momentum

(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 . If more precision is required, the dependency of the potential energy could be expressed in terms of the norm of , as in

(625)

in which case differentiating with respect to is equivalent to the definition of directional derivative

(626)

The same computation, this time with respect to the coordinates where is the corresponding velocity,

(627)

(628)

(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 of the angular momentum .

$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)

(631)

(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)

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)

(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,

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

(635)

(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

(637)

(638)

(639)

(640)

The right-hand side is the component of the variation of the potential energy with respect to . A graphics analysis of and algebraic derivation as done in the problem of the section Angular velocity results in

(641)

(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 , resulting in the second equation of movement of a rigid body, for the time derivative of the angular momentum

(643)

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 with regards to an inertial reference system .

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 is the translational acceleration of the frame K' as seen from .

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

	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 . Denote the position vector and velocity in the non-inertial frame J as and ,

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

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

where is the Coriolis force and is the centrifugal force.

Solution

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

(658)

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

(659)

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

(660)

where is the translational velocity of K' viewed from . Inserting this relation (660) into gives , 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))

(661)

In turn, the velocities and in the frames K' and J are related by

(662)

where is the angular velocity of the frame J viewed from K'. Also, since K' and J have the same origin, and the Lagrangian in J is

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)

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

(665)

(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

(668)

Evaluate these derivatives and replace

(669)

Isolating and collecting vector products we get the expected result

(670)

(671)

Part II (forthcoming)

	The Hamiltonian and equations of motion; Poisson brackets

	Canonical transformations

	The Hamilton-Jacobi equation

	References

Download Mechanics_(2022).mw

Download Mechanics_(2022).pdf

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

Meet the Co-ops! - Winter 2023

by: jthomson 90

February 03 2023

5 0

Have you ever wondered who the students are that help create Maplesoft’s family of products?

In this blog, we thought that it was fitting to introduce ourselves and give the MaplePrimes community some insight into the students who are committed to helping Maplesoft improve its products and who believe that Math Matters!

I’ll begin. My name is Jack Thomson and I’m in my second year of the Mathematics (Waterloo) and Business Administration (Laurier) Double Degree. This term I am the Product Management Co-op at Maplesoft where I will be helping support the development of Maplesoft's academic market products, including Maple Learn and Maple Calculator. My favorite areas of math are statistics and probability. These areas are my favorite since I like to be able to draw conclusions from data and predict the future with past trends. I am also fascinated with probability since it allows us to make more educated decisions about real-life events. This ties into my belief of why Math Matters, since it is hidden in every aspect of life and helps us understand the world around us. Besides my love for the world of mathematics, I love the outdoors, more specifically, mountain biking, backcountry camping, and skiing. I also enjoy taking photos, watching Formula 1, playing hockey, and improving my skills in the kitchen.

Continue reading below to find out more about my fellow Co-op students!

Development:

I’m Zhengmao (he/him), and I’m a third year in Software Engineering at the University of Waterloo. I’ll be working until the end of April here at Maplesoft as a Software Developer, where I’ll be working to fix bugs, add new features, and improve existing ones for our Maple Learn as well as Maple Calculator products. By the way, if you ever have any suggestions or ideas about them, don’t hesitate to reach out to me!

I’ve always been curious about working at a math company because I’ve always been so interested in math. In fact, Maplesoft is the only company I’ve consistently applied to every time I’ve gone through the Co-op application cycle! However, there’s not really any particular reason why I enjoy the subject. I find math to be beautiful in and of itself, almost like an art, and I find the kinds of math that are more discrete or algebraic tend to be a little nicer. As long as there aren’t decimals, I’m pretty happy. So, my ideal kind of math is just that: ideal! Exact values, unrealistic ideas, and as few numbers as possible. In terms of my university career, I’ve always enjoyed linear algebra much more than calculus.

Overall, I’m quite excited for this term at Maplesoft. I’ve never worked in web or mobile development before, so I’m looking forward to learning a lot of new things!

Content Creation:

Hi, I’m Paige (she/her). I am a second-year Honors Mathematics student at the University of Waterloo. This term, I am creating content for the Maple Learn document gallery. My favorite area of math is calculus because I love visualizing functions. Math matters because it is a universal language. All the math concepts we know are naturally occurring; people have observed and documented them, but no one invented them. Because of this, people from a wide range of cultures have come to the same conclusions (ex: defining pi). Math is universally understandable, which is why it can be used to connect everyone on earth (and maybe on other planets too!?!?!?!). In my free time, I like doing hand embroidery, playing video games, and cuddling my cat Licky.

A cat sitting on a window sill

Description automatically generated with low confidence

My name is Laura (she/her) and I’m a second year in the math program at the University of Waterloo. This term, I am working as a ‘Math Content Developer’ at Maplesoft; I’ll be creating and scripting documents for Maple Learn’s Example Gallery, updating older content, and handling customer requests. My favorite areas of math are probability, since I find questions like the Birthday Problem interesting, and biostatistics because I enjoy learning about biology and how biological experiments can be analyzed mathematically. I believe math matters because mathematics is essential to sending equipment and people into outer space; we will never meet aliens without using math.

Quality Assurance:

Hey, my name is Stefan, I'm 19 years old and currently studying Biochemistry in second year at the University of Waterloo. I am a QA analyst here at Maplesoft, working on Maple 2023. Outside of school, a hobby of mine is making digital art. My favorite area of math is definitely calculus & analysis because I found learning the fundamental theorem really intuitive and engaging. Expanding on that, I believe math matters because of its many applications in other fields such as the use of calculus in the research and design aspects of Biochemistry.

Hey! My name is Steven Mou, I'm in CFM at UWaterloo and I'm going into my 2B term after this term. I'm one of the four QA Analyst interns and I'll be testing all things related to Maple. My favorite area of math is anything related to algebra. I just find being able to manipulate variables while maintaining the integrity of the final product, to be very fascinating. I believe math matters because our lives are pretty much completely founded by math; anything from the technology that we use to the logic that is the foundation of our thoughts. I like playing sports, dabbling with different recipes in the kitchen, and discussing any shows that I happen to always finish too quickly.

Hi, my name is Aidan and I'm a 3rd-year mathematical physics student working as a quality assurance analyst for this Co-op term. As a Co-op student working in QA, I will mostly be running tests and reporting bugs to help ensure that Maple 2023 as well as Maple Flow are ready for release. My favorite area of math is vector calculus because as I started learning it I found it very interesting in the ways it applies to things we use in our everyday life. It also combines Linear Algebra and Calculus in a way that I never would have expected before learning about it. I think that math matters because I feel as though everything you can interact with can be described and predicted mathematically and that amazes me.

Hello! My name is Sebastian, I am currently in my second year of physics and astronomy at UWaterloo, and for this Co-op, at Maplesoft, I am working as a quality assurance analyst. In this position, I will be performing tests on Maple and Maple Flow to ensure that when they are released they function as they should and are ready for consumers to use with ease. When I am not focused on my academics I enjoy spending my time playing soccer (also watching it), listening to music, and watching movies. My favorite area of math is calculus because of the interesting and complex problems it provides, and because it is an essential tool needed to understand how the universe works. I believe that math matters because, as teachers always remind us, it is all around us. Math provides the foundation for everything we know and have come to appreciate in our lives, so since it is seen so often in our lives, I believe we should put in the effort to understand it and grasp how cool it is.

Can ChatGpt really help us with Maple programming?...

by: lcz 1039 Product: Maple

February 02 2023

7 5

Recently, ChatGpt has gained a lot of attention and there are mixed reviews about it. Some people think highly of its capabilities, while others are more critical. Regardless, it has been an interesting topic on assistant programming.

I'm going to give two examples to illustrate what's interesting and what's lacking.

First, I want ChatGpt to write a Maple code for the enumeration of all prime numbers up to 100. Then it writes:

primes:= proc(n)
   local i, j, primes;
   primes:= [2];
   for i from 3 to n do
      for j from 2 to i-1 do
         if i mod j = 0 then break;
      end do;
      if j = i then primes:= [op(primes), i]; end if;
   end do;
   return primes;
end proc;

primes(100);

I copy and paste this code, maple shows an error: Error, reserved word `do` unexpected

I report the error to it:

this code cannot run since "Error, reserved word `do` unexpected"

Then it gives: I apologize for the error in the previous code. Here is a corrected version of the code that can be run in Maple:

primes:= proc(n)
   local i, j, primesList;
   primesList:= [2];
   for i from 3 to n do
      for j from 2 to i-1 do
         if i mod j = 0 then 
            break;
         end if;
      end do;
      if j = i then 
         primesList:= [op(primesList), i]; 
      end if;
   end do;
   return primesList;
end proc;

primes(100);

Ok, It's tediously written code, but at least it's correct.

Furthermore, I ask:

Is there a more concise approach that utilizes maple's built-in functions.

It tells me:

Yes, there is a more concise approach that utilizes Maple's built-in functions. Here's an example:

primes := [seq(i, i=2..100)]; primes := [op(select(isprime, primes))];

The seq function generates a list of numbers from 2 to 100. The select function filters out all non-prime numbers using the isprime function, which tests the primality of a number. The result is stored in the primes array.

Ok, this answer meets my expectations.

Another example is the generation of all 6th order non-isomorphic graphs, which pops up as:

# Load the GraphTheory package
with(GraphTheory):

# Generate a list of all non-isomorphic graphs with 6 vertices
graphs := NonIsomorphicGraphs(6):

# Print the list of non-isomorphic graphs
print(graphs);

But NonIsomorphicGraphs(6) as we know is the number of non-isomorphic graphs but not graphs themselves.

So I say:

in the above codes, NonIsomorphicGraphs need the options output = graphs, outputform = graph

ChatGpt gives :

Here is a revised version of the code that takes into account the options output = graphs and outputform = graph:

# Load the GraphTheory package
with(GraphTheory):

# Generate a list of all non-isomorphic graphs with 6 vertices
graphsList := NonIsomorphicGraphs(6, output = graphs, outputform = graph):

# Print the number of non-isomorphic graphs
print("Number of non-isomorphic graphs: ", nops(graphsList));

# Print each non-isomorphic graph in the list
for i in 1 to nops(graphsList) do
  print("Graph ", i, ": ", graphsList[i]);
end do

It still has syntax errors. I made some artificial adjustments.

# Load the GraphTheory package
with(GraphTheory):

# Generate a list of all non-isomorphic graphs with 6 vertices
graphsList := [NonIsomorphicGraphs(6, output = graphs, outputform = graph)]:

# Print the number of non-isomorphic graphs
print("Number of non-isomorphic graphs: ", nops(graphsList));

# Print each non-isomorphic graph in the list
for i from 1 to nops(graphsList) do
  print("Graph ", i, ": ", graphsList[i]);
end do;

I pulled a wry face when it builds a lot of codes that do not exist in Maple.. Here's a little tidbit of my conversation with it.

This software still needs more work. I'm not sure how far it can go, but maybe programming becomes easier. Maybe more interesting examples will be found.

Regular Polyhedra and Schläfli Symbols

by: llayton 120 Product: Maple Learn

January 27 2023

3 1

When introduced to geometry, one of the first things we learn is the definition of the word “polygon”. A polygon is a closed 2-dimensional shape with at least 3 straight sides and angles. A regular polygon is a polygon with congruent sides and equal angles. A regular polygon with n sides has Schläfli symbol {n}. I’m interested in mathematical history, so when I learned that the idea of higher-dimensional spaces was invented in the middle of the nineteenth century I decided to research more about Ludwig Schläfli and the notation he came up with to describe his ideas.

In general, the Schläfli symbol is a notation of the form {p, q, r, ...} for regular polytopes. Polytopes are geometric objects with flat sides. This week, I will be focusing on 3-dimensional polytopes, also called polyhedra.

Similar to regular polygons, regular polyhedra are 3-dimensional shapes whose faces are all the same regular polygon. A regular polyhedron’s Schläfli symbol is of the form {p, q}, where p is the number of edges each face has and q is the number of faces that meet at each of the polyhedron’s vertices.

Below are two regular polyhedra: a cube (also known as a hexahedron) and a great stellated dodecahedron. The cube is one of five Platonic solids, and the great stellated dodecahedron is one of four Kepler-Poinsot polyhedra – all of these can be represented by Schläfli symbols. The cube has Schläfli symbol {4, 3}, since squares have 4 equal sides, and each vertex of a cube is created by the vertices of 3 squares meeting.

Can you figure out the Schläfli symbol for the great stellated dodecahedron?

The great stellated dodecahedron has the Schläfli symbol {5/2, 3}. This is because great stellated dodecahedrons are regular star polygons. As a result, the first number in their Schläfli symbol is an irreducible fraction whose numerator represents a number of sides and whose denominator corresponds to a turning number. The particular fraction 5/2 corresponds to a pentagram – a regular star polygon with 5 points – and great stellated dodecahedrons are composed of 12 of these pentagrams, where 3 pentagrams meet at each vertex of the shape.

One notable example of a regular polytope in pop culture is the tesseract, which has the Schläfli symbol {4, 3, 3}. This is an extension of the cube’s Schläfli symbol, {4, 3}, and the last number indicates that there are three cubes folded together around every edge. Below are two representations of a tesseract: one that uses a Schlegel diagram (left) and one from the 2012 movie Avengers (right).

Try out our Regular Polyhedra Visualization Using Schlafli Symbol Notation! In this document, you can test out your own Schläfli symbols for regular polyhedra. If they are valid Schläfli symbols, you’ll be provided with a 3-D visualization of the shape. If they are invalid, you can check out the logic for finding the specifications for regular polyhedra and this document, which provides all the 3-D regular polyhedra for you to try out.

Simplifying transfer functions with units and...

by: C_R 3442 Product: Maple

January 22 2023

1 3

Transfer functions are normally not used with units. Involving units when deriving transfer functions can help identify unit inconsistencies and reduce the likelihood of unit conversion errors.

Maple is already a great help in not having to do this manually. However, the final step of simplification still requires manual intervention, as shown in this example.

Given transfer function

H(s) = 60.*Unit('m'*'kg'/('s'^2*'A'))/(.70805*s^2*Unit('kg'^2*'m'^2/('s'^3*'A'^2))+144.*s*Unit('kg'^2*'m'^2/('s'^4*'A'^2))+0.3675e-4*s^3*Unit('kg'^2*'m'^2/('s'^2*'A'^2)))

H(s) = 60.*Units:-Unit(m*kg/(s^2*A))/(.70805*s^2*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*s*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(1)

Desired output (derived by hand) where the transfer function is separated in a dimensionless expression and a gain that can be attributed to units with a physical meaning in the context of an application (here displacement per voltage).

(2)

(3)

Units to factor out in the denominator are Unit('kg'^2*'m'^2/('s'^5*'A'^2)) . Quick check:

Unit('m'*'kg'/('s'^2*'A'))/Unit('kg'^2*'m'^2/('s'^5*'A'^2)) = Unit('m'/'V')

Units:-Unit(m*kg/(s^2*A))/Units:-Unit(kg^2*m^2/(s^5*A^2)) = Units:-Unit(m/V)

(4)

(5)

"Simplification" attempts with the denominator

s*(.70805*s*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(6)

collect(s*(.70805*s*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units[Unit](kg^2*m^2/(s^2*A^2))), Unit('kg'^2*'m'^2/('s'^5*'A'^2)))

s*(.70805*s*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(7)

is not effective because all units are wrapped in Unit commands. Example:

Units:-Unit(kg^2*m^2/(s^2*A^2))

(8)

Expand does not expand the argument of Unit commands.

>

Units:-Unit(kg^2*m^2/(s^2*A^2))

Units:-Unit(kg^2*m^2/s^2/A^2)

C1: Expanding Unit command

An expand facility could be a solution that expands a Unit command with combined units to a product of separate Unit commands.

When all units are expanded in a separate Unit command, collect or factor can be used to collect units:

.70805*s*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^3)+144.*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^4)+0.3675e-4*s^2*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^2)

.70805*s*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^3)+144.*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^4)+0.3675e-4*s^2*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^2)

(9)

(.70805*s/Units:-Unit(s)^3+144./Units:-Unit(s)^4+0.3675e-4*s^2/Units:-Unit(s)^2)*Units:-Unit(m)^2*Units:-Unit(kg)^2/Units:-Unit(A)^2

(10)

0.3675e-4*Units:-Unit(kg)^2*Units:-Unit(m)^2*(19266.66666*s*Units:-Unit(s)+3918367.346+.9999999999*s^2*Units:-Unit(s)^2)/(Units:-Unit(A)^2*Units:-Unit(s)^4)

(11)

C2: Using the Natural Units Environment

In this environment, no Unit commands are required and the collection of units should work with Maple commands.
However, for the expressions discussed here, this would lead to a naming conflict with the complex variable s of the transfer function and the unit symbol s for seconds.

C3: A type declaration or unit assumptions on names

A type declaration as an option of commands like in

$Units[TestDimensions](s*(.70805*s*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units[Unit](kg^2*m^2/(s^2*A^2))), {s::(Unit(1/s))})$

(12)

could help Maple in simplification tasks (in its general meaning of making expressions shorter or smaller).
Alternatively, assumptions could provide information of which "unit type" a name is

Error, (in assuming) when calling 'property/ConvertProperty'. Received: 'Units:-Unit(1/s) is an invalid property'

On various occasions (beyond transfer functions) I have looked for such a functionality.

C4: DynamicSystems Package with units

C4.1: The complex variable s could be attributed to the unit 1/s (i.e. Hertz) either by default or as an option. This could enable using units within the dynamic system package which is not possible in Maple 2022. An example what the package provides currently can be found here:
The phase plot shows that the package is already implicitly assuming that the unit of s is Hertz. A logical extension would be to have magnitude plots with units (e.g. m/V, as in this example).

C4.2: A dedicated "gain" command that takes units into account and that could potentially simplify the transfer function to an expression like (2) in SI units. In such a way the transfer function is separated into a dimensionless (but frequency depended) term and a gain term with units.
This would make the transfer of transfer functions to MapleSim easy and avoid unit conversion errors.

Download Collecting_and_expanding_units.mw

Happy Lunar New Year!

by: jthomson 90 Product: Maple Learn

January 20 2023

3 0

Happy Lunar New Year to everyone in the MaplePrimes community, as we enter the Year of the Rabbit. The rabbit symbolizes longevity, positivity, auspiciousness, wittiness, cautiousness, cleverness, deftness and self-protection!

To celebrate, one of our Maple Learn content developers, Laura Layton, made a Lunar New Year Color by Number:

Table

Description automatically generated

In this puzzle, your goal is to simplify the modulo equations in each square, and then fill in the square with the color that corresponds to the answer.

I hope you have fun solving this puzzle and revealing the hidden images and I wish everyone good health and happiness in the coming year!

Installing Maple on the Steam Deck

by: John May 3076 Product: Maple

January 19 2023

12 2

This year, for Christmas, my children^* got a Valve Steam Deck under the Christmas tree. It's a pretty cool device that looks a little like a monsterous Nintendo Switch, but it can run an impressive subset of the Steam video game catalog, games mostly designed to run on Windows PCs. It manages this by sporting a custom x86_64 processor by AMD and running a customized version of Arch Linux that uses Wine via Valve's Proton tool. The key point here, is that it is a tiny x86_64 compatible computer running Linux. So, of course, I needed to install Maple on it. So, I just paired a bluetooth keyboard, rebooted it into desktop mode and with a few small trick, bam, Maple on the Steam Deck:

There were a few small hiccups that required some work. I had absolutely no problems getting the Maple installer onto the device via a USB drive and no problems running it. I only ran into problems durring license activation:

Fortunately, I talked to our crack technical support team and they were able to identify this as a problem with Arch Linux not having full LSB 3.0 support installed by default. The process for fixing that is documented on the Arch Linux Wiki and involves just installing the ld-lsb package via pacman -- with the small additional wrinkle that you need to take the Steam Deck operating system out of 'read-only' mode in order to do that. But once that was done, I had a full version of Maple running well (albeit at 1280x800 resolution on a 7" display).

Since this device is designed for gaming, I was curious how fast it is compared to some other machines I work on. I chose an arbitrary benchmark of exactly solving a random linear system with integer coefficients.

restart;
N := 400;
A := LinearAlgebra:-RandomMatrix(N, N):
b := LinearAlgebra:-RandomVector(N):
v := [seq(cat(v__, i), i = 1 .. N)]:
sys := LinearAlgebra:-GenerateEquations(A, v, b):
CodeTools:-Usage(SolveTools:-LinearSolvers:-Rational(`~`[lhs - rhs](sys), v, dense = false)):

which it solves in decent time:

For comparison, this is 30% faster than the 32 core Xeon e5 workstation I do most of my work on, and only 5% slower than my notebook computer with an 8th gen Intel i7. Not bad for a toy! (please don't make me sad by telling me how much faster this is on a Mac M1 or M2)

Let me know in the comments if you have other benchmarks you want me to run on the Steam Deck. Also, please let me know if you manage to get your employer to buy you a Steam Deck to do scientific computing.

^*Okay, maybe it was a gift for me. Shhhh, don't tell.

The Moving Sofa Problem

by: jthomson 90 Product: Maple Learn

January 17 2023

4 -2

Last week, one of our Maple Learn developers, Valerie McKay-Crites, published a Maple Learn document, based on the very popular Maple application by Highschool Teacher, Jason Schattman called "Just Move It Over There, Dear!".

In the Maple application, Schattman explains the math behind moving a rectangular sofa down a hallway with a 90-degree turn. In the 3D Moving Sofa Problem Estimate, Valerie uses Schattman’s math to determine the largest rectangular sofa that can be taken down a flight of stairs and down a hallway with a 90-degree turn. Both applications reminded me of how interesting the Moving Sofa Problem is, which inspired me to write a blog post about it!

If you’ve ever been tasked with moving a rectangular sofa around a 90-degree turn, you might wonder:

What is the largest sofa that can make the move?

Following these steps as outlined in Schattman’s "Just Move It Over There, Dear!", will guarantee that the sofa will make the turn:

Measure the width of the hallway (h)
Measure the length (L) and width (w) of the sofa.
If L + 2w is comfortably less than triple the width of the hall, you'll make it!

When we work out the math exactly, we see that if the sofa's length plus double its width is less than 2*h*sqrt(2) , the sofa will make the turn!

Chart, line chart

Description automatically generated

This problem is easy if we only consider rectangular sofas, however, the problem becomes significantly more complex if we consider sofas of different shapes and areas. In mathematics, this problem is known as the Moving Sofa Problem, and it is unsolved. If we look at a hallway with a 90-degree turn and legs of width 1 m (i.e. h = 1 above), the largest known sofa that can make the turn is Gerver’s Sofa which has an area of 2.2195 m², this area is known as the Sofa Constant. Gerver’s Sofa, created in 1992, was constructed with 18 curve sections:

Icon

Description automatically generated

Check out this GIF of the sofa moving through the turn. It provides some insight into why Gerver’s sofa is such an interesting shape:

What is fascinating is that no mathematician has yet to prove that Gerver’s sofa is the sofa with the largest area capable of making the 90-degree turn.

The Moving Sofa Problem, is a great example of how math is embedded in our everyday lives. So, don’t stop being curious about the math around you as it can be fascinating and sometimes unproven!

If you are curious to learn more about the moving sofa problem check out this video by Numberphile, featuring Dan Romik from UC Davis: https://www.youtube.com/watch?v=rXfKWIZQIo4&t=1s

Happy New Years greetings to all maple users

by: wswain 240 Products: Maple , MapleSim

January 10 2023

1 0

Hello,

2022 was a wonderful year of progress in using Maple/MapleSim for almost everything my mathematical world.

I just wanted to wish all the Maplesoft user community a very productive and Happy New Year for 2023. I look forward to continue to find great nuggets of capability and insider techniques for using Maple in my endeavors.

Kindest Regards to ALL.
Happy New Year - 2023.
Bill

Maple proof of Poncelet's theorem for the triangle...

by: mmcdara 7896 Product: Maple

January 09 2023

3 0

This post is inspired by a serie of questions from @JAMET.
I wondered if it was possible to prove plane geometry theorems with the geometry package.

Here is an illustration for the Poncelet's theorem for the triangle (French designation), see for instance
https://en.wikipedia.org/wiki/Poncelet%27s_closure_theorem

Are any of you interested in challenging the geometry package to proof other plane geometry theorems?

>

restart:

Poncelet's theorem for the triangle

Let ABC a triangle, c its incircle (center omega, radius r) and C its circumcircle (center Omega, radius R).
Let D the distance between omega and Omega.

then R^2 - D^2 - 2*r*R = 0

Proof

Without loss of generality one sets :

x(A) = y(A) = 0
and y(B) = 0

ABC is a non degerated triangle provided x(B) <> 0 and y(C) <> 0

>	with(geometry):

>	kernelopts(version);

(1)

>	assume(x__B <> 0): assume(y__C <> 0):

>	point(A, 0, 0); point(B, x__B, 0); point(C, x__C, y__C);

(2)

>	triangle(T, [A, B, C])

(3)

>	bisector(bA, A, T); bisector(bB, B, T); eA := isolate(Equation(bA, [x, y]), y): eB := isolate(Equation(bB, [x, y]), y): xc := solve(rhs(eA)=rhs(eB), x): yc := eval(rhs(eA), x=xc): point(omega, xc, yc): r := distance(omega, line(lAB, [A, B]))

(abs(x__B)^2)^(1/2)*abs(x__B*y__C/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))/(x__B^2)^(1/2)

(4)

>	circumcircle(C, T, 'centername' = Omega); R := radius(C);

((1/4)*x__B^2+(1/4)*(-x__B*(x__C^2+y__C^2)+x__C*x__B^2)^2/(x__B^2*y__C^2))^(1/2)

(5)

>	Oo := distance(Omega, omega)

(((1/2)*x__B-(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2))/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))^2+(-(1/2)*(-x__B*(x__C^2+y__C^2)+x__C*x__B^2)/(x__B*y__C)-y__C*(x__B^2)^(1/2)/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))^2)^(1/2)

(6)

>	S := simplify(R^2 - Oo^2 - 2rR) assuming x__B::real, x__C::real, y__C::real

((x__B^2*(-abs(y__C)*signum(y__C)+y__C)*(x__C^2+y__C^2)^(1/2)+x__C^2*(y__C*abs(x__B)-signum(y__C)*abs(x__B*y__C)))*(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(-signum(y__C)*(x__C-x__B)^2*abs(x__B*y__C)+(abs(x__B)^2+x__C*(x__C-2*x__B))*abs(x__B)*y__C)*(x__C^2+y__C^2)^(1/2))/(y__C*((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+abs(x__B))^2)

(7)

>	simplify(S) assuming x__B > 0, y__C > 0; simplify(S) assuming x__B > 0, y__C < 0; simplify(S) assuming x__B < 0, y__C > 0; simplify(S) assuming x__B < 0, y__C < 0;

(8)

>

Download PoTh_proof.mw

Improvements of the geometry package

It already appears that (some) assumptions are not (always) correctly taken into account. This is a weak point which requires, as in the attached mw, to use an indirect approache to construct the incircle.

As a matter of fact, the procedure incircle, whose first lines are

showstat(incircle)

geometry:-incircle := proc(inci, T)
local cname, d, A, B, l1, l2, dis, x, y, tmp, msg;
   1   if nargs < 2 or 3 < nargs then
   2     error "wrong number of arguments"
       end if;
   3   if geometry:-form(T) <> ('triangle2d') then
   4     error "wrong type of arguments"
       end if;
   5   if nargs = 3 and op(1,args[3]) = ('centername') and type(op(2,args[3]),'name') then
   6     cname := op(2,args[3])
       else
   7     cname := cat('center_',inci)
       end if;
   8   if geometry:-method(T) = (':-points') then
   9     d := geometry:-DefinedAs(T);
  10     A := op(1,d);
  11     B := op(2,d);
  12     msg := sprintf("find the bisector of %a at vertex %a",T,A);
  13     userinfo(2,geometry,msg);
  14     geometry:-bisector('l1',A,T);
  15     msg := sprintf("find the bisector of %a at vertex %a",T,B);
  16     userinfo(2,geometry,msg);
  17     geometry:-bisector('l2',B,T);
  18     msg := sprintf("find the intersection of the two bisectors");
  19     userinfo(2,geometry,msg);
  20     geometry:-intersection(cname,l1,l2);

requires that the two bissectors are not parallel(call to geometry:-intersection).

Since the non-parallelism of bisectors is trivial for all non-degenerate triangles, why doesn't incircle inherit this property rather than not being able to decide if the bisectors are parallel or not?)

Here is a detail of what happens and the endless loop in which incircle seems to be caught

kernelopts(version);
                  Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895

AreParallel(bA, bB, 'condition'):
        AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero

assume(lhs(condition) <> 0);
AreParallel(bA, bB);
                             false
intersection(J, bA, bB);
                               J


infolevel[geometry] := 4:
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
intersection: one point of intersection
incircle: find the radius of the incircle
line: define the line from two points
circle: define the circle from its center and radius
circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
Error, (in geometry:-circle) not enough information: the radius might not be positive
assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0):

assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0): 
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero
intersection: two given lines intersect each other if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0
Error, (in geometry:-intersection) not enough information

assume(-y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0):
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
intersection: one point of intersection
incircle: find the radius of the incircle
line: define the line from two points
circle: define the circle from its center and radius
circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__C^2+y__C^2)^(1/2)+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
Error, (in geometry:-circle) not enough information: the radius might not be positive

Minimal edge cuts of a graph using fundamental...

by: dharr 8355 Product: Maple

January 08 2023

3 2

Recently here Mapleprimes member @Icz asked about generating all minimal edge cuts of a graph. I gave a brute force algorithm based on testng all bipartitions of the vertices. I then looked into improving the method by using fundamental cutsets to generate the cutsets. A description of the method and the ideas behind it is given here, together with a comparison of the two methods.

[Edit - see below for a newer version.]

Find all minimal edge cuts of a connected undirected graph using fundamental cutsets. (Assume no self-loops or multiple edges.)
The MinimalEdgeCuts procedure is in the startup code edit region (See Edit -> Startup Code.) (and is reproduced at the end of the worksheet).

>

Choose a graph.

>

G := Graph({{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}); DrawGraph(G, size = [200, 200])

The objective is to find all minimal edge cuts. If we partition the vertices of a graph in two non-empty sets (parts), then a cut is defined as a set of edges of the graph that have one end in the first part and the other end in the second part. Removal of these edges breaks the graph into two (or more) disconnected components.

For example the cut with edges {1,2} and {2,3} leaves two components with vertices [1,3,4] in one component and vertex [2] in the other component.

>

cut := {{1, 2}, {2, 3}}; IsCutSet(G, cut); DrawGraph(DeleteEdge(G, cut, inplace = false), size = [200, 200])

A cut is minimal if adding back any of the edges makes it connected again, i.e, the cut leave only two components. The above cut is minimal. Sometimes the term cutset is used to mean a minimal cut, but sometimes, as in Maple's IsCutSet, cutset means any removal of edges that increases the number of components.

This graph has 6 minimal cuts:

>

Vector[column](%id = 36893490386482949164)

This partitioning process suggests a brute force algorithm to find all minimal cuts. Find all partitions of the vertices into two (non-empty) sets. The two induced subgraphs and have some of the edges of the graph. All edges of the graph that are not in the two subgraphs form a cut. For example for the above partition we have:

>

partition := [[1, 3, 4], [2]]; G__1 := InducedSubgraph(G, partition[1]); 'Edges(G__1)' = Edges(G__1); G__2 := InducedSubgraph(G, partition[2]); 'Edges(G__2)' = Edges(G__2); cut := `minus`(`minus`(Edges(G), Edges(G__1)), Edges(G__2))

$GraphTheory:-Edges(G__1) = {{1, 3}, {1, 4}, {3, 4}}$

$GraphTheory:-Edges(G__2) = {}$

Now we have to check that each component is connected, or that there are a total of two components. Since there are in this case, the cut is minimal.

>

Note that we do have to test that there are only two components. For example for the path graph 1--2--3, the partition [[2],[1,3]] has three components (the three individual vertices). The associated cut {{1,2},{2,3}} is not minimal because we could put edge {1,2} back and the graph would still be disconnected.

The number of partitions to be tested is , where n is the number of vertices. Any subset of the vertices could be in the first set (except the empty set and the set of all vertices), and then the other vertices must be in the second set. However that will give double the possibilities, since swapping set 1 and set 2 makes no difference. In the present case there are possibilities, as seen below, where 0 indicates that the vertex in that position is in set 1 and 1 indicates that it is in set 2. For example, possibility 4 is the above partition .

>

1 2 3 4
1: 0 0 0 1
2: 0 0 1 0
3: 0 0 1 1
4: 0 1 0 0
5: 0 1 0 1
6: 0 1 1 0
7: 0 1 1 1

The algorithm described is implemented as procedure MinimalEdgeCutsPartition in the startup code. (The minimal cuts are not produced in the same order as above.)

However, a better approach is by building up the cuts fron a set of fundamental cuts (or cutsets, using the word in the restrictive sense). These are associated with a spanning tree. We choose any spanning tree, which we highlight with red edges. This is a tree graph (connected with no cycles) that has the same vertices as the graph.

>

Any spanning tree graph has edges. This is also the number of vertices minus the number of connected components (=1), which is the rank of the graph, which we will call .
We will find one fundamental cutset for each edge of the tree, so there will be fundamental cutsets.

>

NumberOfVertices(G)-1, NumberOfEdges(treeG), GraphRank(G)

Select one of the tree edges, say {1,2}. Removal of this edge partitions the tree into two components. We find the vertex partition and the corresponding cut by the procedure above. (We do not have to test for only two components because that follows from the tree structure.). This is the example we did above. The first three cuts in mincuts are the fundamental ones. Each one has one of the tree edges.

>

It is possible to generate all the cuts of the graph involving the tree edges by taking all subsets of the set of fundamental cutsets (this includes the empty set), and adding the elements of each subset with a ring sum operation.

A ring sum of two edge sets includes all edges in both sets except those that are common to both (the symmetric difference or disjunctive union):

>

(Maple has this as the symmdiff command, but not in infix form.)
We can also represent a set of edges by a vector, with entry 1 if the corresponding edge is in the set and 0 otherwise. The edge order we choose here is the order in . Then the ring sum operation is the same as addition modulo 2 of the corresponding vectors, though we won't make much use of this in implementing the algorithm. Set up to convert to vectors.

>

edges := Edges(G); e := NumberOfEdges(G); edgetable := table(`~`[`=`](edges, [`$`(1 .. e)])); tovec := proc (edges) options operator, arrow; Vector(e, `~`[`=`](map(proc (x) options operator, arrow; edgetable[x] end proc, edges), 1)) end proc

Let's try the ring sum on f[1] and f[2]. It is a cut, and since it now has two tree edges it is not a fundamental cut. (It happens to be minimal because it splits the graph into two components.)

>

f[1]*`&oplus;`*f[2] = `&oplus;`(f[1], f[2]); IsCutSet(G, rhs(%))

Do the same thing in the vector representation. The same logic about the tree edges (first three entries) shows the fundamental vectors are basis vectors and the three vectors are linearly independent (with respect to the addition mod 2 operation).

>

V[1], V[2], V[3] := Vector(5, {(1) = 1, (2) = 0, (3) = 0, (4) = 1, (5) = 0}), Vector(5, {(1) = 0, (2) = 1, (3) = 0, (4) = 1, (5) = 1}), Vector(5, {(1) = 0, (2) = 0, (3) = 1, (4) = 0, (5) = 1})

`mod`(V[1]+V[2], 2) = Vector[column](%id = 36893490386597633668)

The ring sum of any two cuts is another cut, which may be minimal (a cutset) or a disjoint union of cutsets (= two cutsets together). Combining the elements of the 8 subsets gives the eight cuts for this graph associated with the chosen tree graph, i.e., it has all combinations of tree edges. In vector form they are:

>

Vector[column](%id = 36893490386597646332), Vector[column](%id = 36893490386597647052), Vector[column](%id = 36893490386597647772), Vector[column](%id = 36893490386597648252), Vector[column](%id = 36893490386597648972), Vector[column](%id = 36893490386597649452), Vector[column](%id = 36893490386597649932), Vector[column](%id = 36893490386597650172)

Although we constructed these cuts relative to a specific spanning tree, they can be shown to be all the cuts of the graph. Consider the set of all edges of the graph. It might seem that this is a cut that has been missed, because it separates the graph into two or more (actually 4) components. However, looking at the edges

>

and recalling the definition of a cut, we see that this is not actually a cut. We need the edges in a cut all to have one end in one set of vertices and the other end in the set of other vertices, and this is not possible here. (Maple's IsCutSet is not this strict and returns true.) We can show it is not possible with IsBipartite.

>

Although we have found all cuts, we now need to test them to find which ones are minimal. The first case with no edges does not cut the graph into two components and we reject it. (It is usually included as a cut to make the cuts into a vector space.) The 6th one is not minimal. Its first and third entries are 1, meaning it is the ring sum of the first and third fundamental cuts. It splits the graph into three components:

>

This is a case where the ring sum is a disjoint union (not minimal), and arises because f[1] and f[3] have no edges in common:

>

For the present graph, all but the empty set and are minimal, and so there are 6 minimal cuts. In vector form:

>

Vector[column](%id = 36893490386597690532), Vector[column](%id = 36893490386597690652), Vector[column](%id = 36893490386597690772), Vector[column](%id = 36893490386597690892), Vector[column](%id = 36893490386597691012), Vector[column](%id = 36893490386597691132)

So the algorithm is to generate all cuts and test them for minimality. If we neglect the empty cut, then we have to generate cuts for testing. This is the same as the number of partitions that we generated in the partition algorithm. However, for the fundamental cuts we do not need to test that the components are connected, and more importantly, the ring sum is a more efficient way of finding the edges than partitioning and then using InducedSubgraphs.

The tests for minimality can be made more efficient by finding cases that can be classified without having to check the number of components. Two simple cases with modest savings (implemented here) are:

1. The fundamental cuts are minimal by construction.

2. If two cuts have no edges in common, then the ring sum is a disjoint union and is not minimal. This is the case for example for f[1] and f[3].

There may be other ways to classify without doing the more difficult minimality test, for example variations of (2.) for ring sums with more than two cuts. Hovever, the cutset algorithm is already significantly more efficient than the partition algorithm. Perhaps working with the vectors would be more efficient. Other more efficient algorithms may be known that I am not aware of.

The MinEdgeCuts algorithm in the startup code is reproduced here:

>

MinimalEdgeCuts:=proc(G::GRAPHLN)
  local mincutsets,n,edges,treeedges,treeG,
    i,j,r,vertexpartition,partitionedges,nf,fc,
    fsets,ref,fcutsets,cutsetedges;
  uses GraphTheory;
  if not IsConnected(G) then return Vector([]) end if;
  n:=NumberOfVertices(G);
  edges:=Edges(G);
  # choose a spanning tree and find corresponding fundamental cutsets
  treeG:=SpanningTree(G);
  treeedges:=Edges(treeG);
  r:=n-1;    # =nops(treeedges) = GraphRank(G)
  mincutsets:=table(); # first r ones are the fundamental ones
  for j to r do
    vertexpartition:=ConnectedComponents(DeleteEdge(treeG,treeedges[j],'inplace'=false));
    partitionedges:=`union`(map(Edges,map2(InducedSubgraph,G,vertexpartition))[]);
    mincutsets[j]:=edges minus partitionedges;
  end do;
  j:=r;
  # now find ringsums of all subsets of the set of fundamental cutsets
  fsets := Iterator:-BinaryGrayCode(r,'rank'=2); #skip empty set
  for ref in fsets do;
    fcutsets:=seq(ifelse(ref[i]=1,mincutsets[i],NULL),i=1..r);
    nf:=nops([fcutsets]);
    if nf=1 then
      next      # fundamental cutsets already in mincutsets
    elif nf=2 and (fcutsets[1] intersect fcutsets[2] = {}) then
      next      # pair of disjoint cutsets not minimal
    # todo other detections of non-minimal cutsets
    else        # rest the hard way
      cutsetedges:=fcutsets[1];
      for i from 2 to nf do
              fc:=fcutsets[i];
              cutsetedges:=(cutsetedges union fc) minus (cutsetedges intersect fc);
      end do;
      vertexpartition:=ConnectedComponents(DeleteEdge(G,cutsetedges,'inplace'=false));
      if nops(vertexpartition)=2 then # cutset is minimal
        j:=j+1;
        mincutsets[j]:=cutsetedges
      end if
    end if;
  end do;
  Vector(j,mincutsets);
end proc:

>

Download CutSets.mw

One of the basics of the Draghilev method

by: one man 1355 Product: Maple

January 06 2023

0 1

Does everyone have a good idea of the work of the Draghilev method? For example, in this answer https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions#answer291268
( https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions )
there was a very successful attempt by Rouben Rostamian to calculate the line of intersection of surfaces without applying the Draghilev method.
Let now not 3d, but 8d. And how will the solve command work in this case? Imagine that aij ((i=1..7,j=1..8)) are partial derivatives, and xj (,j=1..8) are derivatives, as in the above example. f8 is responsible for the parametrization condition.

restart;
 f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
 f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
 f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
 f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8;
 f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8; 
 f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
 f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;
 f8 := x1^2+x2^2+x3^2+x4^2+x5^2+x6^2+x7^2+x8^2-1; 
allvalues(solve({f1, f2, f3, f4, f5, f6, f7, f8}, {x1, x2, x3, x4, x5, x6, x7, x8}))

And this is how the Draghilev method works in this case.

restart; with(LinearAlgebra):
f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8; 
f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8;
f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;

n := 7; 
x := seq(eval(cat('x', i)), i = 1 .. n+1):
 F := [seq(eval(cat('f', i)), i = 1 .. n)]: 
A := Matrix(nops(F), nops(F)+1):
 for j to nops(F) do for i to nops(F)+1 do A[j, i] := op(1, op(i, op(j, F))) 
end do:
           end do: 

# b[i] and b[n+1] are solutions of a linear homogeneous system and at the 
# same time they serve as the right-hand sides of an autonomous ODE.
for i to n do

 b[i] := Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1)) 
                                                                     end do:
 b[n+1] := -Determinant(DeleteColumn(A, n+1)):

Only the original seven linear homogeneous equations with eight variables are needed. We solve them according to Cramer's rule, and in order to have uniqueness when solving the ODE, we use a point on the curve (according to the theory). (This point is sought in any convenient way.)
If we want to get a parameterization, then additionally, directly in dsolve, we can add the following:

for i to n do 
b[i] := simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))) 
end do:
b[n+1] := simplify(-Determinant(DeleteColumn(A, n+1))); 
deqs := seq(diff(x[i](s), s) = b[i]/(b[1]^2+b[2]^2+b[3]^2+b[4]^2+b[5]^2+b[6]^2+b[7]^2+b[8]^2)^.5, i = 1 .. n+1):

Maple should still try harder

by: sursumCorda 1254 Product: Maple 2022

January 06 2023

0 4

In an old post, vv reported a bug in simpl/max, which has been "fixed" in Maple 2018. However, it seem that such repairs are not complete enough.
For example, suppose it is required to find the (squared) distance between the origin and a point on x³ - x + y² = ⅓ which is closest to the origin. In other words, one needs to minimize x²＋y² among the points on this curve, i.e.,

extrema(x^2 + y^2, {x^3 + y^2 - x = 1/3}, {x, y}, 's'); # in exact form

Unfortunately, an identical error message appears again:

>

restart;

>	$extrema(x^2+y^2, {x^3+y^2-x = -2/(3*sqrt(3))}, {x, y})$

(1)

>	$extrema(x^2+y^2, {x^3+y^2-x = 1/3}, {x, y})$

Error, (in simpl/max) complex argument to max/min: 1/36*((36+12*I*3^(1/2))^(2/3)+12)^2/(36+12*I*3^(1/2))^(2/3)

>	$`~`[`^`](extrema(sqrt(x^2+y^2), {x^3+y^2-x = 1/3}, {x, y}), 2)$

(2)

>	$extrema(x^2+1/3-x^3+x, {x^3+y^2-x = 1/3}, {x, y})$

(3)

>	$MTM[limit](extrema(x^2+y^2, {x^3+y^2-x = a}, {x, y}), 1/3)$

(4)

Download tryHarder.mws

How about changing the values of parameter ?

for a from -3 by 3/27 to 3 do
    try
        extrema(x^2 + y^2, {x^3 + y^2 - x = a}, {x, y}); 
    catch:
        print(a); 
    end;
od;
                               -1
                               --
                               3 

                               -2
                               --
                               9 

                               -1
                               --
                               9 

                               1
                               -
                               9

                               2
                               -
                               9

                               1
                               -
                               3

By the way, like extrema, Student[MultivariateCalculus]:-LagrangeMultipliers also executes the Lagrange Multiplier method, but strangely,

Student[MultivariateCalculus][LagrangeMultipliers](y^2 + x^2, [x^3 + y^2 - x - 1/3], [x, y], output = plot):

does not cause any errors.

Integral Vector Calculus and parametrization of...

by: ecterrab 14540 Product: Maple

December 31 2022

9 6

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

Kudos to the people working on plots

by: rcorless 740 Product: Maple 2022

December 20 2022

2 0

I have been making animated 3d plots recently; the last time was perhaps three years ago, and I had some problems then. If I recall correctly, I couldn't make an animated 3d plot that was plotted in non-Cartesian coordinates.

I am very happy to report that this works very smoothly now in Maple 2022, and it's pretty fast, too. I have a fairly complex function to plot, involving piecewise polynomials on a tensor product grid in the xi and eta variables (actually, I let plot3d pick out the grid; it seems happier to do so) and then plot them on an elliptical base, in coordinates x = d*cosh(xi)*cos(eta) and y=d*sinh(xi)*sin(eta) (d is just a numerical constant, giving the location of the foci at (d,0) and (-d,0)), for 0 <= xi <= xi[0] (the outer elliptical boundary) and 0 <= eta <= 2Pi. The straightforward command works, and building a sequence of plots and using plots[display] works. I put option remember into my procedure w(xi,eta) and because the sample points are consistent for the time-dependent function exp(I*omega*t)*w(xi,eta) the xi-eta grid needs only to be done once and then one can compute (basically) as many frames as one wants in rapid succession.

Works great. Thanks, folks!

for k to nplots do
t := evalf(2*Pi*(k - 1)/nplots);
plts[k] := plot3d([(xi, eta) -> focus*cosh(xi)*cos(eta), (xi, eta) -> focus*sinh(xi)*sin(eta), (xi, eta) -> Re(exp(omega*t*I)*w(xi, eta))], 0 .. xi[0], 0 .. 2*Pi, colour = ((xi, eta) -> Re(exp(omega*t*I)*w(xi, eta))), style = surfacecontour, lightmodel = "none");
end do;
plots[display](seq(plts[k], k = 1 .. nplots), insequence = true);

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
Latest Posts Feed

Physics Courseware Support: Mechanics

Meet the Co-ops! - Winter 2023

Can ChatGpt really help us with Maple programming?...

Regular Polyhedra and Schläfli Symbols

Simplifying transfer functions with units and...

Happy Lunar New Year!

Installing Maple on the Steam Deck

The Moving Sofa Problem

Happy New Years greetings to all maple users

Maple proof of Poncelet's theorem for the triangle...

Minimal edge cuts of a graph using fundamental...

One of the basics of the Draghilev method

Maple should still try harder

Integral Vector Calculus and parametrization of...

Kudos to the people working on plots

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post Latest Posts Feed

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Latest Post
Latest Posts Feed