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 018package org.apache.commons.math3.analysis.function; 019 020import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; 021import org.apache.commons.math3.analysis.FunctionUtils; 022import org.apache.commons.math3.analysis.ParametricUnivariateFunction; 023import org.apache.commons.math3.analysis.UnivariateFunction; 024import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; 025import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; 026import org.apache.commons.math3.exception.DimensionMismatchException; 027import org.apache.commons.math3.exception.NullArgumentException; 028import 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 */ 036public class HarmonicOscillator implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction { 037 /** Amplitude. */ 038 private final double amplitude; 039 /** Angular frequency. */ 040 private final double omega; 041 /** Phase. */ 042 private final double phase; 043 044 /** 045 * Harmonic oscillator function. 046 * 047 * @param amplitude Amplitude. 048 * @param omega Angular frequency. 049 * @param phase Phase. 050 */ 051 public HarmonicOscillator(double amplitude, 052 double omega, 053 double phase) { 054 this.amplitude = amplitude; 055 this.omega = omega; 056 this.phase = phase; 057 } 058 059 /** {@inheritDoc} */ 060 public double value(double x) { 061 return value(omega * x + phase, amplitude); 062 } 063 064 /** {@inheritDoc} 065 * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)} 066 */ 067 @Deprecated 068 public UnivariateFunction derivative() { 069 return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative(); 070 } 071 072 /** 073 * Parametric function where the input array contains the parameters of 074 * the harmonic oscillator function, ordered as follows: 075 * <ul> 076 * <li>Amplitude</li> 077 * <li>Angular frequency</li> 078 * <li>Phase</li> 079 * </ul> 080 */ 081 public static class Parametric implements ParametricUnivariateFunction { 082 /** 083 * Computes the value of the harmonic oscillator at {@code x}. 084 * 085 * @param x Value for which the function must be computed. 086 * @param param Values of norm, mean and standard deviation. 087 * @return the value of the function. 088 * @throws NullArgumentException if {@code param} is {@code null}. 089 * @throws DimensionMismatchException if the size of {@code param} is 090 * not 3. 091 */ 092 public double value(double x, double ... param) 093 throws NullArgumentException, 094 DimensionMismatchException { 095 validateParameters(param); 096 return HarmonicOscillator.value(x * param[1] + param[2], param[0]); 097 } 098 099 /** 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 */ 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 */ 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 */ 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 */ 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}