1 /*
2 * Licensed to the Apache Software Foundation (ASF) 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 ASF 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 * http://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.apache.commons.math4.legacy.analysis.function;
19
20 import org.apache.commons.math4.legacy.analysis.ParametricUnivariateFunction;
21 import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure;
22 import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateDifferentiableFunction;
23 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
24 import org.apache.commons.math4.legacy.exception.NullArgumentException;
25 import org.apache.commons.math4.core.jdkmath.JdkMath;
26
27 /**
28 * <a href="http://en.wikipedia.org/wiki/Harmonic_oscillator">
29 * simple harmonic oscillator</a> function.
30 *
31 * @since 3.0
32 */
33 public class HarmonicOscillator implements UnivariateDifferentiableFunction {
34 /** Amplitude. */
35 private final double amplitude;
36 /** Angular frequency. */
37 private final double omega;
38 /** Phase. */
39 private final double phase;
40
41 /**
42 * Harmonic oscillator function.
43 *
44 * @param amplitude Amplitude.
45 * @param omega Angular frequency.
46 * @param phase Phase.
47 */
48 public HarmonicOscillator(double amplitude,
49 double omega,
50 double phase) {
51 this.amplitude = amplitude;
52 this.omega = omega;
53 this.phase = phase;
54 }
55
56 /** {@inheritDoc} */
57 @Override
58 public double value(double x) {
59 return value(omega * x + phase, amplitude);
60 }
61
62 /**
63 * Parametric function where the input array contains the parameters of
64 * the harmonic oscillator function. Ordered as follows:
65 * <ul>
66 * <li>Amplitude</li>
67 * <li>Angular frequency</li>
68 * <li>Phase</li>
69 * </ul>
70 */
71 public static class Parametric implements ParametricUnivariateFunction {
72 /**
73 * Computes the value of the harmonic oscillator at {@code x}.
74 *
75 * @param x Value for which the function must be computed.
76 * @param param Values of norm, mean and standard deviation.
77 * @return the value of the function.
78 * @throws NullArgumentException if {@code param} is {@code null}.
79 * @throws DimensionMismatchException if the size of {@code param} is
80 * not 3.
81 */
82 @Override
83 public double value(double x, double ... param)
84 throws NullArgumentException,
85 DimensionMismatchException {
86 validateParameters(param);
87 return HarmonicOscillator.value(x * param[1] + param[2], param[0]);
88 }
89
90 /**
91 * Computes the value of the gradient at {@code x}.
92 * The components of the gradient vector are the partial
93 * derivatives of the function with respect to each of the
94 * <em>parameters</em> (amplitude, angular frequency and phase).
95 *
96 * @param x Value at which the gradient must be computed.
97 * @param param Values of amplitude, angular frequency and phase.
98 * @return the gradient vector at {@code x}.
99 * @throws NullArgumentException if {@code param} is {@code null}.
100 * @throws DimensionMismatchException if the size of {@code param} is
101 * not 3.
102 */
103 @Override
104 public double[] gradient(double x, double ... param)
105 throws NullArgumentException,
106 DimensionMismatchException {
107 validateParameters(param);
108
109 final double amplitude = param[0];
110 final double omega = param[1];
111 final double phase = param[2];
112
113 final double xTimesOmegaPlusPhase = omega * x + phase;
114 final double a = HarmonicOscillator.value(xTimesOmegaPlusPhase, 1);
115 final double p = -amplitude * JdkMath.sin(xTimesOmegaPlusPhase);
116 final double w = p * x;
117
118 return new double[] { a, w, p };
119 }
120
121 /**
122 * Validates parameters to ensure they are appropriate for the evaluation of
123 * the {@link #value(double,double[])} and {@link #gradient(double,double[])}
124 * methods.
125 *
126 * @param param Values of norm, mean and standard deviation.
127 * @throws NullArgumentException if {@code param} is {@code null}.
128 * @throws DimensionMismatchException if the size of {@code param} is
129 * not 3.
130 */
131 private void validateParameters(double[] param)
132 throws NullArgumentException,
133 DimensionMismatchException {
134 if (param == null) {
135 throw new NullArgumentException();
136 }
137 if (param.length != 3) {
138 throw new DimensionMismatchException(param.length, 3);
139 }
140 }
141 }
142
143 /**
144 * @param xTimesOmegaPlusPhase {@code omega * x + phase}.
145 * @param amplitude Amplitude.
146 * @return the value of the harmonic oscillator function at {@code x}.
147 */
148 private static double value(double xTimesOmegaPlusPhase,
149 double amplitude) {
150 return amplitude * JdkMath.cos(xTimesOmegaPlusPhase);
151 }
152
153 /** {@inheritDoc}
154 * @since 3.1
155 */
156 @Override
157 public DerivativeStructure value(final DerivativeStructure t)
158 throws DimensionMismatchException {
159 final double x = t.getValue();
160 double[] f = new double[t.getOrder() + 1];
161
162 final double alpha = omega * x + phase;
163 f[0] = amplitude * JdkMath.cos(alpha);
164 if (f.length > 1) {
165 f[1] = -amplitude * omega * JdkMath.sin(alpha);
166 final double mo2 = - omega * omega;
167 for (int i = 2; i < f.length; ++i) {
168 f[i] = mo2 * f[i - 2];
169 }
170 }
171
172 return t.compose(f);
173 }
174 }