View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.hipparchus.ode;
19  
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.Field;
22  import org.hipparchus.special.elliptic.jacobi.FieldCopolarN;
23  import org.hipparchus.special.elliptic.jacobi.FieldJacobiElliptic;
24  import org.hipparchus.special.elliptic.jacobi.JacobiEllipticBuilder;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.MathArrays;
27  
28  /**
29   * This class is used in the junit tests for the ODE integrators.
30  
31   * <p>This specific problem correspond to torque-free motion of a solid body
32   * with moments of inertia I₁, I₂, and I₃ with respect to body axes x, y, and z.
33   * We use here the notations from Landau and Lifchitz Course of Theoretical Physics,
34   * Mechanics vol 1.
35   * </p>
36   * <p>
37   * The equations of torque-free motion are given in the solid body frame by
38   * equation 36.5:
39   * <pre>
40   *    I₁ dΩ₁/dt + (I₃ - I₂) Ω₂ Ω₃ = 0
41   *    I₂ dΩ₂/dt + (I₁ - I₃) Ω₃ Ω₁ = 0
42   *    I₃ dΩ₃/dt + (I₂ - I₁) Ω₁ Ω₂ = 0
43   * </pre>
44   * <p>
45   * The moments of inertia and initial conditions are: I₁ = 3/8, I₂ = 1/2, I₃ = 5/8,
46   * Ω₁ = 5, Ω₂ = 0, Ω₃ = 4. This corresponds to a motion with angular velocity
47   * describing a large polhode around Z axis in solid body frame. The motion is almost
48   * unstable as M² is only slightly greater than 2EI₂. Increasing Ω₁ to √(80/3) ≈ 5.16398
49   * would imply M²=2EI₂, and the polhode would degenerate to two intersecting ellipses.
50   * </p>
51   * <p>
52   * The torque-free motion can be solved analytically using Jacobi elliptic functions
53   * (in the rotating body frame)
54   * <pre>
55   *   τ      = t √([I₃-I₂][M²-2EI₁]/[I₁I₂I₃])
56   *   Ω₁ (τ) = √([2EI₃-M²]/[I₁(I₃-I₁)]) cn(τ)
57   *   Ω₂ (τ) = √([2EI₃-M²]/[I₂(I₃-I₂)]) sn(τ)
58   *   Ω₃ (τ) = √([M²-2EI₁]/[I₃(I₃-I₁)]) dn(τ)
59   * </pre>
60   * </p>
61   * <p>
62   * This problem solves only solves the rotation rate part, whereas {@code FieldTestProblem8}
63   * solves the full motion (rotation rate and rotation).
64   * </p>
65   * @param <T> the type of the field elements
66   */
67  public class TestFieldProblem7<T extends CalculusFieldElement<T>>
68      extends TestFieldProblemAbstract<T> {
69  
70      /** Moments of inertia. */
71      final T i1;
72      final T i2;
73      final T i3;
74  
75      /** Twice the angular kinetic energy. */
76      final T twoE;
77  
78      final T m2;
79  
80      /** Time scaling factor. */
81      final T tScale;
82  
83      /** State scaling factors. */
84      final T o1Scale;
85      final T o2Scale;
86      final T o3Scale;
87  
88      /** Jacobi elliptic function. */
89      final FieldJacobiElliptic<T> jacobi;
90  
91      /**
92       * Simple constructor.
93       * @param field field to which elements belong
94       */
95      public TestFieldProblem7(Field<T> field) {
96          super(convert(field, 0.0),
97                convert(field, new double[] { 5.0, 0.0, 4.0 }),
98                convert(field, 4.0),
99                convert(field, new double[] { 1.0, 1.0, 1.0 }));
100         i1 = convert(field, 3.0 / 8.0);
101         i2 = convert(field, 1.0 / 2.0);
102         i3 = convert(field, 5.0 / 8.0);
103 
104         final T[] s0 = getInitialState().getPrimaryState();
105         final T o12 = s0[0].multiply(s0[0]);
106         final T o22 = s0[1].multiply(s0[1]);
107         final T o32 = s0[2].multiply(s0[2]);
108         twoE    =  i1.multiply(o12).add(i2.multiply(o22)).add(i3.multiply(o32));
109         m2      =  i1.multiply(i1.multiply(o12)).add(i2.multiply(i2.multiply(o22))).add(i3.multiply(i3.multiply(o32)));
110         tScale  = FastMath.sqrt(i3.subtract(i2).multiply(m2.subtract(twoE.multiply(i1))).divide(i1.multiply(i2).multiply(i3)));
111         o1Scale = FastMath.sqrt(twoE.multiply(i3).subtract(m2).divide(i1.multiply(i3.subtract(i1))));
112         o2Scale = FastMath.sqrt(twoE.multiply(i3).subtract(m2).divide(i2.multiply(i3.subtract(i2))));
113         o3Scale = FastMath.sqrt(m2.subtract(twoE.multiply(i1)).divide(i3.multiply(i3.subtract(i1))));
114 
115         final T k2 = i2.subtract(i1).multiply(twoE.multiply(i3).subtract(m2)).divide(i3.subtract(i2).multiply(m2.subtract(twoE.multiply(i1))));
116         jacobi = JacobiEllipticBuilder.build(k2);
117     }
118 
119     @Override
120     public T[] doComputeDerivatives(T t, T[] y) {
121 
122         final T[] yDot = MathArrays.buildArray(getField(), getDimension());
123 
124         // compute the derivatives using Euler equations
125         yDot[0] = y[1].multiply(y[2]).multiply(i2.subtract(i3)).divide(i1);
126         yDot[1] = y[2].multiply(y[0]).multiply(i3.subtract(i1)).divide(i2);
127         yDot[2] = y[0].multiply(y[1]).multiply(i1.subtract(i2)).divide(i3);
128 
129         return yDot;
130 
131     }
132 
133     @Override
134     public T[] computeTheoreticalState(T t) {
135 
136         final T[] y = MathArrays.buildArray(getField(), getDimension());
137 
138         final FieldCopolarN<T> valuesN = jacobi.valuesN(t.multiply(tScale));
139         y[0] = o1Scale.multiply(valuesN.cn());
140         y[1] = o2Scale.multiply(valuesN.sn());
141         y[2] = o3Scale.multiply(valuesN.dn());
142 
143         return y;
144 
145     }
146 
147 }