WOLFRAM LANGUAGE TUTORIAL

Rigid Body Solvers

Introduction

The equations of motion for a free rigid body whose center of mass is at the origin are given by the following Euler equations (see [MR99]).

Two quadratic first integrals of the system are:

The first constraint effectively confines the motion from to a sphere. The second constraint represents the kinetic energy of the system and, in conjunction with the first invariant, effectively confines the motion to ellipsoids on the sphere.

Numerical experiments for various methods are given in [HLW02] and a variety of NDSolve methods will now be compared.

Manifold Generation and Utility Functions

Load some useful packages.
In[6]:=
Click for copyable input
Define Euler's equations for rigid body motion together with the invariants of the system.
In[8]:=
Click for copyable input
The equations of motion evolve as closed curves on the unit sphere. This generates a three-dimensional graphics object to represent the unit sphere.
In[13]:=
Click for copyable input
This function superimposes a solution from NDSolve on a given manifold.
In[14]:=
Click for copyable input
This function plots the various solution components.
In[15]:=
Click for copyable input

Method Comparison

Various integration methods can be used to solve Euler's equations and they each have different associated costs and different dynamical properties.

Adams Multistep Method

Here an Adams method is used to solve the equations of motion.
In[21]:=
Click for copyable input
This shows the solution trajectory by superimposing it on the unit sphere.
In[22]:=
Click for copyable input
Out[22]=
The solution appears visually to give a closed curve on the sphere. However, a plot of the error reveals that neither constraint is conserved particularly well.
In[23]:=
Click for copyable input
Out[23]=

Euler and Implicit Midpoint Methods

This solves the equations of motion using Euler's method with a specified fixed step size.
In[16]:=
Click for copyable input
This solves the equations of motion using the implicit midpoint method with a specified fixed step size.
In[17]:=
Click for copyable input
This shows the superimposition on the unit sphere of the numerical solution of the equations of motion for Euler's method (left) and the implicit midpoint rule (right).
In[19]:=
Click for copyable input
Out[21]=
This shows the components of the numerical solution using Euler's method (left) and the implicit midpoint rule (right).
In[30]:=
Click for copyable input
Out[32]=

Orthogonal Projection Method

Here the method is used to solve the equations.
In[33]:=
Click for copyable input
Only the orthogonal constraint is conserved so the curve is not closed.
In[34]:=
Click for copyable input
Out[34]=
Plotting the error in the invariants against time, it can be seen that the orthogonal projection method conserves only one of the two invariants.
In[35]:=
Click for copyable input
Out[35]=

Projection Method

The method takes a set of constraints and projects the solution onto a manifold at the end of each integration step.

Generally all the invariants of the problem should be used in the projection; otherwise the numerical solution may actually be qualitatively worse than the unprojected solution.

The following specifies the integration method and defers determination of the constraints until the invocation of NDSolve.
In[36]:=
Click for copyable input

Projecting One Constraint

This projects the first constraint onto the manifold.
In[37]:=
Click for copyable input
Out[39]=
Only the first invariant is conserved.
In[40]:=
Click for copyable input
Out[40]=
This projects the second constraint onto the manifold.
In[41]:=
Click for copyable input
Out[43]=
Only the second invariant is conserved.
In[44]:=
Click for copyable input
Out[44]=

Projecting Multiple Constraints

This projects both constraints onto the manifold.
In[45]:=
Click for copyable input
Out[47]=
Now both invariants are conserved.
In[48]:=
Click for copyable input
Out[48]=

"Splitting" Method

A splitting that yields an efficient explicit integration method was derived independently by McLachlan [M93] and Reich [R93].

Write the flow of an ODE as .

The differential system is split into three components, , , and , each of which is Hamiltonian and can be solved exactly.

The Hamiltonian systems are solved and recombined at each integration step as:

This defines an appropriate splitting into Hamiltonian vector fields.
In[49]:=
Click for copyable input
Here is the differential system for Euler's equations.
In[58]:=
Click for copyable input
Out[58]=
Here are the three split vector fields.
In[59]:=
Click for copyable input
Out[59]=
In[60]:=
Click for copyable input
Out[60]=
In[61]:=
Click for copyable input
Out[61]=

Solution

This defines a symmetric second-order splitting method. The coefficients are automatically determined from the structure of the equations and are an extension of the Strang splitting.
In[62]:=
Click for copyable input
This solves the system and graphically displays the solution.
In[63]:=
Click for copyable input
Out[64]=
One of the invariants is preserved up to roundoff while the error in the second invariant remains bounded.
In[65]:=
Click for copyable input
Out[65]=