001    package edu.nrao.sss.math;
002    
003    import java.awt.geom.Point2D;
004    import java.math.BigDecimal;
005    import java.math.MathContext;
006    import java.math.RoundingMode;
007    import java.util.ArrayList;
008    import java.util.List;
009    
010    /**
011     * An interpolated value calculator that uses Newton's divided-difference
012     * formula.
013     * <p>
014     * Under certain conditions this interpolator can have trouble with
015     * the first and last few intervals.  The {@link CubicSpline} may be
016     * a better interpolator for most situations.</p>
017     * <p>
018     * <b>Version Info:</b>
019     * <table style="margin-left:2em">
020     *   <tr><td>$Revision$</td></tr>
021     *   <tr><td>$Date$</td></tr>
022     *   <tr><td>$Author$</td></tr>
023     * </table></p>
024     * 
025     * @author David M. Harland
026     * @since 2007-04-27
027     */
028    public class NewtonDividedDifference
029      extends AbsInterp
030    {
031      private static final MathContext MC = new MathContext(100, RoundingMode.HALF_UP);
032      
033      private List<BigDecimal> coefficients;
034      private List<BigDecimal> xValues;
035    
036      /** Creates a new instance. */
037      public NewtonDividedDifference()
038      {
039        super();
040        
041        coefficients = new ArrayList<BigDecimal>();
042        xValues      = new ArrayList<BigDecimal>();
043      }
044    
045      /* (non-Javadoc)
046       * @see edu.nrao.sss.math.AbsInterp#getComputedValueFor(double)
047       */
048      double getComputedValueFor(double x)
049      {
050        BigDecimal X = new BigDecimal(Double.toString(x), MC);
051        
052        int n = this.size()-1;
053        
054        BigDecimal polynomial = BigDecimal.ZERO;
055        
056        for (int i=0; i <= n; i++)
057        {
058          BigDecimal term = coefficients.get(i);
059          
060          for (int j=0; j < i; j++)
061          {
062            BigDecimal diff = X.subtract(xValues.get(j));
063            term = term.multiply(diff);
064          }
065          
066          polynomial = polynomial.add(term);
067        }
068        
069        return polynomial.doubleValue();
070      }
071    
072      /**
073       * This algorithm was take from section 3.4 of <u>Numerical Analysis, Fourth
074       * Edition</u> by Burden & Faires, p 109.
075       * The variable names in this method are, as much as possible, the same as
076       * those in the book.
077       */
078      @Override
079      void runAlgorithm()
080      {
081        coefficients.clear();
082        xValues.clear();
083        
084        int n = this.size() - 1;
085        
086        BigDecimal[][] Q = new BigDecimal[n+1][n+1];
087        
088        //INPUT
089        for (int i=0; i < n+1; i++)
090        {
091          Point2D point = getPoint(i);
092          
093          xValues.add(new BigDecimal(Double.toString(point.getX()), MC));
094          Q[i][0] = new BigDecimal(Double.toString(point.getY()), MC);
095        }
096        
097        //STEP 1
098        for (int i=1; i <= n; i++)
099        {
100          for (int j=1; j <= i; j++)
101          {
102            BigDecimal numerator   = Q[i][j-1].subtract(Q[i-1][j-1]);
103            BigDecimal denominator = xValues.get(i).subtract(xValues.get(i-j));
104    
105            Q[i][j] = numerator.divide(denominator, MC);
106          }
107        }
108        
109        //STEP 2
110        for (int c=0; c <= n; c++)
111        {
112          coefficients.add(Q[c][c]);
113        }
114      }
115    }