001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.analysis.function;
019
020 import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
021 import org.apache.commons.math3.analysis.FunctionUtils;
022 import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
023 import org.apache.commons.math3.analysis.UnivariateFunction;
024 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
025 import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
026 import org.apache.commons.math3.exception.DimensionMismatchException;
027 import org.apache.commons.math3.exception.NullArgumentException;
028 import org.apache.commons.math3.util.FastMath;
029
030 /**
031 * <a href="http://en.wikipedia.org/wiki/Harmonic_oscillator">
032 * simple harmonic oscillator</a> function.
033 *
034 * @since 3.0
035 * @version $Id: HarmonicOscillator.java 1383441 2012-09-11 14:56:39Z luc $
036 */
037 public class HarmonicOscillator implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction {
038 /** Amplitude. */
039 private final double amplitude;
040 /** Angular frequency. */
041 private final double omega;
042 /** Phase. */
043 private final double phase;
044
045 /**
046 * Harmonic oscillator function.
047 *
048 * @param amplitude Amplitude.
049 * @param omega Angular frequency.
050 * @param phase Phase.
051 */
052 public HarmonicOscillator(double amplitude,
053 double omega,
054 double phase) {
055 this.amplitude = amplitude;
056 this.omega = omega;
057 this.phase = phase;
058 }
059
060 /** {@inheritDoc} */
061 public double value(double x) {
062 return value(omega * x + phase, amplitude);
063 }
064
065 /** {@inheritDoc}
066 * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)}
067 */
068 @Deprecated
069 public UnivariateFunction derivative() {
070 return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative();
071 }
072
073 /**
074 * Parametric function where the input array contains the parameters of
075 * the harmonic oscillator function, ordered as follows:
076 * <ul>
077 * <li>Amplitude</li>
078 * <li>Angular frequency</li>
079 * <li>Phase</li>
080 * </ul>
081 */
082 public static class Parametric implements ParametricUnivariateFunction {
083 /**
084 * Computes the value of the harmonic oscillator at {@code x}.
085 *
086 * @param x Value for which the function must be computed.
087 * @param param Values of norm, mean and standard deviation.
088 * @return the value of the function.
089 * @throws NullArgumentException if {@code param} is {@code null}.
090 * @throws DimensionMismatchException if the size of {@code param} is
091 * not 3.
092 */
093 public double value(double x, double ... param)
094 throws NullArgumentException,
095 DimensionMismatchException {
096 validateParameters(param);
097 return HarmonicOscillator.value(x * param[1] + param[2], param[0]);
098 }
099
100 /**
101 * Computes the value of the gradient at {@code x}.
102 * The components of the gradient vector are the partial
103 * derivatives of the function with respect to each of the
104 * <em>parameters</em> (amplitude, angular frequency and phase).
105 *
106 * @param x Value at which the gradient must be computed.
107 * @param param Values of amplitude, angular frequency and phase.
108 * @return the gradient vector at {@code x}.
109 * @throws NullArgumentException if {@code param} is {@code null}.
110 * @throws DimensionMismatchException if the size of {@code param} is
111 * not 3.
112 */
113 public double[] gradient(double x, double ... param)
114 throws NullArgumentException,
115 DimensionMismatchException {
116 validateParameters(param);
117
118 final double amplitude = param[0];
119 final double omega = param[1];
120 final double phase = param[2];
121
122 final double xTimesOmegaPlusPhase = omega * x + phase;
123 final double a = HarmonicOscillator.value(xTimesOmegaPlusPhase, 1);
124 final double p = -amplitude * FastMath.sin(xTimesOmegaPlusPhase);
125 final double w = p * x;
126
127 return new double[] { a, w, p };
128 }
129
130 /**
131 * Validates parameters to ensure they are appropriate for the evaluation of
132 * the {@link #value(double,double[])} and {@link #gradient(double,double[])}
133 * methods.
134 *
135 * @param param Values of norm, mean and standard deviation.
136 * @throws NullArgumentException if {@code param} is {@code null}.
137 * @throws DimensionMismatchException if the size of {@code param} is
138 * not 3.
139 */
140 private void validateParameters(double[] param)
141 throws NullArgumentException,
142 DimensionMismatchException {
143 if (param == null) {
144 throw new NullArgumentException();
145 }
146 if (param.length != 3) {
147 throw new DimensionMismatchException(param.length, 3);
148 }
149 }
150 }
151
152 /**
153 * @param xTimesOmegaPlusPhase {@code omega * x + phase}.
154 * @param amplitude Amplitude.
155 * @return the value of the harmonic oscillator function at {@code x}.
156 */
157 private static double value(double xTimesOmegaPlusPhase,
158 double amplitude) {
159 return amplitude * FastMath.cos(xTimesOmegaPlusPhase);
160 }
161
162 /** {@inheritDoc}
163 * @since 3.1
164 */
165 public DerivativeStructure value(final DerivativeStructure t) {
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 }