Numerical solution of Euler’s rotation equations for a rigid body about a ﬁxed point

Finding a solution for Euler’s equations is a classic mechanics problem. This study revisits the problem with numerical approaches. For ease of teaching and research, a Maple code comprising 2 lines is written to ﬁnd a numerical solution for the problem. The study’s results are validated by comparing these with previous studies. Our results conﬁrm the correctness of the principle of maximum moment of inertia of the rotating body, which is veriﬁed by thermodynamics. As an essential part of this study, the Maple code is provided.


I. INTRODUCTION
The celebrated Euler equations of motion for a rigid body consists of three nonlinear equations, coupled with differential equations, which are known as one of the famous problem in classical mechanics 1 . Solving the Euler equations has attracted the attention of great scientists for the last 300+ years. Despite great efforts in this respect, a complete general analytical solution has yet to be found [1][2][3][4][5][6][7][8][9][10][11] .
Special cases for which solutions have been found [1][2][3][4][5][6][7][9][10][11] , which include the torque-free rigid body, and the three (or four) famous integrable cases solved by Euler, Lagrange, and Kovalevskaya 3 . The Lagrange and Kovalevskaya case that symmetric rigid rotor's two principal moments of inertia are equal to each other, and double that of the third one, namely I 1 = I 2 = 2I 3 , and the Euler case that all the applied torques are zero (torque-free precession of the rotation axis of the rigid rotor) 1,2,4 . The desire of searching exact solution on the problem has never been ending, recently, Ershkov 14 , who proposed a new and exact solution for the Euler equations.
Besides the above few integrable cases, for the most disintegrable case or the heavy asymmetric case, some approximate analytical solutions have been proposed by, for example, Amer and Abady 15,16 studied analytical solutions for the rigid body motion in the presence of gyrostatic torques in terms of the axes of rotation. Tsiotras and Longuski 12,13 , who considered the time evolution of the angular velocity of a spinning rigid body, subject to the torques of the three axes, and proposed an asymptotic analytic solution.
After more than 300 years of investigations and studying the Euler's equation, it is concluded that no closed solution has been obtained for the general case, while numerical solutions have been used. Although the closed and approximate analytic solutions have great academic value, they have certain limitations, since all of them involve complicated computations and variable transformations. They may be mathematically correct but the difficulty lies in the computational perspective. In order to overcome this situation, it would be natural thinking to use numerical methods for the general case of the Euler equation. If we use numerical methods, the computer code must be written. However, if one studies the literature carefully, it would not be difficult to find that there is no simple computer program, comprising a few line, which has been reported. Therefore, it is highly demanded to have a simple straight-forward numerical program that can be used to compute the problem simply by a click.
This study provides a general Maple code for the numerical solution of the Euler equation, namely a simple Maple code comprising merely 2 lines. In Section 3 we compare our results with that of the benchmark from Tsiotras and Longuski's studies 13 . In section 4 we investigate the torques as functions of time. In Section 5 we compute an asymmetrical top under the action of a gyrostatic moment, and finally, with discussions and a conclusion. For ease of teaching and research, all Maple codes are provided.

II. EULER'S ROTATION EQUATION AND MAPLE CODE
In mechanics, a rigid body may be defined as a system of particles, as the distances between the particles do not vary. To describe the motion of a rigid body, we use two systems of co-ordinates: a "fixed" (i.e. inertial) system, XY Z, and a moving system, x 1 = x, x 2 = y, and x 3 = z which is supposed to be rigidly fixed in the body, and to participate in its motion. The origin of the moving system may conveniently be taken to coincide with the centre of body's mass 7 .
A non-symmetric spinning body in space, subject to constant torques and non-zero initial conditions, using Euler's equations of motion for a rotating rigid body with principal axes at the center of mass, produce the following The principal torques acting on the rigid body, namely K 1 , K 2 and K 3 . The principal angular velocity components in the same frame, namely Ω 1 , Ω 2 and Ω 3 . These quantities are defined in the moving system x 1 = x, x 2 = y, and x 3 = z, namely in the body-fixed frame 1-10 . Euler's equations in 1 are nonlinear differential equations, whose general analytical solutions have not been obtained and might well be impossible to obtain 1-10 . Hence, for general cases, or for the asymmetrical top of the Euler equations, numerical solutions are perhaps the only choice.
To obtain the numerical solutions, one can write codes by using various software and languages such as Fortran, C++, Matlab, MATHEMATICA and so on. To obtain simple and short codes, the symbolic software, Maple, was used. Hence, a Maple code was produced for the Euler equations. The code is shown in Table I below.
The code uses the solver x-rkf45 which is implemented based on the fourth-order Runge-Kutta method. The above code was used to compute all case studies and they are shown in the forthcoming sections for given moment of inertia I k , torques K k and initial condition x0, y0 and z0.

III. TSIOTRAS AND LONGUSKI'S BENCHMARK 13
In 1996 Tsiotras and Longuski 13 proposed a novel approximate solution for the asymmetrical top of Euler's equations, and presented a numerical example. To validate our Maple code, we recomputed the same example and compared it with those of Tsiotras and Longuski 13 . The data that we used is shown in Table II below.
Using our code for the problem, the results for the different ranges of time are shown in the Fig. 1 below.
The kinetic energy T of the asymmetrical top is shown in Table 2 below.
Comparisons of the Ω 1 profiles with those of Tsiotras and Longuski 13 are in shown in Fig.3 below. Our solution concurs with that of Tsiotras and Longuski 13 .   Besides the numerical data is shown in Table IV below, particularly for purpose of future comparative studies.

IV. ELLIPSOID TOP ROTATION
This section demonstrates rotation of ellipsoid top, whose typical shape is shown in Fig. 5 below.
In reality, the torque might depend on the rotation velocity, for instance, in viscous environments. However, there are no results were reported for the torques as functions of time. This section demonstrates three examples in this regards, as shown below.

A. Case with one variable torque
If the torque is a function of rotation velocitiy, the code can provide solutions without any difficulty. In this case study, we assume that K 1 = −0.1Ω 1 (t) + 0.05. The other data is shown in Table IV below.
The results are shown in Fig. 6 below.

B. Case with two variable torques
Here both K 1 and K 2 are the function of time. All this is presented in Table V below and results are shown in Fig. 7 below.

C. Case with three variable torques
If all three torque are a function of rotation velocitiy, the code can also give solutions without any difficulty. All data are in Table VI.  Table V   TABLE VI. Data for case with three variable torque Semiaxes a = 1, b = 2, c = 3 Total mass µ Moment of initial The results are shown in Fig 8 below.

V. AN ASYMMETRIC TOP UNDER THE ACTION OF A GYROSTATIC MOMENT
If an symmetrical top acted in the action of a gyrostatic moment, all three torques would be function of rotation FIG. 8. Ω profiles of the case in Table VI   TABLE VII. Data for an asymmetric top under the action of a gyrostatic moment FIG. 9. Ω profiles of the case in Table VII velocitiy, the code could provide solutions without any difficulty. The data is presented in Table VII below. The results are shown in Fig.9 below. If we only change one semiaxes, namely c = k, k = 1, , 2, 3, 4, 5 and keep all the others unchanged, then the results will be as follows, shown in Figures below, namely 10, 11 and 12.

VI. THERMODYNAMICAL CONSIDERATION OF ROTATING BODY
The above numerical studies show that longer axes of the rigid body are sensitive to any change of the torque and that little variation leads to a dramatic change in the rotation velocity. In other words, the maximum moment of inertia is the key factor in controlling rigid body rotation. Rotation about the maximum principal moment of inertia represents the minimum possible kinetic energy for a given angular momentum that a system can possibly have. This principle has been confirmed from thermodynamical perspectives. Since the entropy S of a body is a function of its internal energy E in = E − M 2 /(2I), namely S = S(E in ) = S(E − M 2 /(2I)), where E is total energy, M is moment and I is moment of inertia. Because the body is a closed system, thus its total energy and angular momentum are conserved, and hence the entropy must have the maximum value possible for the given M and E. Therefore, the equilibrium rotation of the body takes place about the axis with respect to which the moment of inertia has the greatest possible value 17 . This is called the principle of maximum moment of inertia of the rotating body.

VII. CONCLUSIONS
In conclusion, the moment of inertia is the key factor involved in controlling rigid body rotation 17 . Proper design of the body shape can achieve the desired performance of a top. The comparison confirmed that Tsiotras and Longuski's results 13 were highly accurate. For ease of teaching and research, we have written a Maple code to find a solution for the Euler equations. The provided Maple code can be used generally to solve any problem related to Euler's equations, while, importantly, it is also simple and user-friendly.
Availability of data: The data that supports the findings of this study is available from the corresponding author, upon reasonable request.
Conflict of interest: The author declares that he/she has no known competing financial interests or personal relationships that may have influenced the work reported in this paper.
B.-H. Sun: Conceptualization, Methodology, Formulations, Formal analysis, Funding acquisition, Investigation, Writing Original draft preparation, Writing, Reviewing and Editing, and all relevant associated works.