# Numerical torque calculation

Suppose I can compute interaction energy of two rigid bodies as a function of their coordinates of centers of masses and Euler rotation angles (total 6 + 6 degrees of freedom). Now I can numerically compute force acting on the center of mass of the body by calculating numerical derivatives e.g. $F_x = (E(x + dx) - E(x - dx)) / (2 * dx)$. But if you do the same for Euler angles this doesn't give you torques. So how do I convert numerical derivatives of energy by Euler angles to the resulting torque on a body?

-
Can you expand on why this doesn't work? It seems like it should, if you're just taking the derivative for an object at a particular orientation. (By the way you can just use regular derivatives in the question, we don't need to deal with numeric derivatives.) – David Z May 27 '12 at 5:35
I need this for testing. See my answer below. – abc360 May 27 '12 at 15:19

## 2 Answers

OK. I found the answer:

$$\partial V/\partial \theta = N_x \cos \psi - N_y \sin \psi$$ $$\partial V/\partial \phi = N_x \sin \theta \sin \psi + N_y \sin \theta \cos \psi + N_z \cos \theta$$ $$\partial V/\partial \psi = N_z$$

Where $\theta, \psi, \phi$ are Euler angles and $N_x, N_y, N_z$ are Torque components.

-
You should specify that you are talking about torques for each body with the origin at the center of mass. Also, the formalism you are using (Euler angles) is suboptimal, because it isn't very simple to rotate arbitrarily--- the best way to describe the orientation of the body is by an SO(3) (rotation) matrix R that tells you how the body axes is rotated relative to the coordinate axes. The Euler angles are just one set of coordinates for SO(3), and they aren't the best. The matrix itself is better, especially in a numerical implementation. – Ron Maimon May 27 '12 at 22:45
Yes, the torques are with the origin at the COM. I know about rotation matrices and unit quaternions to describe rotations. In fact Euler angles are never used in numerical integration of dynamics because there are terms where the division by increment of Euler angle is performed which causes numerical noise. – abc360 May 27 '12 at 23:46
BTW- The trig functions in are written \sin, \cos and so on in LaTeX. I'll fix it for you on this one. – dmckee May 27 '12 at 23:57

The mathematically correct method is to use a Jacobian, which is a matrix of partial derivatives. It works by the chain rule, which for a gradient vector looks like: $$\begin{pmatrix} \frac{\partial V(\theta,\psi,\phi)}{\partial x} \\ \frac{\partial V(\theta,\psi,\phi)}{\partial y} \\ \frac{\partial V(\theta,\psi,\phi)}{\partial z} \end{pmatrix} = \begin{pmatrix} \frac{\partial \theta}{\partial x} & \frac{\partial \psi}{\partial x} & \frac{\partial \phi}{\partial x} \\ \frac{\partial \theta}{\partial y} & \frac{\partial \psi}{\partial y} & \frac{\partial \phi}{\partial y} \\ \frac{\partial \theta}{\partial z} & \frac{\partial \psi}{\partial z} & \frac{\partial \phi}{\partial z} \end{pmatrix}\begin{pmatrix} \frac{\partial V(\theta,\psi,\phi)}{\partial \theta} \\ \frac{\partial V(\theta,\psi,\phi)}{\partial \psi} \\ \frac{\partial V(\theta,\psi,\phi)}{\partial \phi} \end{pmatrix}$$ The trouble is expressing the derivatives of $\theta$, $\psi$, and $\phi$ in terms of x, y, and z. Rotations about the fixed x, y, and z axes cannot form the Euler angles, so how do you find the derivatives? Seriously, I'm asking, I've been trying to figure it out for two weeks now.

-