I'm trying to simulate a simple electric motor from scratch. In order to do this I need to be able to apply forces to a body which can only travel around a circular parth, a rotating magnet. If my research is correct I now need something called constraints. I found this article which in the beginning describes pretty much exactly what I want to do:


However I am quite confused by the formular for the "Lagrange multiplier" as it doesnt enforce any kind of radius. I also tried implementing the constraint force like the article describes but the movement of my magnet now just seems random to me.

Please let me know what I'm missing. If you couldn't already tell by now I'm pretty new to this stuff and there very well could be a much easier solution to my problem.


enter image description here

This is my magnet which is supposed to only be able to rotate around the origin. I've found two ways to make this possible.

  1. To restrict the DOF of it and to only calculate the rotational acceleration produced by the acting forces. This is of course much easier and I think kind of hacky which is not really what I want.

  2. I introduce some kind of constraint force which keeps the position of the magnet fixed.

During my simulation loop I'm adding all the acting forces to my acceleration (a = F/m) to then numerically integrate in order to estimate the next position of my object. I had hoped this constraint force would just be another force I'd be able to add to my acceleration.

  • $\begingroup$ FYI - Search youtube for Lagrange multipliers as there are some really nice explanatory videos on how they enforce constraints. $\endgroup$ Feb 3, 2023 at 16:08
  • $\begingroup$ A nice read for you also here, although it mostly deals with contacts and not constraints. cs.cmu.edu/~baraff/sigcourse/notesd2.pdf $\endgroup$ Feb 3, 2023 at 17:00

2 Answers 2


$\def \b {\mathbf}$

General approach

Newton Equations of Motion \begin{align*} & m\,\ddot{\mathbf{R}}=\b F_a+\b F_c\tag 1 \end{align*} and the constraint equation e.g. for circular path \begin{align*} & x^2+y^2=r^2\quad,2\,x\,\dot{x}+2\,y\dot{y}=0\quad, \underbrace{\begin{bmatrix} x & y \\ \end{bmatrix}}_{\b C_c}\begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}=0 \end{align*} with $~\b F_c=\b C_c^T\,\lambda~$ equation (1) \begin{align*} & m\,\ddot{\mathbf{R}}=\b F_a+\b C_c^T\,\lambda~\tag 2 \end{align*} choose $~x~$ to be the generalized coordinate , thus from the constraint eqaution you obtain

\begin{align*} &y=\sqrt{r^2-x^2}\quad\Rightarrow \dot{y}=-\frac{x}{y}\dot{x}\\\\ &\dot{\mathbf{R}}=\begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}= \underbrace{\begin{bmatrix} 1 \\ -\frac{x}{y} \\ \end{bmatrix}}_{\b J}\,\dot{x}\quad \Rightarrow \ddot{\mathbf{R}}=\b J\,\ddot{x}+\frac{\partial \dot{\mathbf{R}}}{\partial x}\,\dot{x}= \b J\,\ddot{x}+\frac{\partial \left(\b J\,\dot{x}\right)}{\partial x}\,\dot{x} \end{align*} thus equation (2) \begin{align*} & m\,\b J\,\ddot{x}=\b F_a-\underbrace{m\,\frac{\partial \left(\b J\,\dot{x}\right)}{\partial x}\,\dot{x}}_{\b F_z}+\b C_c^T\,\lambda~\tag 3 \end{align*} multiply Eq. (3) from the left with $~\b C_c~$ you obtain \begin{align*} & \left[\b C_c\,\b C_c^T\right]\,\lambda=\b C_c\left(\b F_z-\b F_a\right)\tag 4 \end{align*} from here you get the generalized force $~\lambda~$ , and the constraint forces $~\b F_c=\b C_c^T\,\lambda~$

to obtain the equation of motions , multiply Eq. (3) from the left with $~\b J^T~$ \begin{align*} & m\,\b J^T\,\b J\,\ddot{x}=\b J^T\,\left(\b F_a-\b F_z\right)\tag 5 \end{align*}

notice that $~\b C_c\,\b J=\b 0$

  • $~\b R~$ Position vector
  • $~\b F_a~$ Applied force
  • $~\b F_c~$ Constraint force
  • $~\b \lambda~$ Generalized constraint force


form the time derivative of constraint equation you obtain the constraint matrix $~\b C_c~$.

from the choose of the generalized coordinate $~(x~)~$ and the velocity equation $~\dot{\b{R}}~$ you obtain the Jacobian-Matrix $~\b J~$.

from here go to polar coordinate \begin{align*} &x=r\,\cos(\varphi)\quad,\dot{x}=-r\,\sin(\varphi)\\ &y=r\,\sin(\varphi)\quad,\dot{y}=+r\,\cos(\varphi)\\ \end{align*}

your generalized coordinate is now $~\varphi~$ thus $~\b C_c=\b C_c(\varphi)~$ and $\b J=\b J(\varphi)~$ with \begin{align*} &\b F_z=m\,\frac{\partial \left(\b J\,\dot{\varphi}\right)}{\partial \varphi}\,\dot{\varphi} \end{align*}

you can solve Eq. (4) for $~\lambda~$ and obtain the constraint force $\b F_c=\b C_c^T\,\lambda$

\begin{align*} \b F_c=\left[ \begin {array}{c} - \left( \cos \left( \varphi \right) {\it F_x}+\sin \left( \varphi \right) {\it Fy}+m\,r{\dot\varphi }^{2} \right) \cos \left( \varphi \right) \\ - \left( \cos \left( \varphi \right) {\it F_x}+\sin \left( \varphi \right) {\it F_y }+m\,r{\dot\varphi }^{2} \right) \sin \left( \varphi \right) \end {array} \right]\\ \end{align*}

with Eq. (5) you obtain the equations of motion \begin{align*} &m\,r^2\,\ddot\varphi=r\,(F_y\cos(\varphi)-F_y\sin(\varphi)) \end{align*}

The constraint force components $~(x,y)~$ system \begin{align*} &\b F_c= \left[ \begin {array}{c} -{\frac { \left( {\it F_x}\,x+{\it F_y}\,y+m{{ \dot x}}^{2}+m{{\dot y}}^{2} \right) x}{{r}^{2}}}\\ -{\frac { \left( {\it F_x}\,x+{\it F_y}\,y+m{{\dot x}}^{2}+m{{\dot y}}^{ 2} \right) y}{{r}^{2}}}\end {array} \right] \end{align*}

and the EQM's

\begin{align*} &\begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \end{bmatrix}=\b F_c \end{align*}

the initial conditions must fulfilled the constraint equation $~x_0^2+y_0^2=r^2~$

The simulation with $~r=1~,F_x=F_y=1,x_0=0.3~,y_0=\sqrt{1-0.3^2}$

enter image description here

enter image description here

  • $\begingroup$ Thank you so much for your very detailed answer :) I'm still going through your explanation trying to understand as much as possible. Could you clarify what "Cc" and "CTc" are? Also is this constraint force now the force that pushes the body back on its path because when I tried implementing it the movement again just seemed random and glitchy? $\endgroup$
    – 3m1l
    Feb 3, 2023 at 23:38
  • $\begingroup$ $~C_c~$ is constraint matrix $~C_c^T~$ is the matrix transpose. $\endgroup$
    – Eli
    Feb 4, 2023 at 8:29
  • $\begingroup$ I will put you more information for the implementation $\endgroup$
    – Eli
    Feb 4, 2023 at 9:41
  • $\begingroup$ Again thank you for your patience with me :) It is still not really working as intended. Just to be clear I treat Fc like any other force and add it to my acceleration at the end right? Also for mass equal to 1 Fx and Fy are equal to x'' and y'' right? $\endgroup$
    – 3m1l
    Feb 4, 2023 at 22:45
  • $\begingroup$ Not clear what you are looking for. If you are looking for the equation of motion, you don’t need the constraint forces at all. $\endgroup$
    – Eli
    Feb 5, 2023 at 11:26

There is an alternative way to the excellent answer @Eli provided.

Instead of using $(x,y,z)$ coordinates to track the position of each magnet, write the position in terms of coordinates that describe the motion. In this case of each magnet is constrained in the xy plane, but at a fixed radius $r$ from the center then you write

$$ \text{position vector} \\ \boldsymbol{p}(\theta) = \pmatrix{ r \cos \theta \\ r \sin \theta \\ 0} $$

where your position vector $\boldsymbol{p}$ is described in terms of ONE variable only, the angle $\theta$. This is called a generalized coordinate.

Now apply the derivative chain rule with fixed $r$ and varying $\theta$ to get velocity and acceleration

$$ \text{velocity vector} \\ \boldsymbol{v}(\theta, \dot{\theta}) = \pmatrix{ -r \dot{\theta} \sin \theta \\ r \dot{\theta} \cos \theta \\ 0} $$

$$ \text{acceleration vector} \\ \boldsymbol{a}(\theta, \dot{\theta}, \ddot{\theta}) = \pmatrix{ -r \ddot{\theta} \sin\theta - r \dot{\theta}^2 \cos\theta \\ r \ddot{\theta}\cos\theta - r \dot{\theta}^2 \sin \theta \\ 0 } $$

Now you decompose the velocity vector into two parts, based on your generalized coordinate speed, as so

$$ \boldsymbol{v} = \boldsymbol{J}\, \dot{\theta} = \pmatrix{ -r \sin \theta \\ r \cos \theta \\ 0} \dot{\theta} $$

where $\boldsymbol{J}$ is the jacobian which describes how much each velocity component changes with each change in the generalized coordinate.

Finally, a constraint force $\boldsymbol{F}_c$ must do no work in the system, so the power at any instant should be zero, $P = \boldsymbol{F}_c \cdot \boldsymbol{v} = 0$. This sets up the constraint equations for the forces

$$ \boldsymbol{F}_c \cdot \boldsymbol{J} = 0 $$

where $\cdot$ is the dot product, and in linear algebra world it is equal to $\boldsymbol{J}^\top \boldsymbol{F}_c =0$ where ${}^\top$ is the transpose operator.

So your system of equations are

$$\begin{gathered} \boldsymbol{F}_{\rm ext} + \boldsymbol{F}_c = m \boldsymbol{a} \\ \boldsymbol{J}^\top \boldsymbol{F}_c = 0 \end{gathered} $$

where $\boldsymbol{F}_{\rm ext}$ is any applied external forces.

In terms of components the above is

$$ \begin{aligned} Fx_c + Fx_{\rm ext} & = - m r \dot{\theta}^2 \cos \theta - m r \ddot{\theta} \sin \theta \\ Fy_c + Fy_{\rm ext} & = - m r \dot{\theta}^2 \sin \theta + m r \ddot{\theta} \cos \theta \\ Fz_c + Fz_{\rm ext} & = 0 \\ Fy_c r \cos \theta - Fx_c r \sin \theta & = 0 \end{aligned} $$

The system is to be solved in terms of the constraint forces $Fx_c$, $Fy_c$, $Fz_c$ and motion $\ddot{\theta}$.

The above process can easily be generalized for multiple bodies with multiple degrees of freedom, all formulated from your choice of generalized coordinates to track all the centers of masses.

$$ \boldsymbol{q} = \pmatrix{ \vdots \\ \theta_i \\ \vdots} $$

$$ \boldsymbol{v} = \pmatrix{ \vdots \\ \boldsymbol{v}_i \\ \vdots} = \mathbf{J}\, \boldsymbol{ \dot{q} } $$

$$ \boldsymbol{F}_c = \pmatrix{ \vdots \\ {\boldsymbol{F}_c}_i \\ \vdots} $$

$$ \boldsymbol{a} = \mathbf{J}\,\ddot{q} + \mathbf{\dot{J}} \dot{q} $$ $$ \mathbf{M} = \begin{bmatrix} \ddots & & \\ & m_i & \\ & & \ddots \end{bmatrix}$$ $$ \boldsymbol{F}_{\rm ext} + \boldsymbol{F}_c = \mathbf{M}\,\boldsymbol{a}$$ $$ \mathbf{J}^\top \boldsymbol{F}_c = \boldsymbol{0}$$

  • $\begingroup$ The constraint force in your approach is depending on the generalized acceleration $~\ddot{\mathbf{q}}~$ this means that if you want to simulate friction (that depending on the constraint force) ,you get "algebraic loop" in the equations , which cause numerical simulation problem. $\endgroup$
    – Eli
    Feb 4, 2023 at 14:34
  • $\begingroup$ other way is : from the Jacobian matrix $~\mathbf J(\mathbf q)~$ obtain the constraint matrix $~\mathbf C~$ by solving the matrix equation $~\mathbf C\,\mathbf J=\mathbf 0~$ the solution $~\mathbf C(\mathbf q)~$ is not uniqe but with need jest one solution. from here you can use "my" equation (4). for this example $~\mathbf C=[1,\tan(\varphi)]~$ is one solution that fulfilled the equation $\endgroup$
    – Eli
    Feb 4, 2023 at 16:17
  • $\begingroup$ to make it clear . with different solution for $~\mathbf C~$ we obtain deferent generalized forces $~\mathbf \lambda~$ but the constraint force $~\mathbf F_c~$ is always the same $\endgroup$
    – Eli
    Feb 4, 2023 at 16:28
  • 1
    $\begingroup$ @Eli - friction can be a problem indeed. But there are ways to include its effect as the constraint does work. You will end up with a DAE system in the end, regardless of the approach if friction is included. Since sliding friction depends on normal forces, you need to include an algebraic expression between forces (similar to the Painleve problem). $\endgroup$ Feb 4, 2023 at 21:30
  • $\begingroup$ John ,this is your DAE ? \begin{align*} \begin{bmatrix} \mathbf J & -\mathbf I_3\, \\ \mathbf 0 & \mathbf J^T \,\\ \end{bmatrix}\, \begin{bmatrix} \ddot{\mathbf{q}} \\ \mathbf F_c \\ \end{bmatrix}= \begin{bmatrix} \mathbf F_a-\mathbf F_z \\ \mathbf 0\\ \end{bmatrix} \end{align*} the solution for $~\mathbf F_c~$ is equal zero. $\endgroup$
    – Eli
    Feb 5, 2023 at 11:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.