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 }