The Fictitious Force Awakens

Inverse Dynamics: Prologue B

Falk Mielke, University of Antwerp, July 3, 2020

forces dynamics calculations

You can download this notebook here.

This is one part of a series of blog posts on Inverse Dynamics.

import numpy as NP # numerics
import sympy as SYM # symbolic operations
import sympy.physics.mechanics as MECH # some physics/mechanics tools
import matplotlib.pyplot as MPP # plotting

import WrenchToolbox as WT # contains an improved "Wrench" object

# printing of formulas
SYM.init_printing(use_latex='mathjax', pretty_print = False)
MECH.init_vprinting(use_latex='mathjax', pretty_print = False)

# to disable a current sympy/lambdify/numpy hiccup
from warnings import filterwarnings
filterwarnings('ignore', category=NP.VisibleDeprecationWarning)

# shorthand to get the basis of a reference frame
GetCoordinates = lambda rf: [rf[idx] for idx in rf.indices]

# transforming a sympy vector to a matrix and back
MatrixToVector = lambda mat, rf: sum([elm * base for elm, base in zip(mat, [rf.x, rf.y, rf.z]) ])

The custom Wrench Toolbox can be downloaded here.

The Fictitious Force Awakens

In the previous notebook, I introduced my notation to solve inverse dynamics for systems of rigid bodies.

Based on a simple example, I demonstrated the general procedure, calculated the equations involved in Balance of Forces and Balance of Moments, and introduced Wrenches to naturally combine both balances. These balances directly emerge from the conservation of linear and angular momentum, so they can be regarded as a physical first principle.

I closed the example by pointing at a mismatch between the expected Coriolis Force (zero) and the calculation according to a naïve application of the Wikipedia formula. I will continue this here by calculating Balances in the rotating reference frame.

The main purpose of this notebook series is to see how fictitious forces enter the balance equations in inverse dynamics. My idea is to use the information about the inertial properties of the system under study to dissect the net joint force and moment into (i) biological components and (ii) components that stem from passive, fictitious forces. I want to know whether my piglets are on the dark side (greedily exploiting fictitious forces to reduce energy demand of their movements) or on the light side (resisting the fictitious force). [update 2020/09/12] As I understood with the fourth example, the fictitious forces are not forces, but have a different character, so the initial question formulated here turned out to be irrelevant.

There are three fictitious forces that arise from reference frame rotation: the Euler force, the Centrifugal force, and the Coriolis force. [update 2020/09/12] There is a fourth one due to linear acceleration of a reference frame, sometimes called D’Alambert Force, which I initially missed although it might be the one most often used to introduce fictitious forces in physics courses (cars, elevators).

Starting in the next notebook with the Euler force, I intend to construct elementary examples in which one can see the fictitious forces in isolation.

Fictitious forces only emerge when executing the balance equations in a non-static reference frame. Although “fictitious” sounds a bit like “fiction”, the consequences of fictitious forces are real. They can be felt and will cause trouble, not only when you try to walk drunk over a merry-go-round. I will elaborate more on this assessment later.

The balance equations I previously demonstrated were all in the static reference frame \({S}\). So to get to the fictitious forces, we’ll have to repeat the procedure in the “body” frame \({B}\).

The following line will reload most of the preparation so that we do not have to repeat it.

# %load

Recall the problem:

This time, we’ll move into the body frame \({B}\). But to be able to do that, we’ll revisit and extend the kinematics.

Kinematics Ambiguity

We’ll use shortcuts for selecting which derivative of position vectors applies.

But these are only indices to lists that are assembled below.

# shorthands for the derivative levels
p = 0 # zero'th derivative: the position vector
v = 1 # first derivative: the velocity
a = 2 # second derivative: the acceleration

Linear Kinematics

The situation is more complex than before, where we only had to use the three derivatives for a single \(x\) and \(\omega\). Instead of constructing longer and longer variable names for many variables (e.g. eightteen of the kind x_S_joint_ddot), I’ll give this some structure by using multi-level dictionaries and lists.

# prepare an empty dictionary
x = {}

# start with the position vectors
x['S'] = {pt: [pointmass.pos_from(refpoint).express(world).to_matrix(world)] \
          for pt, refpoint in [('O', origin), ('J', joint), ('P', pointmass)] \
x['B'] = {pt: [pointmass.pos_from(refpoint).express(body).to_matrix(body)] \
          for pt, refpoint in [('O', origin), ('J', joint), ('P', pointmass)] \

# automatically get time derivatives, in loops
# loop reference frames
for frame in x.keys():
    # loop reference points
    for refpoint in x[frame].keys():
        # loop differentials
        for diff_nr in range(2): 

One can access these structures easily, for example by writing x['B']['O'][p], which gives the point mass (x) position vector (p) relative to the origin (O) in the body frame (B).

# printout of position vectors in different frames
[0  1.0  0]
[0  1.0  0]

Note that the position vectors are equal in both reference frames, except for the mandatory rotation (transformation). But look at their derivatives:

# in contrast, the time derivatives depend on reference frame
[0  0  0]
[-1.0  0  0]

There is a first interesting observation [1]:

Hence, there is an ambiguity in how to get kinematics vectors in the body frame: one can use zero vectors which are correct “per definition”, or one could take the equally correct vectors from the static frame and transform them.

The latter one, for distinction, will be denoted by a tilde:

\( \tilde x_{B} = R_{BS}x_{S}\)

with \(R_{BS}\) being the transformation matrix. The time derivatives of \(\tilde x_{B}\) are \(\dot{\tilde x}_{B}\) and \(\ddot{\tilde x}_{B}\), accordingly.

Angular Kinematics

To a certain extent, this is analogous to linear kinematics. However, angular kinematics do not require a reference point. If the frame turns, it turns, and all points fixed to the frame turn the same way (i.e. if they don’t move relative to each other).

# the same dictionary magic as above
ω = {}

# angular positions
ω['S'] = [SYM.Matrix([[0],[0],[φ]])]
ω['B'] = [SYM.Matrix([[0],[0],[0]])]

# time derivatives
for frame in ω.keys():
    for diff_nr in range(2): 

# quick check: which of these matches the settings we entered for our reference frames:
SYM.pprint(PlugAllIn(ω['S'][v] - body.ang_vel_in(world).to_matrix(world)).T)
SYM.pprint(PlugAllIn(ω['B'][v] - body.ang_vel_in(body).to_matrix(body)).T)
# answer: both!
[0  0  0]
[0  0  0]

Now we can also check for the same ambiguity as above:

# In contrast, the time derivatives depend on reference frame
[0  0  0]
[0  0  1.0]

Second observation [2]:

As with linear measures, I keep \(\tilde \omega_B\) , \(\dot{\tilde\omega}_B \) and \(\ddot{\tilde\omega}_B\) for distinction:

\(\tilde \omega_{B} = R_{BS}\omega_{S}\)

The \(\omega_B\) (the zero one) might be completely artificial. Most textbooks I have seen just write \( \omega_S \) (and actually mean: \( \dot{\tilde\omega}B = R\{BS}\dot\omega_S \) ) for the kinematic equations.

But which ones of the two should we take? And is the choice identical to that for the linear motion parameters?

Chosing the Correct Vectors

Let’s birefly check with the Coriolis Force, as defined by Wikipedia (but with explicit dots for derivatives):

\(F_{cor} = -2m \dot{\omega} \times \dot{x}\)

We can use \(\dot{{x}}_{B}\) and \(\dot{{\omega}}_{B}\):

F_cor = -2*m*(ω['B'][v]).cross(x['B']['O'][v])

\(\left[\begin{matrix}0 & 0 & 0\end{matrix}\right]\)

Or we use \(\dot{\tilde x}_{B}\) and \(\dot{\tilde \omega}_{B}\):

F_cor = -2*m*((trafo['B']['S'] * ω['S'][v]).cross(trafo['B']['S'] * x['S']['O'][v]))

\(\left[\begin{matrix}0 & 2.0 & 0\end{matrix}\right]\)

There are a third and fourth option: the “mixed ones”.

One with \(\dot{{x}}_{B}\) but \(\dot{\tilde \omega}_{B}\)

F_cor = -2*m*(trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])

\(\left[\begin{matrix}0 & 0 & 0\end{matrix}\right]\)

And yet another with \(\dot{\tilde x}_{B}\) and \(\dot{{\omega}}_{B}\)

F_cor = -2*m*(ω['B'][v]).cross(trafo['B']['S'] * x['S']['O'][v])

\(\left[\begin{matrix}0 & 0 & 0\end{matrix}\right]\)

Looks like we can rule out the second option. However, some of the others are certainly zero by coincidence, because our example is simple.

This shows that there is definitely more to work through. Spoiler: number three (the first mixed strategy, \(\dot{{x}}_{B}\) and \(\dot{\tilde \omega}_{B}\)) is most likely the correct option. But I will only demonstrate that three notebooks ahead.

I focused here on Coriolis because it is the decisive one: Euler- and Centrifugal Force are indifferent to the observation [1] above because they only require a position vector for \(P\).

But there is a general way of how to derive fictitious forces, which others have explained better than I could.

The Five-Term Acceleration Equation

Obviously, the faulty Coriolis force gave me a trouble. I did some more web research and found this notebook by Frank Owen, who runs a company offering “engineering consultancy”.

That document derives the accelerations which appear in a rotating reference frame with sufficient mathematical rigor. He takes a position of an arbitrary point on a rotating, not necessarily rigid body. The author then deconstructs the position vector into

Deriving that combined vector in all its elements with the product rule gives all the fictitious accelerations. Because that procedure is sufficiently rigorous, it is obvious that there can’t be any other fictitious forces (answering a question that I’ve been asking myself early one, why are these three or four fictitious forces supposed to be all there are?; more discussion in the fourth notebook on D’Alambert force). The outcome of the maths are vectorial formulas that can be easily applied in the sympy setting we established above. There are five components:

(1) The acceleration of the rigid body COM (abbreviated C, with distance \(OC\) from the origin) where the reference frame is attached.

\(a_{A} = a_{com, S} = \frac{d^2(OC)}{dt^2}\)

(2) The movement of a focus point in that reference frame.

\(a_{xy} = a_{com, B} = = \frac{d^2(CP)}{dt^2}\)

(where \(P\) is a point different from the COM; this term is zero in our simplified example)

(3) Coriolis acceleration

\(a_{C} = -2m\cdot (R\_{BS} \dot{\omega}_{S}) \times \dot{x}_{S} = -2m\dot{\tilde \omega}_{B} \times \dot{x}_{B}\)

(4) tangential acceleration (corresponding to Euler Force)

\(a_{t} = -m\cdot (R_{BS} \ddot{\omega}_{S}) \times x_{S} = -m\ddot{\tilde \omega}_{B} \times {x}_{B}\)

(5) normal acceleration (corresponding to Centrifugal Force)

\(a_{n} = -m\cdot (R_{BS} \dot{\omega}_{S}) \times \left( (R_{BS} \dot{\omega}_{S}) \times \dot{x}_{S} \right) = -m\dot{\tilde \omega}_{B} \times \left( \dot{\tilde \omega}_{B} \times \dot{{x}}_{B} \right)\)

There is more explanation worth reading on the particularly non-intuitive Coriolis force in a separate notebook.

What is the difference between Owen’s formulas and the Wikipedia one?

The difference is that Owen keeps track of all the reference frames and definitions. From his document, it is clear which “\(r\)”, “\(x\)” and “\(\omega\)” he uses. He also adds accurate sketches that visualize the outcome.

I wrote about these two different possible \(x_B\) [1] and \(\omega_B\) [2] above: one that is zero due to the definition of the reference frame, and one that is the transformed \(\tilde \omega_{B} = R_{BS} \omega_{S}\). Seeing them in the acceleration formulas now, you might be tempted to forget about \(\omega_{B}\), for example, and use \(\tilde \omega_{B}\) all over the place (yes, I actually did that mistake, and it cost me hours to find out). However, we actually need both: recall that, for the Balance Equation, we have to use \(x_B\) and \(\omega_{B}\), and not the transformed \(R_{BS} x_{S}\)/\(R_{BS} \omega_{S}\).

Take home message: Make sure, in this and all that follows, in what reference frame parameters are defined and whether they were transformed. And oh, yes, never trust Wikipedia.

[update 2020/09/12] Although I must admit I should have read Wikipedia (and Owen) more carefully from the start. I completely missed out on term (1) from above, the acceleration of the attachment point.

Overview: Fictitious Forces

We’ll assume for now (and try to disprove) that the third option from above is correct. Also, because the choice should be consistent for all three fictitious forces, I will apply the same logic to the others as well.

This gives:

The Euler force


\(F_e = -m \ddot{\tilde\omega}_{B} \times x_B\)

F_e_B = -m*(trafo['B']['S'] * ω['S'][a]).cross(x['B']['J'][p])

\(\left[\begin{matrix}0 & 0 & 0\end{matrix}\right]\)

the Centrifugal Force:

\(F_{cf} = -m\dot{\tilde \omega}_{B} \times \left(\dot{\tilde \omega}_{B} \times {x}_{B} \right)\)

F_cf_B = -m*(trafo['B']['S'] * ω['S'][v]).cross((trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][p]))

\(\left[\begin{matrix}0 & 1.0 & 0\end{matrix}\right]\)

the Coriolis Force: (never forget the factor two!)

\(F_{cor} = -2m\dot{\tilde \omega}_{B} \times \dot{{x}}_{B}\)

F_cor_B = -2*m*(trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])

\(\left[\begin{matrix}0 & 0 & 0\end{matrix}\right]\)

These are the values we expect to see for our initial example. And for those of my readers who noticed: indeed, I intentionally haven’t discussed the reference point for linear motion parameters here, because 'J' and 'O' coincide. I’ll save that for a later notebook.

Balance Of Wrenches

In the Inertial Frame

We repeat the calculation of Balance of Wrenches for the inertial frame, because (if correct) that result must be correct, independent of fictitious forces. It is always a good starting point.

Because the balance equations are based on a first principle, they must give the same result, no matter which reference frame was used for calculation. We can and should compare later.

(1) Assemble all external wrenches that change the momentum of the rigid body. Here, this will only be the joint wrench we are interested in.

force_components_S = SYM.symbols('f_{JS1:4}', real = True)
moment_components_S = SYM.symbols('m_{JS1:4}', real = True)

joint_wrench_S = WT.Wrench.FromComponents(world, joint, force_components_S, moment_components_S)

[f_{JS1}  f_{JS2}  f_{JS3}  m_{JS1}  m_{JS2}  m_{JS3}]

(2) Collect the dynamic wrench.

dynamic_force_S = m*x['S']['O'][a]
dynamic_moment_S = I_S['P']*ω['S'][a] + ω['S'][v].cross(I_S['P']*ω['S'][v])

# assembling the dynamic wrench
dynamic_wrench_S = WT.Wrench.FromMatrices(world, pointmass, dynamic_force_S, dynamic_moment_S \
# display the outcome

\(\left[\begin{matrix}sin\left(t\right) & - cos\left(t\right) & 0 & 0 & 0 & 0\end{matrix}\right]\)

(3) Solve the balance equations for the unknown joint wrench.

# get the equations (one per component)
equations_S = dynamic_wrench_S.Equate(joint_wrench_S)

# desired outcome variables
components_S =  [*force_components_S, *moment_components_S]

# extract the solution
solutions = {}
for param, sol in SYM.solve(equations_S \
                           , components_S \
    solutions[param] = sol

# transform the result into a wrench
W_JS = WT.Wrench.FromMatrix(world, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_S]) )
[sin(t)  -cos(t)  0  0  0  0]

Same result as before. Keep it stored, e.g. by writing it on your hand or forehead.

Naïve Balance in the Body Frame

The procedure is the same in the body frame. As mentioned somewhere, always double check that you get the reference points and reference frames right!

As you will see, it’s not enough to change all ‘S’ to ‘B’, the outcome will not match.

force_components_B = SYM.symbols('f_{JB1:4}', real = True)
moment_components_B = SYM.symbols('m_{JB1:4}', real = True)

joint_wrench_B = WT.Wrench.FromComponents(body, joint, force_components_B, moment_components_B)

[f_{JB1}  f_{JB2}  f_{JB3}  m_{JB1}  m_{JB2}  m_{JB3}]
dynamic_force_B = m*x['B']['O'][a]
dynamic_moment_B = I_B['P']*ω['B'][a] + ω['B'][v].cross(I_B['P']*ω['B'][v])

# assembling the dynamic wrench
dynamic_wrench_B = WT.Wrench.FromMatrices(body, pointmass, dynamic_force_B, dynamic_moment_B \
# display the outcome

\(\left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]\)

# get the equations (one per component)
equations_B = dynamic_wrench_B.Equate(joint_wrench_B)

# desired outcome variables
components_B =  [*force_components_B, *moment_components_B]

# extract the solution
solutions = {}
for param, sol in SYM.solve(equations_B \
                           , components_B \
    solutions[param] = sol

# transform the result into a wrench
W_JB = WT.Wrench.FromMatrix(body, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_B]) )
[0  0  0  0  0  0]

Our Centrifugal Wrench is gone.

Adding Fictitious Wrenches

Recall the definition of fictitious forces: they have to be introduced when going to the non-static reference frame in order for the balance equations to be correct.

Here you go. This is exactly what we have to do. Of course we use wrenches, and as elaborated above we use certain variants of the kinematic elements (namely \(\tilde \omega_{B}\) and \(x_{B}\)).

Here comes Euler, the “tangential” one (which actually means “not normal”, as he certainly was):

euler_wrench = WT.Wrench(body, pointmass \
                    , MatrixToVector(-m*((trafo['B']['S'] * ω['S'][a]).cross(x['B']['O'][p])), body) \

[0  0  0  0  0  0]

Here comes the normal component, also called centrifugal, so hold your seat:

centrifugal_wrench = WT.Wrench(body, pointmass \
                    , MatrixToVector(-m*( \
                             (trafo['B']['S'] * ω['S'][v]).cross( \
                                                 (trafo['B']['S'] * ω['S'][v]).cross(x['B']['O'][p]) \
                                            ) \
                                        ), body) \

[0  1.0  0  0  0  0]

And here the one that everyone remembers from, but noone understands at school, Coriolis:

coriolis_wrench = WT.Wrench(body, pointmass \
                    , MatrixToVector(-2*m*((trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])), body) \

[0  0  0  0  0  0]

Note a little computational trick that I used above in initializing the wrenches:

But pay attention: defining wrenches by only the forces requires you to set the correct point of application and the correct reference frame. That is easier said than done.

In real problems, one will be less interested in the ones above. We are actually after the joint wrench:

force_components_B = SYM.symbols('f_{JB1:4}', real = True)
moment_components_B = SYM.symbols('m_{JB1:4}', real = True)

joint_wrench_B = WT.Wrench.FromComponents(body, joint, force_components_B, moment_components_B)

[f_{JB1}  f_{JB2}  f_{JB3}  m_{JB1}  m_{JB2}  m_{JB3}]

Then there is the left handside with the dynamic wrench:

dynamic_force_B = m*x['B']['O'][a]

dynamic_moment_B = I_B['P']*ω['B'][a] + ω['B'][v].cross(I_B['P']*ω['B'][v])

# assembling the dynamic wrench
dynamic_wrench_B = WT.Wrench.FromMatrices(body, pointmass, dynamic_force_B, dynamic_moment_B \


\(\left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]\)

And, finally, the balance equations:

# balance equations
equations_B = dynamic_wrench_B.Equate(joint_wrench_B + euler_wrench + centrifugal_wrench + coriolis_wrench)

# the usual solving procedure
components_B =  [*force_components_B, *moment_components_B]

solutions = {}
for param, sol in SYM.solve(equations_B \
                           , components_B \
    solutions[param] = SYM.simplify(sol)

# and here comes the body frame joint wrench
W_JB = WT.Wrench.FromMatrix(body, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_B]) )
[0  -1.0  0  0  0  0]

Still different from the result above … … because we need to transform!


\(\left[\begin{matrix}sin\left(t\right) & - cos\left(t\right) & 0 & 0 & 0 & 0\end{matrix}\right]\)


\(\left[\begin{matrix}sin\left(t\right) & - cos\left(t\right) & 0 & 0 & 0 & 0\end{matrix}\right]\)

# because the "PlugAllIn" could cheat us, here is the ultimate check:
PlugAllIn(( - W_JB ).Matrix()).T

\(\left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]\)



I have encountered several types of definitions of fictitious forces during my initial attempts.

The first type is a practical one, which goes along the lines of “fictitious forces arise because an observer who is captured in the rotating frame is blind to the movement of the frame”. Yes, that is a valid description, but a disappointing definition: What else might the observer be blind to? How can we calculate something we are blind for? And why can’t I not always switch perspective in cases where I’m in charge of the data?

The second type of definition is a technical one: “fictitious forces need to be introduced to make balance equations work in the rotating reference frame”. That is also a valid point, yet it leaves similar questions as above: Which ones do I have to include when? Where do the formulas come from?

I have hinted at these aspects above. They certainly help to develop an intuition about fictitious forces, but for me they were not satisfactory. You can find these two categories scattered all over the Wikipedia articles about fictitious forces. Despite my rants herein, I’m not disappointed about Wikipedia, it is just sad that this confirmed my prior on that “Encyclopedia”. But I also see myself as a source of the problem: whereas for trained physicists it is often obvious which and how vectors add up, my humble, untrained mind needs some more explicit guidance. [update 2020/09/12] And although I’m using myself in a tautologous way here, those of you who read my first version of the notebook are witness that I could confidently ignore the D’Alambert force without noticing.

I found a third definition most convincing, yet it is more like a prescription than a definition. It is the one which I found in the workbooks by Frank Owen (five-term acceleration, see references).

[update 2020/09/12] Finally, after having finished my work on this series of notebooks, I would like to attempt another definition, which I think captures the inverse dynamics standpoint: Instead of being forces, Fictitious Forces represent a change of the momentum which is not captured by kinematics in the observer’s frame of reference. Instead of writing them on the side of the external forces (of which I think as the cause of a momentum change), I tend to put them on the “effect” side and would rather call them “unbalanced change of momentum” (“unbalanced” as in “not part of the balance”).

I had to write a couple of notebooks to come to that conclusion, and would be glad if you like to follow through. The fourth one will explain this fourth definition in more detail.

This preparational notebook has solved some puzzles which I encountered when starting to learn about fictitious forces. Make sure to keep track of your reference frames and reference points. Good (i.e. explicit, structured, well commented) programming style can facilitate this. And remember that, even if you think it is obvious which parameters you took in which transformation, make sure that others can reproduce your line of thought.

With these tools at our hands, we can advance to elementary examples. I’ll start by isolating Euler, i.e. with an example that should show the Euler force in isolation, in the next notebook. If you don’t mind the explanations, feel free to skip ahead to the code from the application to a general (n-link) limb.