C_R

3717 Reputation

21 Badges

6 years, 197 days

MaplePrimes Activity


These are Posts that have been published by C_R

From a discussion about expanding unit expressions with compound units I concluded that expanding derived units such as Newton, Watt, Volt, Tesla,... to SI base units is difficult in Maple.

Unintentionally, I came across a rather simple solution for SI units.

toSIbu := x -> x = Units:-Unit(simplify(x/Unit('kg'))*Unit('kg'));

converts derived SI units to SI base units. It’s the inverse of what the units packages and simplify do (i.e. simplification to derived units).

What makes it maybe more interesting: It also works, again unintentionally, on other units than SI units. If, one day, you come along an erg or a hartree or or a kyne and you cannot guess the SI units convert/units needs, try

toSIbu(Unit('pound'));
toSIbu(Unit('hp'));
toSIbu(Unit('electron'));
toSIbu(Unit('hartree'));
toSIbu(Unit('bohr'));
toSIbu(Unit('barye'));
toSIbu(Unit('kyne'));
toSIbu(Unit('erg'));
toSIbu(Unit(mile/gal(petroleum)));

Maybe handy one day when you do not trust AI or the web.


 

NULL 

toSIbu := x -> x = Units:-Unit(simplify(x/Unit('kg'))*Unit('kg')):
toSIbu(Unit('N'));
toSIbu(Unit('J'));
toSIbu(Unit('W'));
toSIbu(Unit('Pa'));
toSIbu(Unit('C'));
toSIbu(Unit('F'));
toSIbu(Unit('S'));
toSIbu(Unit('H'));
toSIbu(Unit('T'));
toSIbu(Unit('V'));
toSIbu(Unit('Wb'));
toSIbu(Unit('Omega'));
toSIbu(Unit('lx'));
toSIbu(Unit('lm'));
toSIbu(Unit('degC'));
toSIbu(Unit('rad'));
toSIbu(Unit('sr'));

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

 

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

 

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

 

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

 

Units:-Unit(C) = Units:-Unit(A*s)

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

Units:-Unit(lx) = Units:-Unit(cd/m^2)

 

Units:-Unit(lm) = Units:-Unit(cd)

 

Units:-Unit(`°C`) = Units:-Unit(K)

 

Units:-Unit(rad) = Units:-Unit(m/m(radius))

 

Units:-Unit(sr) = Units:-Unit(m^2/m(radius)^2)

(1)

NULL


 

Download toSIbu.mw


(All done with Maple 2024 without loading any package)

 

 

 

This post summarizes links for those who have not studied numerical integration methods from scratch and are interested in simulation settings in MapleSim (like me).

The MapleSim help pages simulation settings and advanced simulation settings give first guidance for the trained user but do not provide explanations or links for the terms used in the description of the settings (as for example: stiffness, constraint stabilization, constraint projection, events and event iteration,...).

It can easily be overlooked that Maple help pages provide further information for most of the terms. Under the assumption that MapleSim uses the same terminology as Maple, I recommend to first have a look at Maple help topics before consulting the web or other resources. Since searching and retrieving can be time consuming, I made a list of helpful links.

There are still some open points. I would be happy for more links and help in filling these gaps.

 

How Maple simulates

?MapleSimUserGuide,Chapter04:
section 4.1 How MapleSim Simulates a Model

?tasks,generatingCode

Ein Bild, das Text, Screenshot, Diagramm, Design enthält.

Automatisch generierte Beschreibung

 

Solvers

An overview of solvers: ?dsolve,numeric

Differential Algebraic Equation introduction: ?MaplePortal,DAE

Overview of numeric differential-algebraic equation solvers (index reduction, constraint drift, projection):
 ?examples/numeric_DAE and ?dsolve,numeric,DAE_extension

Stiffness and stiff solvers

Stiffness and stiff IVPs: ?dsolve,Stiffness

Events

?dsolve,numeric,Events

Time events and state events

Event handling:

?MapleSimUserGuide,Chapter04:
section 4.1 How MapleSim Simulates a Model

Event iteration:

?MapleSimUserGuide,Chapter05:
section 5.5 Selecting the Code Generation Options

Iteration, hysteresis, Intermediate steps: ?tasks,generatingCode

Hysteresis:

Hysteresis in value or also in time?

Do variable solvers adapt the value of event hysteresis during runtime?

 

Baumgarte constraint stabilization, unconstrained dynamics, constrained dynamics

?MapleSim,Multibody,Dynamic_Exports
(in combination with ?MapleSim,Multibody,Kinematic_Exports)

?examples/numeric_DAE

?tasks,generatingCode

?MapleSimUserGuide,Chapter05:
section 5.5 Selecting the Code Generation Options

Error control

              ?dsolve,numeric,Error_Control

              Absolute error: ?dsolve,numeric,IVP

              Relative error: (relative to what?)

Index1 error control and Index1 Tollerance: see solvers

Scaling

scalemethod (this does not seem to exist in Maple)

 

Examples (Multibody)

Events

                            Catapult
                            (from MapleSim>Help>Examples>Physical Domains>Multibody)
                            contact events

                          Catapult_-_Events.msim

                            Throwing a ball
                            (from MapleSim>Help>Examples>Physical Domains>Multibody)

                            conditional events (with boolean logic)

                          Throwing_a_Ball_-_Events.msim

              Solvers

              Conservation of energy of a pendulum depends on solvers.
                           Euler increases energy, implict Euler dissipates energy.

             Pendulum_for_solver_comparision.msim

           

Constraint dirft/projection

              2-d rigid slider crank

               (from MapleSim>Help>Examples>Physical Domains>Multibody)

              projection off leads to assembly desintegration after 2000 s simulation

             2D_Rigid_Slider_Crank_-_constraint_projection.msim

                         A stiff solver improves constraint drift, but only delays desintegration

                         Baumgarte constraint stabilization prevents simulation error but shows dislocated rigid body frames

 

It can happen when an operation is interrupted by  that Maple does not return to  and still shows .

This can give the false impression that the Maple server in charge of the evaluation did not get the message to stop whatever it was doing.

By giving Maple an impossible task to solve analytically

f1 := x1 - x1*sin(x1 + 5*x2) - x2*cos(5*x1 - x2);
f2 := x2 - x2*sin(5*x1 - 3*x2) + x1*cos(3*x1 + 5*x2);
solve({f1, f2});

I have noticed in the Windows Task Manager that freeing allocated memory can take much longer than one might think.

In one case it took 30 minutes to free 24 Gb of total allocated memory (21 Gb of it in RAM/physical memory). In this case the interrupt button became active (turned from grey to red ) two times and memory continued piling up  again.

Lessons learned for me:

  • The task manager is not only a valuable indicator for task activity but also for the interruption/memory freeing process.
  • Before killing a whole Maple session and potentially losing the last state of a worksheet it can pay off to wait and repeatedly interrupt an operation.

 

Suggestion: When the maple server gets an interrupt request, it could report to the GUI that it is in an interruption state and is no longer evaluating input. For example changing the message in the status bar from Evaluating... to Interrupting...

A ball on a turntable can move in circles instead of falling off the edge (provided the initial conditions are appropriate). The effect was demonstrated in a video and can be simulated with MapleSim. The amination below shows a simulation of a frictionless case (falling off the table) and the case with a coefficient of friction of one.

Also demonstrated in the video: Tilting the table leads to a sideward (not a downhill) movement of the ball.

The presenter of the video noted that in the untilted state, the ball eventually drifts off the table, attributing this to slippage. This drift is also observable in the animation above, where the ball starts moving in a spiral, whereas in a Maple simulation (below to the left), the ball follows a perfect circle. Only after optimizing contact and initial conditions, MapleSim produced a result (using the same parameters) that exhibits a similar circle, with a slight difference in angular orientation after completing two revolutions about the center of the circle.

 

Some observations on the MapleSim model:

  • The results only slightly depend on the solvers. Numerical inaccuracies do not seem to be the reason for the difference in orientation. (Edit: see update below for the reason).
  • The ball bounces up and down in the MapleSim simulation (0.0025 of the balls radius). The bouncing is caused by the fact that the initial position of the ball is above the elastic equilibrium position with respect to the table (the elastic contact makes the ball sink in a bit). Adding damping in the settings of the contact element attenuates this effect and reduces the drift.
  • Drift is not observable anymore if in the contact element setting for "kmu" (smoothness coefficient of sliding friction) is set to larger values (above 10 in this example). This is an indication that sliding friction occurs during the simulation.
  • Choosing the equilibrium position as initial condition for the ball does not prevent initial bouncing because MapleSim sets an initial velocity for the ball that is directed away from the table. I did not manage to enforce strictly zero velocity. Maybe someone can tell why that is and how to set MapleSim to start the simulation without vertical velocity. (Edit: see update below for the reason).
  • Assuming an initial velocity towards the contact to cancel the initial vertical velocity set by MapleSim substantially reduces bouncing and increases the diameter of the circle. This finally leads to a diameter that matches the Maple simulation. Therefore the initial bouncing combined with slippage seems to dissipate energy which leads to smaller circles. Depending on the contact settings and initial conditions for vertical movement the diameter of the circle varied moderately by about 15%.

In summary, MapleSim can be parametrized to simulate an ideal case without slippage. From there it should be possible to tune contact parameters to closely reproduce experiments where parameters are often not well known in advance.  

Some thoughts for future enhancement of MapleSim:

  • In the model presented here, a patterned ball would have been helpful to visualize the tumbling movement of the ball. Marking two opposite sides of the ball with colored smaller spheres is a fiddly exercise that does not look nice.
  • A sensor component that calculates the energy of a moving rigid body would have helped identifying energy loss of the system. Ideally the component could have two ports for the rotational and translational energy components. I see professional interest in such calculations and not only educational value for toy experiments.
  • A port for slippage on the contact elements would have helped understanding where slippage occurs. Where slippage is, there is wear and this is also of interest for industrial applications.

Turntable_Paradox.msim (contains parameter sets for the above observations)

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.

1 2 3 4 5 6 Page 3 of 6