Applications, Examples and Libraries

3D plot of robot and Wankel engine using a single...

by: Mahardika Mathematics 35 Product: Maple 2016

April 17 2023

1 2

TODAY I GOT AN INSPIRATION TO CREATE 3D GRAPH EQUATION OF WALKING ROBOT (ED-209) IN CARTESIAN SPACE USING ONLY WITH SINGLE IMPLICIT EQUATION.

ENJOY...

How to Create Graph Equation of Wankel Engine on Cartesian Plane using Single Implicit Function run by Maple Software

Enjoy...

DMS degrees minute seconds to decimal conversion

by: Christopher2222 6035 Product: Maple

April 11 2023

1 2

It seems once in a while someone asks about converting from degrees minutes seconds format to decimal or the other way around.

I created a little procedure for just that purpose.

A little procedure to convert the form of degrees, minutes, seconds to decimal and vice versa.

$dms := proc (d, m := 0, s := 0) local con, d1, d2, d3; if 0 < frac(d) then d1 := floor(d); d2 := floor(60*frac(d)); d3 := 60*frac(60*frac(d)); con := cat(d1, `° `, d2, `' `, d3, `"`) else con := d+m/60.+s/3600. end if; print(con) end proc$

Examples of use.

(1)

(2)

Download dms.mw

edit - a quick realization is to remove the decimals that changes the fractions to floating decimals and change print(con) to print(evlaf(con)) to avoid rounding issues.

Solving a least-squares problem with QR decomposit...

by: mmcdara 7896 Product: Maple

April 11 2023

2 0

This post is motivated by a question asked by @vs140580 ( The program is making intercept zero even though There is a intercept in regression Fit (A toy code showing the error attached) ).

The problem met by @vs140580 comes from the large magnitudes of the (two) regressors and the failure to Fit/LinearFit to find the correct solution unless an undecent value of Digits is used.
This problem has been answerd by @dharr after scaling the data (which is always, when possible, a good practice) and by
myself while using explicitely the method called "Normal Equations" (see https://en.wikipedia.org/wiki/Least_squares).

The method of "Normal Equations" relies upon the inversion of a symmetric square matrix H whose dimension is equal to the number of coefficients of the model to fit.
It's well known that this method can potentially lead to matrices H extremely ill-conditionned, and thus face severe numerical problems (the most common situation being the fit of a high degree polynomial).

About these alternative methods:

In English: http://www.math.kent.edu/~reichel/courses/intr.num.comp.1/fall09/lecture4/lecture4.pdf
In French: https://moodle.utc.fr/pluginfile.php/24407/mod_resource/content/5/MT09-ch3_ecran.pdfI

The attached file illustrates how the QR decomposition method works.
The test case is @vs140580's.

Maybe the development team could enhance Fit/LinearFit in future versions by adding an option which specifies what method is to be used?

>

restart:

>	with(Statistics):

>	interface(version)

(1)

>

Data := Matrix([[4.74593554708566, 11385427.62, 2735660038000], [4.58252830679671, 25469809.77, 12833885700000], [4.29311160501838, 1079325200, 11411813200000000], [4.24176959154225, 1428647556, 18918585950000000], [5.17263072694618, 1428647556, 18918585950000000], [4.39351114955735, 1877950416, 30746202150000000], [4.39599006758777, 1428647556, 18918585950000000], [5.79317412396815, 2448320309, 49065217290000000], [4.48293612651735, 2448320309, 49065217290000000], [4.19990181982522, 2448320309, 49065217290000000], [5.73518217699046, 1856333905, 30648714900000000], [4.67943831980476, 3071210420, 75995866910000000], [4.215240105336, 2390089264, 48670072110000000], [4.41566877563247, 3049877383, 75854074610000000], [4.77780395369828, 2910469403, 74061327950000000], [4.96617430604669, 1416936352, 18891734280000000], [4.36131111330988, 1416936352, 18891734280000000], [5.17783192063198, 1079325200, 11411813200000000], [4.998266287191, 1067513353, 11402362980000000], [4.23366152474871, 2389517120, 48661380410000000], [4.58252830679671, 758079709.3, 5636151969000000], [6.82390874094432, 1304393838, 14240754750000000], [4.24176959154225, 912963601.2, 8621914602000000], [4.52432881167557, 573965555.4, 3535351888000000], [4.84133601918601, 573965555.4, 3535351888000000], [6.88605664769316, 732571773.2, 5558875538000000], [4.35575841415627, 1203944381, 13430693320000000], [4.42527441640593, 955277678, 8795128298000000], [6.82390874094432, 997591754.9, 8968341995000000], [4.35144484433733, 143039477.1, 305355143300000]]):

>

# Direct use of LinearFit.
#
# As far as I know LinearFit is based on the resolution of the "Normal Equations"
# (see further down), a system of equations that is known to be ill-conditioned
# when regressors have large values (in particular when polynomial regression
# is used).

X := Data[.., [2, 3]]:
Y := Data[.., 1]:

LinearFit(C1+C2*v+C3*w, X, Y, [v, w]);

Warning, model is not of full rank

(2)

>

# For roundoff issues the 3-by-3 matrix involved in the "Normal Equations" (NE)
# appears to of rank < 3.
# The rank of this matrix is rqual to 1+rank(X) and one can easily verify that
# the 2 columns of X are linearly independent:

LinearAlgebra:-LinearSolve(X, Vector(numelems(Y), 0));
LinearAlgebra:-Rank(X);

(3)

>	# Solve the least squares problem by using explicitely the NE. # # To account for an intercept we augment X by a vector column of "1" # traditionally put in column one. Z := `<\|>`(Vector(numelems(Y), 1), X): A := (Z^+ . Z)^(-1) . Z^+ . Y; # Normal Equations

A := Vector(3, {(1) = 4.659353816079307, (2) = 0.5985084089529947e-9, (3) = -0.27350964718426345e-16}, datatype = float[8])

(4)

>	# What is the rank of Z? # Due to the scale of compared to "1", Rank fails to return the good value # of rank(Z), which is obviously equal to rank(X)+1. LinearAlgebra:-LinearSolve(Z, Vector(numelems(Y), 0)); LinearAlgebra:-Rank(Z);

(5)

A WORKAROUND : SCALING THE DATA

>	model := unapply( LinearFit(C1+C2v+C3w, Scale(X), Scale(Y), [v, w]), [v, w] );

(6)

>	mX, sX := (Mean, StandardDeviation)(X); mY, sY := (Mean, StandardDeviation)(Y);

(7)

>	MODEL := model((x1-mX[1])/sX[1], (x2-mX[2])/sX[2]) * sY + mY

(8)

>	# Check that the vector of regression coefficients is almost equal to A found above # relative error lesst than 10^(-14) A_from_scaling := < coeffs(MODEL) >: Relative_Discrepancy := (A_from_scaling - A) /~ A

$Relative_Discrepancy := Vector(3, {(1) = 0.5718679809000842e-15, (2) = 0.14166219140514066e-13, (3) = 0.14308415396983588e-13}, datatype = float[8])$

(9)

THE QR DECOMPOSITION (applied on raw data)

The QR decomposition, as well as Given's rotation method, are two alternatives to the the NE method
to find the vector of regression coefficients.
Both of them are known to be less sensitive to the magnitudes of the regressors and do nt require (not
always) a scaling of the data (which can be quite complex with polynomial regression or when some
transformation is used to liearize the statistical model, for instanc Y=a*exp(b*X) --> log(Y)=log(a)+b*X).

>	N := numelems(Y); P := numelems(Z[1]);

(10)

>	# Perform the QR decomposition of Z. Q, R := LinearAlgebra:-QRDecomposition(Z, fullspan);

$Q, R := Vector(4, {(1) = ` 30 x 30 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}), Vector(4, {(1) = ` 30 x 3 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*triangular[upper], (4) = `Order: `*Fortran_order})$

(11)

>	# Let C the column vector of length P defined by: C := (Q^+ . Y)[1..P];

C := Vector(3, {(1) = -26.6044149698075, (2) = -.517558353158176, (3) = -.881007519371895})

(12)

>	# Then the vector of regression coefficients is given by: A_QR := (R[1..P, 1..P])^(-1) . C; Relative_Discrepancy := (A_QR - A) /~ A

$A_QR := Vector(3, {(1) = 4.65935381607931, (2) = 0.5985084090e-9, (3) = -0.2735096472e-16})$

$Relative_Discrepancy := Vector(3, {(1) = 0.3812453206e-15, (2) = 0.1105656128e-13, (3) = 0.1216778632e-13})$

(13)

>	# The matrix H = Z^+ . Z writes H := Z^+ . Z: H_QR := R^+ . Q^+ . Q . R: Relative_Discrepancy := (H_QR - H) /~ H

$Relative_Discrepancy := Matrix(3, 3, {(1, 1) = -0.1184237893e-15, (1, 2) = 0., (1, 3) = -0.3491343943e-15, (2, 1) = 0., (2, 2) = 0.1930383052e-15, (2, 3) = 0.3369103254e-15, (3, 1) = -0.1745671971e-15, (3, 2) = -0.5053654881e-15, (3, 3) = -0.1366873891e-15})$

(14)

>	# H_QR expression is required to obtain the covariance matrix of the regression coefficients.

Download LeastSquares_and_QRdecomposition.mw

plot of dna helix using as parametric surface

by: Mahardika Mathematics 35 Product: Maple 2016

April 09 2023

1 1

CREATING GRAPH EQUATION OF "DNA" IN CARTESIAN SPACE USING PARAMETRIC SURFACE EQUATION RUN ON MAPLE SOFTWARE

ENJOY...

Easter Eggs

by: pstone 205 Product: Maple Learn

April 05 2023

3 0

Happy Springtime to all in the MaplePrimes Community! Though some in our community may not live in the northern hemisphere where flowers are beginning to bloom, many will be celebrating April holidays like Ramadan, Passover, and Easter.

One of my favorite springtime activities is decorating eggs. Today, the practice is typically associated with the Christian holiday of Easter. However, painted eggs have roots in many cultures.

For over 3,000 years, painting eggs has been a custom associated with the holiday of Nowruz, or Persian New Year, during the spring equinox. Furthermore, in the Bronze Age, decorated ostrich eggs were traded as luxury items across the Mediterranean and Northern Africa. Dipped eggs have also played an important role in the Jewish holiday of Passover since the 16th century.

To celebrate this tradition, I would like to invite all of the Maplesoft community to create a decorated egg of their own with the Easter Egg Art Maple Learn document. In this document, an ovoid egg equation is used to define the shape of an egg.

The ovoid egg equation mimics the shape of a typical hen’s egg. Each bird species lays differently shaped eggs. For example, an ostrich’s egg is more oblong than an owl’s, and an owl’s egg is rounder than a goose’s. Surprisingly, every egg can be described by a single equation with four parameters:

Learn more about this equation and others like it with John May’s Egg Formulas Maple Learn document.

The Easter Egg Art document includes 9 different decorative elements; users can change the color, position, and size of each in order to create their own personal egg! The egg starts out looking like this:

In just a couple of minutes, you can create a unique egg. Have fun exploring this document and share a screenshot of your egg in the comments below! Here’s one I made:

Creating Graph Equation of Water Drop Wave in...

by: Mahardika Mathematics 35 Product: Maple 2016

April 03 2023

1 0

How to Create Graph Equation of Water Drop Wave in Cartesian Space using single Implicit Function only run by Maple Software

The Equation is: z = - cos( (x²+y²)^0.5 - a) with paramater a is moving form 0 to 2pi

Enjoy...

Plese Click link below to see full equation in Maple software

Water_Drop_Wave.mw

plots using only single implicit functions

by: Mahardika Mathematics 35 Product: Maple 2016

March 29 2023

1 5

Creating Graph Equation of An Apple in Cartesian Space using single Implicit Function only run by Maple software

Enjoy...

Please click the link below to see full equation on Maple file:

2._Apel_3D_A.mw

Creating Graph Equation of A Candle on Cartesian Plane using single Implicit Function only run by Maple software

Enjoy...

3D_Candle.mw

Sea_Shells.mw

Today I'm very greatfull to have Inspiration to create Graph Equation of 3D Candle in Cartesian Space using single 3D Implicit Function only, run by Maple software.

Enjoy...

Candle_1.mw

Today I got an inspiration to create graph equation of "Petrol Truck" using only with Single Implicit Equation in Cartesian space run by Maple Software

Maple software is amazing...

Enjoy...

CREATING 3D GRAPH EQUATION OF BACTERIOPHAGE USING ONLY WITH SINGLE IMPLICIT EQAUTION IN CARTESIAN SPACE RUN BY MAPLE SOFTWARE

MAPLE SOFTWARE IS AMAZING...

ENJOY...

GRAPH EQUATION OF A FEATHER

AND THE EQUATION IS:

ENJOY...

I like this Equation and post it because it is so beautiful...

Click this link below to see full equation and download the Maple file:

Bulu_Angsa_3.mw

GRAPH EQUATION OF "383" CREATED BY DHIMAS MAHARDIKA

ENJOY...

Download 383.mw

Drawing Eifel Tower using Implicit Equation in Cartesian Space

What is the true equation your time and space...

by: mmcdara 7896 Product: Maple

March 23 2023

3 0

This post is closely related to this one Centered Divided Difference approximations whose purpose is to build Finite-Difference (FD) approxmation schemes.
In this latter post I also talked, without explicitely naming it, about Truncation Error, see for instance https://en.wikiversity.org/wiki/Numerical_Analysis/Truncation_Errors.

I am foccusing here on a less known concept named "Equivalent Equation" (sometimes named "Modified Equation").
The seminal paper (no free acccess, which is surprising since it was first published 50 years ago) is this one by Warming and Hyett https://www.sciencedirect.com/science/article/pii/0021999174900114.
For a scholar example you can see here https://guillod.org/teaching/m2-b004/TD2-solution.pdf.
An alternative method to that of Warming and Hyett to derive the Equivalent Equation is given here (in French)
http://www.numdam.org/item/M2AN_1997__31_4_459_0.pdf (Carpentier et al.)

I never heard of the concept of Modified Equation applied to elliptic PDEs ; it's main domain of application is advection PDEs (parabolic PDEs in simpler cases).

Basically, any numerical scheme for solving an ODE or a PDE has a Truncation Error: this means there is some discrepancy between the exact solution of this PDE and the solution the scheme provides.
This discrepancy depends on the truncation error, in space alone for ODEs or elliptic PDEs, or in space and time for hyperbolic or advection PDEs.

One main problem with the Truncation Error is that it doesn't enable to understand the impact it will have on the solution it returns. For instance, will this sheme introduce diffusion, antidiffusion, scattering, ...
The aim of the Modified Equation is to answer these questions by constructing the continuous equation a given numerical scheme solves with a zero truncation error.
This is the original point of view Warming and Hyett developped in their 1974 paper.

This is a subject I work on 30 years ago (Maple V), and later in 2010.
It is very easy with Maple to do the painstaking development that Warming and Hyett did by hand half a century ago. And it is even so easy that the trick used by Carpentier et al. to make the calculations easier has lost some of its interest.

Two examples are given in the attched file
EquivalentEquation.mw

If a moderator thinks that this post should be merged with Centered Divided Difference approximations, he can do it freely, I won't be offended.

Centered Divided Difference approximations

by: mmcdara 7896 Product: Maple

March 21 2023

3 9

This code enables building Centered Divided Difference (CDD) approximations of derivatives of a univariate function.
Depending on the stencil we choose we can get arbitrary high order approximations.

The extension to bivariate functions is based upon what is often named tensorization in numerical analysis: for instance diff(f(x, y), [x, y] is obtained this way (the description here is purely notional)

Let CDD_x the CDD approximation of diff( f(x), x) ) .
CDD_x is a linear combination of shifted replicates of f(x)
Let s one of this shifted replicates
Let CDD_y(s) the CDD approximation of diff( s(y), y) ) .
Replace in CDD_x all shifted replicates by their corresponding expression CDD_y(s)

REMARKS:

When I write for instance "approximation of diff(f(x), x)", this must be intended as a short for "approximation of diff(f(x), x) at point x=a"
I agree that a notation like, for instance, diff(f(a), a) is not rigourous and that something like a Liebnitz notation would be better. Unfortunately I don't know how to get it when I use mtaylor.

>

restart:

CDDF stands for Cendered Divided Difference Formula

>

CDDF := proc(f, A, H, n, stencil)
  description "f = target function,\nA = point where the derivatives are approximated,\nH = step,\nn = order of the derivative,\nstencil = list of points for the divided differenceCDDF\n";
  local tay, p, T, sol, unknown, Unknown, expr:

  tay := (s, m) -> convert(
                     eval(
                       convert(
                         taylor(op(0, f)(op(1, f)), op(1, f)=A, m),
                         Diff
                       ),
                       op(1, f)=A+s*H),
                     polynom
                   ):

  p   := numelems(stencil):
  T   := add(alpha[i]*tay(i, p+1), i in stencil):
  T   := convert(%, diff):

  if p > n+1 then
    sol := solve([seq(coeff(T, h, i)=0, i in subsop(n+1=NULL, [$0..p]))], [seq(alpha[i], i in stencil)])[];
  else
    sol := solve([seq(coeff(T, H, i)=0, i in subsop(n+1=NULL, [$0..n]))], [seq(alpha[i], i in stencil)])[];
  end if:

  if `and`(is~(rhs~(sol)=~0)[]) then
    WARNING("no solution found"):
    return
  else
    unknown := `union`(indets~(rhs~(sol))[])[];
    Unknown := simplify(solve(eval(T, sol) = Diff(op(0, f)(A), A$n), unknown));
    sol     := lhs~(sol) =~ eval(rhs~(sol), unknown=Unknown);
    expr    := normal(eval(add(alpha[i]*op(0, f)(A+i*H), i in stencil), sol));
  end if:

  return expr
end proc:

>	Describe(CDDF)

# f = target function,
# A = point where the derivatives are approximated,
# H =
# step,
# n = order of the derivative,
# stencil = list of points for the divided
# differenceCDDF
#
CDDF( f, A, H, n, stencil )

>	# 2-point approximation of diff(f(x), x) \| x=a CDDF(f(x), a, h, 1, [-1, 1]); convert(simplify(mtaylor(%, h=0, 4)), Diff);

Diff(f(a), a)+(1/6)*(Diff(Diff(Diff(f(a), a), a), a))*h^2

(1)

>	# 3-point approximation of diff(f(x), x$2) \| x=a CDDF(f(x), a, h, 2, [-1, 0, 1]); convert(simplify(mtaylor(%, h=0)), Diff);

Diff(Diff(f(a), a), a)+(1/12)*(Diff(Diff(Diff(Diff(f(a), a), a), a), a))*h^2

(2)

>	# 5-point pproximation of diff(f(x), x$2) \| x=a CDDF(f(x), a, h, 2, [$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);

Diff(Diff(f(a), a), a)-(1/90)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a))*h^4

(3)

>	# 7-point approximation of diff(f(x), x$2) \| x=a CDDF(f(x), a, h, 2, [$-3..3]); # simplify(taylor(%, h=0, 10)); convert(simplify(mtaylor(%, h=0, 10)), Diff);

Diff(Diff(f(a), a), a)+(1/560)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a), a))*h^6

(4)

>	# 4-point staggered approximation of diff(f(x), x$3) \| x=a CDDF(f(x), a, h, 3, [seq(-3/2..3/2, 1)]); convert(simplify(mtaylor(%, h=0, 6)), Diff);

-(f(a-(3/2)*h)-3*f(a-(1/2)*h)+3*f(a+(1/2)*h)-f(a+(3/2)*h))/h^3

Diff(Diff(Diff(f(a), a), a), a)+(1/8)*(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a))*h^2

(5)

>	# 6-point staggered approximation of diff(f(x), x$3) \| x=a CDDF(f(x), a, h, 3, [seq(-5/2..5/2, 1)]); # simplify(taylor(%, h=0, 8)); convert(simplify(mtaylor(%, h=0, 8)), Diff);

(1/8)*(f(a-(5/2)*h)-13*f(a-(3/2)*h)+34*f(a-(1/2)*h)-34*f(a+(1/2)*h)+13*f(a+(3/2)*h)-f(a+(5/2)*h))/h^3

Diff(Diff(Diff(f(a), a), a), a)-(37/1920)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a))*h^4

(6)

>	# 5-point approximation of diff(f(x), x$4) \| x=a CDDF(f(x), a, h, 4, [$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);

Diff(Diff(Diff(Diff(f(a), a), a), a), a)+(1/6)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a))*h^2

(7)

>	# 7-point approximation of diff(f(x), x$4) \| x=a CDDF(f(x), a, h, 4, [$-3..3]); convert(simplify(mtaylor(%, h=0, 10)), Diff);

Diff(Diff(Diff(Diff(f(a), a), a), a), a)-(7/240)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a), a))*h^4

(8)

A FEW 2D EXTENSIONS

>

# diff(f(x, y), [x, y]) approximation over a (2 by 2)-point stencil

stencil := [-1, 1]:

# step 1: approximate diff(f(x, y), x) over stencil < stencil >

fx := CDDF(f(x), a, h, 1, stencil):
fx := eval(% , f=(u -> f[u](y))):
ix := [indets(fx, function)[]]:

# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.

fxy := add(map(u -> CDDF(u, b, k, 1, stencil)*coeff(fx, u), ix)):

# step 3: rewrite fxy in a more convenient form

[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );

convert(mtaylor(fxy, [h=0, k=0]), Diff)

(9)

>

# Approximation of diff(f(x, y), [x, x, y, y] a (3 by 3)-point stencil

stencil := [-1, 0, 1]:

# step 1: approximate diff(f(x, y), x) over stencil < stencil >

fx := CDDF(f(x), a, h, 2, stencil):
fx := eval(% , f=(u -> f[u](y))):
ix := [indets(fx, function)[]]:

# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.

fxy := add(map(u -> CDDF(u, b, k, 2, stencil)*coeff(fx, u), ix)):

# step 3: rewrite fxy in a more convenient form

[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );

convert(mtaylor(fxy, [h=0, k=0], 8), Diff)

-(2*f(a, b-k)-4*f(a, b)+2*f(a, b+k)-f(a-h, b-k)+2*f(a-h, b)-f(a-h, b+k)-f(a+h, b-k)+2*f(a+h, b)-f(a+h, b+k))/(h^2*k^2)

(10)

>

# Approximation of diff(f(x, y), [x, x, y] a (3 by 2)-point stencil

stencil_x := [-1, 0, 1]:
stencil_y := [-1, 1]:

# step 1: approximate diff(f(x, y), x) over stencil < stencil >

fx := CDDF(f(x), a, h, 2, stencil_x):
fx := eval(% , f=(u -> f[u](y))):
ix := [indets(fx, function)[]]:

# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.

fxy := add(map(u -> CDDF(u, b, k, 1, stencil_y)*coeff(fx, u), ix)):

# step 3: rewrite fxy in a more convenient form

[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );

convert(mtaylor(fxy, [h=0, k=0], 6), Diff)

(11)

>

# Approximation of the laplacian of f(x, y)

stencil := [-1, 0, 1]:

# step 1: approximate diff(f(x, y), x) over stencil < stencil >

fx := CDDF(f(x), a, h, 2, stencil):
fy := CDDF(f(y), b, k, 2, stencil):

fxy := simplify( eval(fx, f=(u -> f(u, b))) + eval(fy, f=(u -> f(a, u))) );

convert(mtaylor(fxy, [h=0, k=0], 6), Diff)

(12)

>

Download CDD.mw

Shortest walks using LinearAlgebra:-Generic

by: dharr 8597 Product: Maple

March 21 2023

3 4

Recently Mapleprimes user @vs140580 asked here about finding the shortest returning walk in a graph (explained more below). I provided an answer using the labelled adjacency matrix. @Carl Love pointed out that the storage requirements are significant. They can be reduced by storing only the vertices in the walks and not the walks themselves. This can be done surprisingly easily by redefining how matrix multiplication is done using Maple's LinearAlgebra:-Generic package.

(The result is independent of the chosen vertex.)

>

Consider the following graph. We want to find, for a given vertex , the shortest walk that visits all vertices and returns to . A walk traverses the graph along the edges, and repeating an edge is allowed (as distinct from a path, in which an edge can be used only once).

>

As is well known, the () entry of gives the number of walks of length between vertices and . The labelled adjacency matrix shows the actual walks.

>

Matrix(%id = 36893490455680637396)

For example, the (3,3) entry of has 3 terms, one for each walk of length 2 that starts and ends at vertex 3. The factors in a term give the edges visited.

So the algorithm to find the shortest walk is to raise to successively higher powers until one of the terms in the diagonal entry for the vertex of interest has indices for all the vertices.

>

Matrix(%id = 36893490455709504684)

The expressions for higher powers get very large quickly, so an algorithm that only retains the sets of vertices in each term as a set of sets will use less storage space. So we can consider the following modified labelled adjacency matrix.

>

Matrix(%id = 36893490455703852204)

Now we need to modify how we do matrix multiplication, but Maple has the LinearAlgebra:-Generic package to do this. We can redefine addition and multiplication to operate correctly on the sets of sets.

Consider two sets of sets of vertices for walks.

>

Addition is just combining the possibilities, and set union will do this. And addition of "zero" should add no possibilities, so we take {} as zero.

>

Multiplication is just combining all pairs by union. Notice here that {7} union {1,3,5} and {1,5} union {3,7} give the same result, but that we do not get duplicates in the set.

>

The unit for multiplication should leave the set of sets unchanged, so {{}} can be used

>

And we should check that zero multiplied by is zero

>

Define these operations for the LinearAlgebraGeneric package:

>

F[`=`] := `=`; F[`/`] := `/`; F[`0`] := {}; F[`1`] := {{}}; F[`+`] := `union`; F[`*`] := proc (x, y) options operator, arrow; {seq(seq(`union`(i, j), `in`(i, x)), `in`(j, y))} end proc

Warning, (in F[*]) `j` is implicitly declared local

Warning, (in F[*]) `i` is implicitly declared local

Compare with . We have lost information about the details of the walks except for the vertices visited.

>

Matrix(%id = 36893490455680647164)

So here is a procedure for finding the length of the shortest walk starting and ending at vertex v.

>

WalkLength:=proc(G::Graph,v)
  uses GraphTheory;
  local x,y,i,j,vv,A,B,F,n,vertset;
  if IsDirected(G) then error "graph must be undirected" end if;
  if not member(v,Vertices(G),'vv') then error "vertex not in graph" end if;
  if not IsConnected(G) then return infinity end if;
  F[`=`]:=`=`:F[`/`]:=`/`: # not required but have to be specified
  F[`0`]:={}:
  F[`1`]:={{}}:
  F[`+`]:=`union`;
  F[`*`]:=(x,y)->{seq(seq(i union j, i in x), j in y)};
  n:=NumberOfVertices(G);
  vertset:={$(1..n)};
  A:=Matrix(n,n, (i, j)-> if AdjacencyMatrix(G)[i, j] = 1 then {{i, j}} else {} end if);
  B:=copy(A);
  for i from 2 do
    B:=LinearAlgebra:-Generic:-MatrixMatrixMultiply[F](B,A);
  until vertset in B[vv,vv];
  i;
end proc:

>

Download WalksGenericSetOfSets.mw

Simple Interest Collection

by: 2teachis2learn4ever 60 Product: Maple Learn

February 16 2023

3 0

Hello everyone! Alex, Sarah, and I decided to create this collection of financial literacy documents as we noticed a lack of resources for this strand in mathematics. With many curricula around the world implementing financial literacy concepts, we thought it might be useful not just for Ontario, but for many jurisdictions around the world.

There are 4 documents in the Simple Interest collection; Introduction, Equation Generator, Mental Calculations, and Reflection. The Introduction is designed for intermediate and advanced level students as it introduces students to the concept of interest and how to calculate it. Students get to fill in the table by filling in the calculations on the right. This provides enough scaffolding so students of various grades can participate in this activity.

The Equation Generator document uses sliders to help students investigate linear equations in the form of y=mx+b. It also relates the simple interest equation (I=Prt) to the linear equation by asking students to compare interest rates. The idea behind this document is to bridge concepts outlined in the 2021 grade 9 destreamed math curriculum; in particular, the financial literacy, and linear relations strands. The document provides some reflection questions for students to think about the relationship between the variables.

The third document in the collection is the mental calculations document which presents a series of questions in increasing difficulty designed to help students compare interest rates. Students are intended to choose which scenario they think is more appropriate without using a calculator. There are hints provided on the right side if students wish for a hint, as well as explanations further to the right of the hints and answers below the main questions. Through our analysis of the curricula around the world, we noticed that many jurisdictions focus on mental math as a skill that their students should develop. Students may not always have access to a calculator and it is important for them to know how to make financially sound decisions or analyze advertisements that they may see around their neighbourhood.

Lastly, the last document is the reflection page where students are able to analyze their findings. In particular, “interest” may carry a negative connotation for students such that we want them to think of the potential benefits of interest as well. The reflection questions are designed to help students consolidate their learnings and can be further expanded on by the teacher. Such possibilities can include scenario-based questions.

May you find these documents helpful!

Physics Courseware Support: Mechanics

by: ecterrab 14630 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

Installing Maple on the Steam Deck

by: John May 3081 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.

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 8597 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

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