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