    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    }