C_R

3412 Reputation

21 Badges

5 years, 313 days

MaplePrimes Activity


These are Posts that have been published by C_R

Plots of physical quantities has significantly improved with Maple 2022. The updated useunits option makes unit conversion errors in plots very unlikely. A lot of time is saved when creating plots of physical quantities where values and units must be correct.

One final source of user errors remains: The manual entry of incorrect units in labels.

Below is a way to avoid such errors by computing labels with units for three prevalent axis labeling schemes.


Other desireable labels are given as a suggestion for future plot label enhancements where plot commands could provide formating functionality.

The rendering on this website adds double brakets ⟦ ⟧ wherever units are used. You have to open the document to see how Maple renders.

NULL

Simple plot example: Solar irradiance in space

G__0 := 1361*Unit('W'/'m'^2)

1361*Units:-Unit(W/m^2)

(1)

G__0*sin(2*Pi*t/(24*Unit('h')))

1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h))

(2)

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'))

 

This plot has inconsistent axis labeling:

• 

The vertical axis has units but no name

• 

The horizontal axis has a name and units but they are not easily distinguishable. Misinterpretation is possible. Due to the close spacing the label could be read as a product of the dimension "time squared" (the time t times hours h is of the dimension time squared). Or the reader confounds name and units. (The use of italic fonts for names and roman fonts for units might not be noticeable and is a convention that is not used everywhere.)

 

The above labeling should be improved for communication, documentation or publication purposes.

 

A quick attempt using strings and the options useuints and labels.

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = ['d', kW/m^2], labels = ["Time t in days", "Exposure G in kV/m^2"])

 

Axes are now consistent and can be interpreted unambiguously. Formatting can still be improved.

 

Unfortunately, using the options useunits (for unit conversion) and labels this way introduces a new source of user error when labels are entered with the wrong units.

 

A way to address this and to ensure unit error-free plotting of expressions of physical quantities is the following:

 

Step1: Define two lists, one for the units to display and the other for the names to display

a := [Unit('s'), Unit('W'/'cm'^2)]; b := [t, G]

[t, G]

(3)

Step2: Compute labels from the lists

This step avoids the labeling error: No manual entry of units in labels required.

c := [b[1]/a[1], typeset(b[2]/a[2])]; d := [typeset(b[1], "  ", "⟦", a[1], "⟧"), typeset(b[2], "  ⟦", a[2], "⟧")]; e := [typeset(b[1], "  ", "(", a[1], ")"), typeset(b[2], "  (", a[2], ")")]

[typeset(t, "  ", "(", Units:-Unit(s), ")"), typeset(G, "  (", Units:-Unit(W/cm^2), ")")]

(4)

NULL

 

Dimensionless labels

 Double brackets

Parenthesis

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = c)

 

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = d)

 

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = e)

 

The axis values equal physical quantities divided by their units. The algebraic equation G*cm^2/W = 0.8e-1, for example, is physically speaking correct. Most functions of Maple can process dimensionless expression of the kind G*cm^2/W if G is given with appropriate physical units.

This way of using physical quantities is consistent with ISO 80000.  

Used in Maple to enter units in 2D-Math input mode

Can be confounded with functional notation. Units are therefore often written as a whole word (e.g. seconds instead of s).

 

 

NULL

The time to produce the above three plots was about 10 Minutes. The most part was spent to get the typesetting of the second and third plot correct.

 

What takes significant more time (more a question of hours when Typesetting is used for the first time) are

 

Labels with "/ cm^(2) "or 1/cm^2 formatting.

 

This formatting might be preferred but is unfortunately again not free from user errors. (I would probably use it if there was a simple and safe way).

f := [b[1]/a[1], b[2]/`#mrow(mo("W "),mo(" "),mo(" / "),msup(mo("cm"),mn("2")))`]; g := [typeset(b[1], "  ", "⟦", a[1], "⟧"), typeset(b[2], "  ⟦", (`@`(`@`(Units:-Unit, numer), op))(a[2]), "/", (`@`(`@`(Units:-Unit, denom), op))(a[2]), "⟧")]; h := [typeset(b[1], "  ", "(", Unit('s'), ")"), typeset(b[2], "  (", `#mrow(mo("W"),mo(" "),msup(mo("cm"),mn("-2")))`, ")")]

[typeset(t, "  ", "(", Units:-Unit(s), ")"), typeset(G, "  (", `#mrow(mo("W"),mo(" "),msup(mo("cm"),mn("-2")))`, ")")]

(5)

 

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = f)

 

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = g)

 

plot(1361*Units:-Unit(W/m^2)*sin((1/12)*Pi*t/Units:-Unit(h)), t = 0 .. 12*Unit('h'), useunits = a, labels = h)

 

NULL

 

 

 

Remarks

• 

For two reasons, I have not given an example with the often used square brackets [ ] because:
    
    Maple uses square brackets already for lists and indexing purposes,
    and ISO 80000 uses square brackets as an operator that extracts the unit from a physical quantity (e.g.       [G] = Unit('W'/'cm'^2)).

• 

Adding a unit to each value at axes ticks would definitely be a nice labeling feature for simple units.

• 

Programmatically analyzing the units defined in list a above and converting them in a generic way to a typesetting structure is not possible with a few high-level commands.

 

• 

For inline quotients like in 1/2, an additional backslash must be entered in 2D-Math: \/  

Unit('W')/Unit('cm')^2

Units:-Unit(W)/Units:-Unit(cm)^2

(6)

     This will not prevent evaluation to a normal quotient but the input can be used to create an atomic variable (select with mouse -> 2-D Math -> Atomic Variable)

`#mrow(mfenced(mi("W",fontstyle = "normal"),open = "⟦",close = "⟧"),mo("/"),mo("⁢"),msup(mfenced(mi("cm",fontstyle = "normal"),open = "⟦",close = "⟧"),mn("2")),mo("⁢"))`

`#mrow(mfenced(mi("W",fontstyle = "normal"),open = "⟦",close = "⟧"),mo("/"),mo("⁢"),msup(mfenced(mi("cm",fontstyle = "normal"),open = "⟦",close = "⟧"),mn("2")),mo("⁢"))`

(7)

     This makes labeling much easier as compared to typesetting commands (compare to the above statements).

f := [b[1]/a[1], b[2]/`#mrow(mfenced(mi("W",fontstyle = "normal"),open = "⟦",close = "⟧"),mo("/"),mo("⁢"),msup(mfenced(mi("cm",fontstyle = "normal"),open = "⟦",close = "⟧"),mn("2")),mo("⁢"))`]

[t/Units:-Unit(s), G/`#mrow(mfenced(mi("W",fontstyle = "normal"),open = "⟦",close = "⟧"),mo("/"),mo("⁢"),msup(mfenced(mi("cm",fontstyle = "normal"),open = "⟦",close = "⟧"),mn("2")),mo("⁢"))`]

(8)

In any case it is a good idea to read ?plot,typesetting before experimenting with typesetting.

 

Axes_with_unit_labels.mw

My personal preference is for dimensionless labels.

Note:

The solution to avoid labeling errors works only for Maple 2022 and higher.

Some plot commands do not support plotting with units, or they do not fully support it yet.

This is an successfull attempt to simulate space frames in MapleSim using the relatively new Rod component.

At t=3s, a lateral force component is applied to make the simulation more interesting.

The structure collapses/folds in an origami style fashion.

To build the model, MapleSim needs additional components.

For example, an equilateral triangle

requires the addition of rigid body frames at the connection of two rods.

Additionally, the rigid body frames must have initial position conditions (ICs) that match the intended structure.

Interestingly the ICs do not have to be set to Treat as Guess. It is only required to put an approximate coordinate. Leaving the ICs on Ignore was sufficient for the attached model.

Rod components can be replaced by Flexible Beam components which require considerably more simulation time and either Revolute joints at their ends or a rather complex connection with Rigid Body frames (of zero length) to adjacent Flexible Beam components.

 

Spaceframe_2.msim

 

The pendulum and the cantilever share simple-looking ordinary differential equations (ODEs), but they are challenging to solve:

This post derives solutions from Lawden and Bisshopp by Maple commands, which (to the best of my knowledge) have not been published providing not only results but also the accompanying computer algebra techniques. A tabulated format has been chosen to better highlight similarities and differences.

Both solutions have in common that in a first step, the unknown function is integrated and then in a second step the inverse of the unknown function (i.e., the independent variable) is integrated. Only in combination with a well-chosen set of initial/boundary conditions solutions are possible. This makes these two cases difficult to handle by generic integration methods.

Originally, I was not looking for this insight. I was more interested in an exact solution for nonlinear deformations to benchmark numerical simulation results.  Relatively quickly, I was able to achieve this with the help of this forum, but after that I was left with some nagging questions:

Why does Maple not provide a solution for the pendulum although one exists?

Why isn’t there an explicit solution for the cantilever when there is one for the pendulum?

Why is it so difficult to proof that elliptic expressions are equal?

Repeatedly, whenever there was time, I came back to these questions and got more and more a better understanding of the two problems and the overall context. It also required me to learn more of Maple, and I had to revisit fundamentals of functions, differential equations, and integration, which was entirely possible within Maples help system. Today, I am satisfied to the point that I think it is too much to expect Maple to provide a high-level general integration method for such problems.

I am also satisfied that I was able to combine all my findings scattered across many documents and Mapleprime questions into a single executable textbook-style document with hidden Maple code that:

Exclusively uses and manipulates equation expressions (no assignment operators := were required),

Avoids differentials that are often used in textbooks but (for good reasons) are not supported by Maple,

Exclusively applies high-level commands (i.e. no extraction of subexpression, manipulation
           and subsequent re-assembly of expressions was needed).

The solutions for the pendulum and the cantilever are substantially different although the ODEs and essential derivation steps are similar. I think that an explicit solution for the cantilever, as it exists for the pendulum, is impossible (using elliptic integrals and functions). I leave it open to comments: whether this is correct and whether it is attributable to the set of initial and boundary conditions, the different symmetry of the sine and cosine functions in the ODEs, or both. I hope that the tabular presentation will provide an easy overview, allowing to form an own opinion.

If you are patient enough to work through the table, you will find a link between the cantilever and the pendulum that you are probably not expecting. 

Finally, I have to give credit to Bisshopp, who was probably the first to provide a solution for the cantilever. The clarity and compactness on only 3 pages and the way how the inverse of functions was determined before the age of computers makes this paper worth studying. Also, Lawden has to be mentioned, who did the same on 3 pages for the pendulum in a marvelous book on elliptic functions and applications. It happens that he is overlooked in more recent publications and it’s unclear to me if he was the first who published an explicit solution. His book might be one of the last of its kind in this age of computers, and for that reason alone, it is worth enjoying as he enjoyed writing it.

 

The Pendulum and the Cantilever Side by Side

C_R, Summer 2023

• 

To better compare the pendulum and the cantilever, the symbol `ϕ` was chosen for the angle of the pendulum for the simple reason that this comparison started with bending theory, where `ϕ` is often used to denote a deflection angle.

• 

Leibniz and Newton notation was not used to make functional dependencies of variables visible. Instead functional notation  `ϕ` = `ϕ`(t) and `ϕ` = `ϕ`(s)is used.

• 

To create an executable document that maintains a clear representation, it is necessary to use functional notation for differential equations and remove functional notation for integration. To avoid using the same symbol for both the integration variable and the upper limit of integration, this document uses two ways to express when the upper limit of integration varies (i.e., depends on the dependent variable of the functions being searched, namely `ϕ` = `ϕ`(t) and `ϕ` = `ϕ`(s)). Both ways have their pros and cons.

 

Typesetting:-EnableTypesetRule({"EllipticE", "EllipticE2", "EllipticF", "EllipticK", "InverseJacobiAM", "InverseJacobiSN", "JacobiSN"})

 

©_®

Pendulum

Cantilever

Independent variable

Time t

Arclength s

Dependent variable

Angle of the pendulum with respect to direction of gravity `ϕ`(t)

varphi(t)

(1)

The slope of the cantilever with respect to the unbend state `ϕ`(s)

varphi(s)

(2)

Parameters

• 

Length l

• 

Gravitational constant g 

• 

Length L

• 

A force F at the free end

• 

A bending moment M at the free end

• 

The bending stiffness EI 

ODE

diff(varphi(t), t, t)+C*sin(varphi(t)) = 0

diff(diff(varphi(t), t), t)+C*sin(varphi(t)) = 0

(3)

(for derivation see for example Wikipedia [1] or Lawden [2])

diff(varphi(s), s, s)+C*cos(varphi(s)) = 0

diff(diff(varphi(s), s), s)+C*cos(varphi(s)) = 0

(4)

(for derivation see for example Bisshopp [3] or Beléndez [4])

Definitions

`ϕ`(0) = 0, `ϕ`((1/4)*T) = `ϕ__0`

varphi(0) = 0, varphi((1/4)*T) = varphi__0

(5)

`ϕ`(t)is periodic with the oscillation period T (i.e., the movement is bounded):

0 < abs(`&varphi;__0`) and abs(`&varphi;__0`) < Pi

0 < abs(varphi__0) and abs(varphi__0) < Pi

(6)

`&varphi;`(L) = `&varphi;__0`, Eval(diff(varphi(s), s), s = L) = 1/rho

varphi(L) = varphi__0, Eval(diff(varphi(s), s), s = L) = 1/rho

(7)

For a downward force:

0 < `&varphi;__0` and `&varphi;__0` < (1/2)*Pi

0 < varphi__0 and varphi__0 < (1/2)*Pi

(8)

Parameter C

"C=omega^(2)", where omegais the angular frequency of the pendulum for small anglular excursions

C = g/l

C = g/l

(9)

"Specific" Load

C = F/EI

C = F/EI

(10)

Initial/

boundary
conditions

Eval(varphi(t), t = -(1/4)*T) = -`&varphi;__0`, Eval(diff(varphi(t), t), t = -(1/4)*T) = 0

Eval(varphi(t), t = -(1/4)*T) = -varphi__0, Eval(diff(varphi(t), t), t = -(1/4)*T) = 0

(11)

Eval(varphi(s), s = 0) = 0, varphi(L) = varphi__0, Eval(diff(varphi(s), s), s = L) = 1/rho

Eval(varphi(s), s = 0) = 0, varphi(L) = varphi__0, Eval(diff(varphi(s), s), s = L) = 1/rho

(12)

Only the second condition is essential.
Additional essential condition: The length L of the cantilever beam is constant (not a boundary condition in its common sense but essential for the solution).

#1 integration step
with the second condition

Method: Integration with an integration factor (and converting to D notation, not shown)

DEtools:-intfactor(diff(diff(varphi(t), t), t)+C*sin(varphi(t)) = 0, `&varphi;`(t)); DETools:-firint((diff(diff(varphi(t), t), t)+C*sin(varphi(t)) = 0)*%, `&varphi;`(t))

-2*C*cos(varphi(t))+(diff(varphi(t), t))^2+c__1 = 0

(13)

Substituting initial conditions (11)

eval(convert(-2*C*cos(varphi(t))+(diff(varphi(t), t))^2+c__1 = 0, D), t = -(1/4)*T); convert(value({Eval(diff(varphi(t), t), t = -(1/4)*T) = 0, Eval(varphi(t), t = -(1/4)*T) = -varphi__0}), D); isolate(eval(`%%`, %), c__1)

c__1 = 2*C*cos(varphi__0)

(14)

and isolating diff(`&varphi;`(t), t)in (13) yields

convert(isolate(subs(c__1 = 2*C*cos(varphi__0), -2*C*cos(varphi(t))+(diff(varphi(t), t))^2+c__1 = 0), diff(`&varphi;`(t), t)), radical)

diff(varphi(t), t) = (2*C*cos(varphi(t))-2*C*cos(varphi__0))^(1/2)

(15)

Alternative method [5]: Integration to an implicit representation

dsolve({(Eval(varphi(s), s = 0) = 0, varphi(L) = varphi__0, Eval(diff(varphi(s), s), s = L) = 1/rho)[3], diff(diff(varphi(s), s), s)+C*cos(varphi(s)) = 0}, `&varphi;`(s), implicit)[1]

Int(1/(-2*C*sin(_a)+(2*C*sin(varphi(L))*rho^2+1)/rho^2)^(1/2), _a = 0 .. varphi(s))-s-c__2 = 0

(16)

and differentiation w.r.t. to s 

diff(Int(1/(-2*C*sin(_a)+(2*C*sin(varphi(L))*rho^2+1)/rho^2)^(1/2), _a = 0 .. varphi(s))-s-c__2 = 0, s); isolate(%, diff(`&varphi;`(s), s)); expand(subs(varphi(L) = varphi__0, Eval(diff(varphi(s), s), s = L) = 1/rho, %))

diff(varphi(s), s) = (-2*C*sin(varphi(s))+2*C*sin(varphi__0)+1/rho^2)^(1/2)

(17)

(This method works only if rho <> infinity; i.e., only with curvature/bending moment at the free end.)

Continuing now without bending momenteval(diff(varphi(s), s) = (-2*C*sin(varphi(s))+2*C*sin(varphi__0)+1/rho^2)^(1/2), rho = infinity)

diff(varphi(s), s) = (-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2)

(18)

#2 integrating the inverse function

With the chain rule

(diff(`&varphi;`(t), t))*(diff(t(`&varphi;`), `&varphi;`)) = 1

(diff(varphi(t), t))*(diff(t(varphi), varphi)) = 1

(19)

isolate((diff(varphi(t), t))*(diff(t(varphi), varphi)) = 1, diff(`&varphi;`(t), t))

diff(varphi(t), t) = 1/(diff(t(varphi), varphi))

(20)

(15) can be written as

isolate((diff(varphi(t), t) = 1/(diff(t(varphi), varphi)))*(1/(diff(varphi(t), t) = (2*C*cos(varphi(t))-2*C*cos(varphi__0))^(1/2))), diff(t(`&varphi;`), `&varphi;`)); subs(`&varphi;`(t) = `&varphi;`, %)

diff(t(varphi), varphi) = 1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2)

(21)

where the functional notation `&varphi;`(t)was removed for integration

map(Int, diff(t(varphi), varphi) = 1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), `&varphi;` = 0 .. `&varphi;`(t), continuous); (`@`(value, lhs) = rhs)(%)

-t(0)+t(varphi(t)) = Int(1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), varphi = 0 .. varphi(t), continuous)

(22)

This relation allows to determine the time t = t(`&varphi;`(t)) the pendulum takes to swing to a certain angle `&varphi;` = `&varphi;`*t__.

Similarly, with

(diff(`&varphi;`(s), s))*(diff(s(`&varphi;`), `&varphi;`)) = 1

(diff(varphi(s), s))*(diff(s(varphi), varphi)) = 1

(23)

isolate((diff(varphi(s), s))*(diff(s(varphi), varphi)) = 1, diff(`&varphi;`(s), s))

diff(varphi(s), s) = 1/(diff(s(varphi), varphi))

(24)

(16) can be written as

isolate((diff(varphi(s), s) = 1/(diff(s(varphi), varphi)))*(1/(diff(varphi(s), s) = (-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2))), diff(s(`&varphi;`), `&varphi;`))

diff(s(varphi), varphi) = 1/(-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2)

(25)

and integrated over s

map(int, diff(s(varphi), varphi) = 1/(-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2), `&varphi;` = 0 .. `&varphi;__s`, continuous)

-s(0)+s(varphi__s) = int(1/(-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous)

(26)

This expression relates the arclength s to the slope `&varphi;__s` at the location s (i.e., `&varphi;`(s) = `&varphi;__s`). It is the inverse of what we intend to solve (i.e., `&varphi;` = `&varphi;`(s)) but it is essential to derive a solution. Unlike for the pendulum, an indexed symbol `#msub(mi("&varphi;",fontstyle = "normal"),mi("s"))`has been chosen to avoid formally correct but uncommon expressions like "t(`&varphi;`(t))."

Special cases

Oscillation period T

When `&varphi;`*t__ = `&varphi;__0` the pendulum has swung a quater of the period. With
t(0) = 0, `&varphi;`(t) = `&varphi;__0`, t(`&varphi;__0`) = (1/4)*T

t(0) = 0, varphi(t) = varphi__0, t(varphi__0) = (1/4)*T

(27)

(22) becomes

subs(t(0) = 0, varphi(t) = varphi__0, t(varphi__0) = (1/4)*T, -t(0)+t(varphi(t)) = Int(1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), varphi = 0 .. varphi(t), continuous))

(1/4)*T = Int(1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), varphi = 0 .. varphi__0, continuous)

(28)

where `&varphi;`(t)was replaced in the integrant by `&varphi;` to create input that can be processes by the int command. After evaluation

isolate(`assuming`([simplify(value((1/4)*T = Int(1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), varphi = 0 .. varphi__0, continuous)))], [C > 0, 0 < `&varphi;__0` and `&varphi;__0` < Pi]), T)

T = 4*2^(1/2)*InverseJacobiAM((1/2)*varphi__0, 2^(1/2)/(1-cos(varphi__0))^(1/2))/(C^(1/2)*(1-cos(varphi__0))^(1/2))

(29)

where 1/am = InverseJacobiAM denotes the inverse Jacobian amplitude function.

Expression for `&varphi;__0`
Calculation of the length L in order to get an expression to determine the unknown slope `&varphi;__0` at the free end of the cantilever as a function of the load parameter C. With
s(0) = 0, s(`&varphi;__s`) = L, `&varphi;__s` = `&varphi;__0`, `&varphi;`(s) = `&varphi;`

s(0) = 0, s(varphi__s) = L, varphi__s = varphi__0, varphi(s) = varphi

(30)

(26) becomessubs(s(0) = 0, s(varphi__s) = L, varphi__s = varphi__0, varphi(s) = varphi, -s(0)+s(varphi__s) = int(1/(-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous))

L = int(1/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__0, continuous)

(31)

where `&varphi;`(s)is replaced by `&varphi;` to create input that can be processes by the int command. After evaluation

`assuming`([value(L = int(1/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__0, continuous))], [C > 0, 0 < `&varphi;__0` and `&varphi;__0` < (1/2)*Pi])

L = EllipticK((1/2)*(2*sin(varphi__0)+2)^(1/2))/C^(1/2)-EllipticF(1/(sin(varphi__0)+1)^(1/2), (1/2)*(2*sin(varphi__0)+2)^(1/2))/C^(1/2)

(32)

where K = EllipticK and F = EllipticF denote the complete and incomplete elliptic integrals of the first kind.

#3 Solutions of particular interest  

Explicit solution for `&varphi;`(t) in bounded motion.
Rearranging (22)

`assuming`([simplify(2*C*(-t(0)+t(varphi(t)) = Int(1/(2*C*cos(varphi)-2*C*cos(varphi__0))^(1/2), varphi = 0 .. varphi(t), continuous)))], [C > 0])/(sqrt(C)*sqrt(2))

-C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = Int(1/(cos(varphi)-cos(varphi__0))^(1/2), varphi = 0 .. varphi(t), continuous)

(33)

and substituting this essential identity (expression 5.1.3 from Lawden [2])

-cos(`&varphi;__0`)+cos(`&varphi;`) = 2*sin((1/2)*`&varphi;__0`)^2-2*sin((1/2)*`&varphi;`)^2

cos(varphi)-cos(varphi__0) = 2*sin((1/2)*varphi__0)^2-2*sin((1/2)*varphi)^2

(34)

yields

subs(cos(varphi)-cos(varphi__0) = 2*sin((1/2)*varphi__0)^2-2*sin((1/2)*varphi)^2, -C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = Int(1/(cos(varphi)-cos(varphi__0))^(1/2), varphi = 0 .. varphi(t), continuous))

-C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = Int(1/(2*sin((1/2)*varphi__0)^2-2*sin((1/2)*varphi)^2)^(1/2), varphi = 0 .. varphi(t), continuous)

(35)

`assuming`([value(-C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = Int(1/(2*sin((1/2)*varphi__0)^2-2*sin((1/2)*varphi)^2)^(1/2), varphi = 0 .. varphi(t), continuous))], [0 < `&varphi;__0` and `&varphi;__0` < Pi])

-C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = 2^(1/2)*InverseJacobiAM((1/2)*varphi(t), csc((1/2)*varphi__0))/sin((1/2)*varphi__0)

(36)

which simplifies further with

t(0) = 0, t(`&varphi;`(t)) = t

t(0) = 0, t(varphi(t)) = t

(37)

to

subs(t(0) = 0, t(varphi(t)) = t, ((-C^(1/2)*2^(1/2)*(t(0)-t(varphi(t))) = 2^(1/2)*InverseJacobiAM((1/2)*varphi(t), csc((1/2)*varphi__0))/sin((1/2)*varphi__0))*(1/sqrt(2)))*sin((1/2)*`&varphi;__0`))

sin((1/2)*varphi__0)*C^(1/2)*t = InverseJacobiAM((1/2)*varphi(t), csc((1/2)*varphi__0))

(38)

Mapping now sn = JacobiSN with the modulus csc((1/2)*`&varphi;__0`) to the expression above (Maple converts InverseJacobiAM to InverseJacobiSN and simplifies automatically)

map(JacobiSN, sin((1/2)*varphi__0)*C^(1/2)*t = InverseJacobiAM((1/2)*varphi(t), csc((1/2)*varphi__0)), csc((1/2)*`&varphi;__0`))

JacobiSN(sin((1/2)*varphi__0)*C^(1/2)*t, csc((1/2)*varphi__0)) = sin((1/2)*varphi(t))

(39)

the angle `&varphi;`as a function of time is obtained explicitly

convert(isolate(JacobiSN(sin((1/2)*varphi__0)*C^(1/2)*t, csc((1/2)*varphi__0)) = sin((1/2)*varphi(t)), `&varphi;`(t)), sincos)

varphi(t) = 2*arcsin(JacobiSN(sin((1/2)*varphi__0)*C^(1/2)*t, 1/sin((1/2)*varphi__0)))

(40)

Bending curve of the cantilever for a given load (i.e., for a given `&varphi;__0`, obtainable from (32)).

To obtain a parametric form x(p), y(p)of the bending curve, the following two ODEs have to be integrated

diff(x(s), s) = cos(`&varphi;`(s)), diff(y(s), s) = sin(`&varphi;`(s))

diff(x(s), s) = cos(varphi(s)), diff(y(s), s) = sin(varphi(s))

(41)

Applying

diff(x(s), s) = (diff(x(`&varphi;`), `&varphi;`))*(diff(`&varphi;`(s), s)), diff(y(s), s) = (diff(y(`&varphi;`), `&varphi;`))*(diff(`&varphi;`(s), s))

diff(x(s), s) = (diff(x(varphi), varphi))*(diff(varphi(s), s)), diff(y(s), s) = (diff(y(varphi), varphi))*(diff(varphi(s), s))

(42)

in the following way to (41)

subs({diff(x(s), s) = (diff(x(varphi), varphi))*(diff(varphi(s), s)), diff(y(s), s) = (diff(y(varphi), varphi))*(diff(varphi(s), s))}, diff(varphi(s), s) = 1/(diff(s(varphi), varphi)), {diff(x(s), s) = cos(varphi(s)), diff(y(s), s) = sin(varphi(s))})[]

(diff(x(varphi), varphi))/(diff(s(varphi), varphi)) = cos(varphi(s)), (diff(y(varphi), varphi))/(diff(s(varphi), varphi)) = sin(varphi(s))

(43)

isolate(((diff(x(varphi), varphi))/(diff(s(varphi), varphi)) = cos(varphi(s)), (diff(y(varphi), varphi))/(diff(s(varphi), varphi)) = sin(varphi(s)))[1], diff(x(`&varphi;`), `&varphi;`)), isolate(((diff(x(varphi), varphi))/(diff(s(varphi), varphi)) = cos(varphi(s)), (diff(y(varphi), varphi))/(diff(s(varphi), varphi)) = sin(varphi(s)))[2], diff(y(`&varphi;`), `&varphi;`))

diff(x(varphi), varphi) = cos(varphi(s))*(diff(s(varphi), varphi)), diff(y(varphi), varphi) = sin(varphi(s))*(diff(s(varphi), varphi))

(44)

yields two ODEs where x and y depend on the variable `&varphi;`(s). Substituting (25) and `&varphi;`(s) = `&varphi;` 

subs(diff(s(varphi), varphi) = 1/(-2*C*sin(varphi(s))+2*C*sin(varphi__0))^(1/2), `&varphi;`(s) = `&varphi;`, [diff(x(varphi), varphi) = cos(varphi(s))*(diff(s(varphi), varphi)), diff(y(varphi), varphi) = sin(varphi(s))*(diff(s(varphi), varphi))])[]

diff(x(varphi), varphi) = cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), diff(y(varphi), varphi) = sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2)

(45)

map(Int, (diff(x(varphi), varphi) = cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), diff(y(varphi), varphi) = sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2))[1], `&varphi;` = 0 .. `&varphi;__s`, continuous), map(Int, (diff(x(varphi), varphi) = cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), diff(y(varphi), varphi) = sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2))[2], `&varphi;` = 0 .. `&varphi;__s`, continuous)

Int(diff(x(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous), Int(diff(y(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous)

(46)

results in a parametric solution with parameter `#mrow(mi("p"),mo("&equals;"),mi("\`&varphi;__s\`"))` where "0<= `&varphi;__s`<=`&varphi;__0`." 
For the x coordinate

subs(`assuming`([x(0) = 0, simplify(value((Int(diff(x(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous), Int(diff(y(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous))[1]))], [0 < `&varphi;` and `&varphi;` < (1/2)*Pi, C > 0]))

x(varphi__s) = (2^(1/2)*sin(varphi__0)^(1/2)-(-2*sin(varphi__s)+2*sin(varphi__0))^(1/2))/C^(1/2)

(47)

For the y coordinate a long expression with the following structure

subs(`assuming`([y(0) = 0, simplify(value((Int(diff(x(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous), Int(diff(y(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous))[2]), radical)], [0 < `&varphi;` and `&varphi;` < (1/2)*Pi, 0 < `&varphi;__0` and `&varphi;__0` < (1/2)*Pi, 0 < `&varphi;__s` and `&varphi;__s` < (1/2)*Pi, C > 0])); y(`&varphi;__s`) = (A__1*EllipticE(z__1, k)+A__2*EllipticF(z__1, k)+A__3*EllipticE(z__2, k)+A__4*EllipticF(z__2, k))/sqrt(C)

y(varphi__s) = (A__1*EllipticE(z__1, k)+A__2*EllipticF(z__1, k)+A__3*EllipticE(z__2, k)+A__4*EllipticF(z__2, k))/C^(1/2)

(48)

is obtained where A__i = A__i(`&varphi;__s`, `&varphi;__0`), z__i = z__i(`&varphi;__s`, `&varphi;__0`) and k = k(`&varphi;__0`).

NULL

Not required in the above: To derive an explicit solution, Lawden performed a change of variable of this kind

sin((1/2)*`&varphi;`) = sin((1/2)*`&varphi;__0`)*sin(u)

sin((1/2)*varphi) = sin((1/2)*varphi__0)*sin(u)

(49)

 

which is not needed with Maple commands.

 

Furthermore: Formally, it would have been nice to have the pendulum start its movement at t=0 at an angle -`&varphi;__0`. However, this leads to an output in (36) with two inverse elliptic functions where `&varphi;`(t) is difficult if not impossible to isolate.

 

Solution for the free end of the cantilever (i.e., `&varphi;__s` = `&varphi;__0`and s = L)

subs(x(0) = 0, `&varphi;__s` = `&varphi;__0`, x(`&varphi;__0`) = x(L), x(varphi__s) = (2^(1/2)*sin(varphi__0)^(1/2)-(-2*sin(varphi__s)+2*sin(varphi__0))^(1/2))/C^(1/2))

x(L) = 2^(1/2)*sin(varphi__0)^(1/2)/C^(1/2)

(50)

subs(`assuming`([y(0) = 0, y(`&varphi;__0`) = y(L), simplify(value(subs(`&varphi;__s` = `&varphi;__0`, (Int(diff(x(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(cos(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous), Int(diff(y(varphi), varphi), varphi = 0 .. varphi__s, continuous) = Int(sin(varphi)/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__s, continuous))[2])))], [0 < `&varphi;` and `&varphi;` < (1/2)*Pi, 0 < `&varphi;__0` and `&varphi;__0` < (1/2)*Pi, 0 < `&varphi;__s` and `&varphi;__s` < (1/2)*Pi, C > 0]))

y(L) = (EllipticK((1/2)*(2*sin(varphi__0)+2)^(1/2))-EllipticF(1/(sin(varphi__0)+1)^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2))-2*EllipticE((1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2))+2*EllipticE(1/(sin(varphi__0)+1)^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2)))/C^(1/2)

(51)

``

Some remarks

NULL

NULL

 

Comparison to the solution from Lawden

Expression (39) compared to sin((1/2)*theta(t)) = sin((1/2)*alpha)*sn(t+K, sin((1/2)*alpha))

sin((1/2)*theta(t)) = sin((1/2)*alpha)*sn(t+K, sin((1/2)*alpha))

(52)

(pasted from DLMF 22.19.2 for the case C = 1). The same with adapted variables  subs(alpha = `&varphi;__0`, sn = JacobiSN, K = 0, theta(t) = `&varphi;__t`, sin((1/2)*theta(t)) = sin((1/2)*alpha)*sn(t+K, sin((1/2)*alpha)))

sin((1/2)*varphi__t) = sin((1/2)*varphi__0)*JacobiSN(t, sin((1/2)*varphi__0))

(53)

Now equating the left hand side of (39) to JacobiSN(z, k)and identifying the parameters z and k

subs(C = 1, lhs(JacobiSN(sin((1/2)*varphi__0)*C^(1/2)*t, csc((1/2)*varphi__0)) = sin((1/2)*varphi(t)))) = JacobiSN(z, k)

JacobiSN(sin((1/2)*varphi__0)*t, csc((1/2)*varphi__0)) = JacobiSN(z, k)

(54)

map(op, JacobiSN(sin((1/2)*varphi__0)*t, csc((1/2)*varphi__0)) = JacobiSN(z, k)); solve([(rhs-lhs)(%)], {k, z})[]

k = csc((1/2)*varphi__0), z = sin((1/2)*varphi__0)*t

(55)

Using the following identity from Maple's FunctionAdvisor and the correspondence in (55)

FunctionAdvisor(identities, JacobiSN(z, 1/k))[5]

JacobiSN(z, k) = JacobiSN(z*k, 1/k)/k

(56)

yields

convert(subs(k = csc((1/2)*varphi__0), z = sin((1/2)*varphi__0)*t, JacobiSN(z, k) = JacobiSN(z*k, 1/k)/k), sincos)

JacobiSN(sin((1/2)*varphi__0)*t, 1/sin((1/2)*varphi__0)) = sin((1/2)*varphi__0)*JacobiSN(t, sin((1/2)*varphi__0))

(57)

Comparing this with (53) confirms that (40) is correct.

Equivalent expressions to determine `&varphi;__0` 

As returned by value:

normal(L = EllipticK((1/2)*(2*sin(varphi__0)+2)^(1/2))/C^(1/2)-EllipticF(1/(sin(varphi__0)+1)^(1/2), (1/2)*(2*sin(varphi__0)+2)^(1/2))/C^(1/2))

L = (EllipticK((1/2)*(2*sin(varphi__0)+2)^(1/2))-EllipticF(1/(sin(varphi__0)+1)^(1/2), (1/2)*(2*sin(varphi__0)+2)^(1/2)))/C^(1/2)

(58)

simplify instead of value:

convert(`assuming`([simplify(L = int(1/(-2*C*sin(varphi)+2*C*sin(varphi__0))^(1/2), varphi = 0 .. varphi__0, continuous))], [0 < `&varphi;__0` and `&varphi;__0` < (1/2)*Pi]), sincos)

L = -I*2^(1/2)*EllipticF(I*sin(varphi__0)^(1/2)/(1-sin(varphi__0))^(1/2), I*(1-sin(varphi__0))^(1/2)/(sin(varphi__0)+1)^(1/2))/(C^(1/2)*(sin(varphi__0)+1)^(1/2))

(59)

With integration tools and change of variables using x = sin(`&varphi;`):

Int(1/sqrt(-2*C*sin(`&varphi;`)+2*C*sin(`&varphi;__0`)), `&varphi;` = 0 .. `&varphi;__0`); L = IntegrationTools:-Change(%, x = sin(`&varphi;`), x); simplify(subs(isolate(x__0 = sin(`&varphi;__0`), `&varphi;__0`), %)); subs(x__0 = sin(`&varphi;__0`), `assuming`([value(%)], [0 < x and x < 1, 0 < x__0 and x__0 < 1]))

L = 2^(1/2)*EllipticF(sin(varphi__0)^(1/2)/(sin(varphi__0)+1)^(1/2), I*(-sin(varphi__0)^2+1)^(1/2)/(-1+sin(varphi__0)))/(C^(1/2)*(1-sin(varphi__0))^(1/2))

(60)

Without having a Maple way: Christian Wolinski has provided 3 additional expressions where one is of particular simplicity [6]

  L = EllipticF(sqrt(1-1/(sin(`&varphi;__0`)+1))*sqrt(2), (1/2)*sqrt(2)*sqrt(sin(`&varphi;__0`)+1))/sqrt(C)

L = EllipticF((1-1/(sin(varphi__0)+1))^(1/2)*2^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2))/C^(1/2)

(61)

 

 

Useful identities

FunctionAdvisor(identities, InverseJacobiSN)[3]

InverseJacobiSN(z, k) = InverseJacobiAM(arcsin(z), k)

(62)

InverseJacobiSN(z, k) = InverseJacobiAM(arcsin(z), k)

InverseJacobiSN(z, k) = EllipticF(z, k)

(63)

FunctionAdvisor(identities, InverseJacobiSN)[1]

JacobiSN(InverseJacobiSN(z, k), k) = z

(64)

FunctionAdvisor(identities, JacobiSN)[5]

JacobiSN(z, 1/k) = k*JacobiSN(z/k, k)

(65)

 

Explicit solution for `&varphi;__0`

Since Maple's solve does not solve (32) for `&varphi;__0`, one could try to isolate `&varphi;__0` in expression (32) by combining "somehow" the two elliptic expression into a single expression and to apply an inverse operation to it.

Maple's simplify or combine do not seem to be able to help in this respect. There might be addition theorems that could be applied but identifying them in Maple or in DLMF requires expertise in this field of special functions.

Easier is to try to evaluate (31) in different ways (see above) and hope for success. This yielded equivalent expressions with only one elliptic integral EllipticF.

Using identities such elliptic integrals can be converted to inverse elliptic functions where elliptic functions can be applied to. Trying this exemplarily for (61)  

 

L/sqrt(C) = InverseJacobiSN(sqrt(1-1/(sin(`&varphi;__0`)+1))*sqrt(2), (1/2)*sqrt(2)*sqrt(sin(`&varphi;__0`)+1))

L/C^(1/2) = InverseJacobiSN((1-1/(sin(varphi__0)+1))^(1/2)*2^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2))

(66)

map(JacobiSN, L/C^(1/2) = InverseJacobiSN((1-1/(sin(varphi__0)+1))^(1/2)*2^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2)), (1/2)*sqrt(2)*sqrt(sin(`&varphi;__0`)+1))

JacobiSN(L/C^(1/2), (1/2)*2^(1/2)*(sin(varphi__0)+1)^(1/2)) = (1-1/(sin(varphi__0)+1))^(1/2)*2^(1/2)

(67)

does not isolate `&varphi;__0`in the same way `&varphi;`(t)could be separated for the pendulum. The reason why such attempts are deemed to fail is simple. Differently to the pendulum, `#msub(mi("&varphi;",fontstyle = "normal"),mi("0"))` is not a constant "input" to the system but, causally speaking, an "output". While C in the case of the pendulum is constant and independent of `&varphi;__0`, C and `&varphi;__0`functionally depend on each other for the cantilever; i.e., `&varphi;__0` = `&varphi;__0`(C). This fundamentally makes the two cases different although many derivation steps are similar.

Student:-ODEs

Student:-ODEs:-ODESteps(diff(diff(varphi(t), t), t)+C*sin(varphi(t)) = 0)

Error, (in Student:-ODEs:-OdeSolveOrder2) ODE is not supported

 

Student:-ODEs:-ODESteps(diff(diff(varphi(s), s), s)+C*cos(varphi(s)) = 0)gives a solution for `&varphi;` = `&varphi;`(s)with two integration constants C1 and C2, but determining the integration constants with the first two boundary conditions of (12) is probably impossible. odetest confirms that the solution is correct but one of the arguments of an elliptic function is not unitfree (which raises more questions): JacobiSN((1/2)*sqrt(2*C-2*C1)*(-s+C2), sqrt(-(C1+C)/(-C1+C))).

Links

https://www.mapleprimes.com/questions/232863-Testing-Maples-Solution-Of-The-Nonlinear

https://www.mapleprimes.com/questions/131520-Animation-Of-A-Simple-Pendulum

 

Applying a horizontal load instead of a vertical bends the cantilever in an arc-like fashion. For this load case the corresponding ODE is that of the pendulum (see [2], chapter 5, exercise 8). The parametric equation (bending curve) of the arc becomes simpler but still no explicit solution `&varphi;__0` = `&varphi;__0`(C)seems possible.

Alternative symbols for 4

θ

phi, theta

References

[1] https://en.wikipedia.org/wiki/Pendulum_(mechanics)
[2] Lawden, Derek F. “Elliptic Functions and Applications.” Acta Applicandae Mathematica 24 (1989): 201-202.
[3] Bisshopp, K.E. and Drucker, D.C. (1945) Large Deflection of Cantilever Beams. Quarterly of Applied Mathematics, 3, 272-275.
[4] BELÉNDEZ, Tarsicio; NEIPP, Cristian; BELÉNDEZ, Augusto. "Large and small deflections of a cantilever beam". European Journal of Physics. Vol. 23, No. 3 (May 2002). ISSN 0143-0807, pp. 371-379
[5] Rouben Rostamian, https://www.mapleprimes.com/questions/236511-How-To-Integrate-This-Ode-And-How-To#answer295192
[6] Christian Wolinski, https://www.mapleprimes.com/questions/233304-How-To-Find-The-Inverse-Function-Of#comment283638

 

 

Download Cantilever_and_pendulum_side_by_side_-_input_hidden.mw

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

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

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

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

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

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

 

 

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

 

 

Some observations on the attached model:

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

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

Ball_bouncing_on_a_flexible_beam.msim

 

 

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

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

Given transfer function

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

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

(1)

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

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

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

(2)

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

true

(3)

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

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

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

(4)

simplify(Units[Unit](m*kg/(s^2*A))/Units[Unit](kg^2*m^2/(s^5*A^2)) = Units[Unit](m/V))

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

(5)

"Simplification" attempts with the denominator

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

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

(6)

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

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

(7)

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

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

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

(8)

Expand does not expand the argument of Unit commands.

expand(Units[Unit](kg^2*m^2/(s^2*A^2))); lprint(%)

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

 

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

 

NULL

C1: Expanding Unit command

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

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

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

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

(9)

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

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

(10)

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

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

(11)

C2: Using the Natural Units Environment

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

NULL

C3: A type declaration or unit assumptions on names

A type declaration as an option of commands like in

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

true

(12)

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

`assuming`([simplify(H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2))))], [s::(Unit(1/s))]); `assuming`([combine(H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2))), 'units')], [s::(Unit(1/s))])

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

 

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

 

C4: DynamicSystems Package with units

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

 

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

 

Download Collecting_and_expanding_units.mw

1 2 3 4 5 6 Page 3 of 6