View Javadoc
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    *      https://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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.interpolation;
23  
24  import java.io.Serializable;
25  
26  import org.hipparchus.analysis.polynomials.PolynomialFunctionLagrangeForm;
27  import org.hipparchus.analysis.polynomials.PolynomialFunctionNewtonForm;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  
30  /**
31   * Implements the <a href=
32   * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
33   * Divided Difference Algorithm</a> for interpolation of real univariate
34   * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
35   * ISBN 038795452X, chapter 2.
36   * <p>
37   * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
38   * this class provides an easy-to-use interface to it.</p>
39   *
40   */
41  public class DividedDifferenceInterpolator
42      implements UnivariateInterpolator, Serializable {
43      /** serializable version identifier */
44      private static final long serialVersionUID = 107049519551235069L;
45  
46      /** Empty constructor.
47       * <p>
48       * This constructor is not strictly necessary, but it prevents spurious
49       * javadoc warnings with JDK 18 and later.
50       * </p>
51       * @since 3.0
52       */
53      public DividedDifferenceInterpolator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
54          // nothing to do
55      }
56  
57      /**
58       * Compute an interpolating function for the dataset.
59       *
60       * @param x Interpolating points array.
61       * @param y Interpolating values array.
62       * @return a function which interpolates the dataset.
63       * @throws MathIllegalArgumentException if the array lengths are different.
64       * @throws MathIllegalArgumentException if the number of points is less than 2.
65       * @throws MathIllegalArgumentException if {@code x} is not sorted in
66       * strictly increasing order.
67       */
68      @Override
69      public PolynomialFunctionNewtonForm interpolate(double x[], double y[])
70          throws MathIllegalArgumentException {
71          /**
72           * a[] and c[] are defined in the general formula of Newton form:
73           * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
74           *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
75           */
76          PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
77  
78          /**
79           * When used for interpolation, the Newton form formula becomes
80           * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
81           *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
82           * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
83           * <p>
84           * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
85           */
86          final double[] c = new double[x.length-1];
87          System.arraycopy(x, 0, c, 0, c.length);
88  
89          final double[] a = computeDividedDifference(x, y);
90          return new PolynomialFunctionNewtonForm(a, c);
91      }
92  
93      /**
94       * Return a copy of the divided difference array.
95       * <p>
96       * The divided difference array is defined recursively by <pre>
97       * f[x0] = f(x0)
98       * f[x0,x1,...,xk] = (f[x1,...,xk] - f[x0,...,x[k-1]]) / (xk - x0)
99       * </pre>
100      * <p>
101      * The computational complexity is \(O(n^2)\) where \(n\) is the common
102      * length of {@code x} and {@code y}.</p>
103      *
104      * @param x Interpolating points array.
105      * @param y Interpolating values array.
106      * @return a fresh copy of the divided difference array.
107      * @throws MathIllegalArgumentException if the array lengths are different.
108      * @throws MathIllegalArgumentException if the number of points is less than 2.
109      * @throws MathIllegalArgumentException
110      * if {@code x} is not sorted in strictly increasing order.
111      */
112     protected static double[] computeDividedDifference(final double x[], final double y[])
113         throws MathIllegalArgumentException {
114         PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
115 
116         final double[] divdiff = y.clone(); // initialization
117 
118         final int n = x.length;
119         final double[] a = new double [n];
120         a[0] = divdiff[0];
121         for (int i = 1; i < n; i++) {
122             for (int j = 0; j < n-i; j++) {
123                 final double denominator = x[j+i] - x[j];
124                 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
125             }
126             a[i] = divdiff[0];
127         }
128 
129         return a;
130     }
131 }