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.special.elliptic.jacobi.CopolarN;
21 import org.hipparchus.special.elliptic.jacobi.JacobiElliptic;
22 import org.hipparchus.special.elliptic.jacobi.JacobiEllipticBuilder;
23 import org.hipparchus.util.FastMath;
24
25 /**
26 * This class is used in the junit tests for the ODE integrators.
27
28 * <p>This specific problem correspond to torque-free motion of a solid body
29 * with moments of inertia I₁, I₂, and I₃ with respect to body axes x, y, and z.
30 * We use here the notations from Landau and Lifchitz Course of Theoretical Physics,
31 * Mechanics vol 1.
32 * </p>
33 * <p>
34 * The equations of torque-free motion are given in the solid body frame by
35 * equation 36.5:
36 * <pre>
37 * I₁ dΩ₁/dt + (I₃ - I₂) Ω₂ Ω₃ = 0
38 * I₂ dΩ₂/dt + (I₁ - I₃) Ω₃ Ω₁ = 0
39 * I₃ dΩ₃/dt + (I₂ - I₁) Ω₁ Ω₂ = 0
40 * </pre>
41 * <p>
42 * The moments of inertia and initial conditions are: I₁ = 3/8, I₂ = 1/2, I₃ = 5/8,
43 * Ω₁ = 5, Ω₂ = 0, Ω₃ = 4. This corresponds to a motion with angular velocity
44 * describing a large polhode around Z axis in solid body frame. The motion is almost
45 * unstable as M² is only slightly greater than 2EI₂. Increasing Ω₁ to √(80/3) ≈ 5.16398
46 * would imply M²=2EI₂, and the polhode would degenerate to two intersecting ellipses.
47 * </p>
48 * <p>
49 * The torque-free motion can be solved analytically using Jacobi elliptic functions
50 * (in the rotating body frame)
51 * <pre>
52 * τ = t √([I₃-I₂][M²-2EI₁]/[I₁I₂I₃])
53 * Ω₁ (τ) = √([2EI₃-M²]/[I₁(I₃-I₁)]) cn(τ)
54 * Ω₂ (τ) = √([2EI₃-M²]/[I₂(I₃-I₂)]) sn(τ)
55 * Ω₃ (τ) = √([M²-2EI₁]/[I₃(I₃-I₁)]) dn(τ)
56 * </pre>
57 * </p>
58 * <p>
59 * This problem solves only solves the rotation rate part, whereas {@code TestProblem8}
60 * solves the full motion (rotation rate and rotation).
61 * </p>
62
63 */
64 public class TestProblem7 extends TestProblemAbstract {
65
66 /** Moments of inertia. */
67 final double i1;
68 final double i2;
69 final double i3;
70
71 /** Twice the angular kinetic energy. */
72 final double twoE;
73
74 final double m2;
75
76 /** Time scaling factor. */
77 final double tScale;
78
79 /** State scaling factors. */
80 final double o1Scale;
81 final double o2Scale;
82 final double o3Scale;
83
84 /** Jacobi elliptic function. */
85 final JacobiElliptic jacobi;
86
87 /**
88 * Simple constructor.
89 */
90 public TestProblem7() {
91
92 super(0.0, new double[] { 5.0, 0.0, 4.0 }, 4.0, new double[] { 1.0, 1.0, 1.0 });
93 i1 = 3.0 / 8.0;
94 i2 = 1.0 / 2.0;
95 i3 = 5.0 / 8.0;
96
97 final double[] s0 = getInitialState().getPrimaryState();
98 final double o12 = s0[0] * s0[0];
99 final double o22 = s0[1] * s0[1];
100 final double o32 = s0[2] * s0[2];
101 twoE = i1 * o12 + i2 * o22 + i3 * o32;
102 m2 = i1 * i1 * o12 + i2 * i2 * o22 + i3 * i3 * o32;
103 tScale = FastMath.sqrt((i3 - i2) * (m2 - twoE * i1) / (i1 * i2 * i3));
104 o1Scale = FastMath.sqrt((twoE * i3 - m2) / (i1 * (i3 - i1)));
105 o2Scale = FastMath.sqrt((twoE * i3 - m2) / (i2 * (i3 - i2)));
106 o3Scale = FastMath.sqrt((m2 - twoE * i1) / (i3 * (i3 - i1)));
107
108 final double k2 = (i2 - i1) * (twoE * i3 - m2) / ((i3 - i2) * (m2 - twoE * i1));
109 jacobi = JacobiEllipticBuilder.build(k2);
110
111 }
112
113 @Override
114 public double[] doComputeDerivatives(double t, double[] y) {
115
116 final double[] yDot = new double[getDimension()];
117
118 // compute the derivatives using Euler equations
119 yDot[0] = y[1] * y[2] * (i2 - i3) / i1;
120 yDot[1] = y[2] * y[0] * (i3 - i1) / i2;
121 yDot[2] = y[0] * y[1] * (i1 - i2) / i3;
122
123 return yDot;
124
125 }
126
127 @Override
128 public double[] computeTheoreticalState(double t) {
129 final CopolarN valuesN = jacobi.valuesN(t * tScale);
130 return new double[] {
131 o1Scale * valuesN.cn(),
132 o2Scale * valuesN.sn(),
133 o3Scale * valuesN.dn()
134 };
135 }
136
137 }