import java.util.ArrayList;

public class Spline { public enum bd_type { second_deriv, first_deriv };

private bd_type m_left;
private bd_type m_right;
private double m_left_value;
private double m_right_value;
private boolean m_force_linear_extrapolation;

private ArrayList<Double> m_x;
private ArrayList<Double> m_y;
private ArrayList<Double> m_a;
private ArrayList<Double> m_b;
private ArrayList<Double> m_c;
private double m_b0;
private double m_c0;

public void set_boundary(bd_type left, double left_value, bd_type right, double right_value, boolean force_linear_extrapolation) {
    assert m_x.size() == 0;
    m_left = left;
    m_right = right;
    m_left_value = left_value;
    m_right_value = right_value;
    m_force_linear_extrapolation = force_linear_extrapolation;
}

public void set_points(ArrayList<Double> x, ArrayList<Double> y, boolean cubic_spline) {
    assert x.size() == y.size();
    assert x.size() > 2;
    m_x = x;
    m_y = y;
    int n = x.size();
    
    // TODO: maybe sort x and y, rather than returning an error
    for (int i = 0; i < n - 1; i++) {
        assert m_x.get(i) < m_x.get(i + 1);
    }
    
    if (cubic_spline) {
        // cubic spline interpolation
        // setting up the matrix and right hand side of the equation system for the parameters b[]
        BandMatrix A = new BandMatrix(n, 1, 1);
        ArrayList<Double> rhs = new ArrayList<Double>(n);
        for (int i = 1; i < n - 1; i++) {
            A.set(i, i - 1, 1.0 / 3.0 * (x.get(i) - x.get(i - 1)));
            A.set(i, i, 2.0 / 3.0 * (x.get(i + 1) - x.get(i - 1)));
            A.set(i, i + 1, 1.0 / 3.0 * (x.get(i + 1) - x.get(i)));
            rhs.add((y.get(i + 1) - y.get(i)) / (x.get(i + 1) - x.get(i)) - (y.get(i) - y.get(i - 1)) / (x.get(i) - x.get(i - 1)));
        }
        
        // boundary conditions
        if (m_left == bd_type.second_deriv) {
            // 2 * b[0] = f''
            A.set(0, 0, 2.0);
            A.set(0, 1, 0.0);
            rhs.set(0, m_left_value);
        } else if (m_left == bd_type.first_deriv) {
            // c[0] = f', needs to be re-expressed in terms of b:
            // (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
            A.set(0, 0, 2.0 * (x.get(1) - x.get(0)));
            A.set(0, 1, 1.0 * (x.get(1) - x.get(0)));
            rhs.set(0, 3.0 * ((y.get(1) - y.get(0)) / (x.get(1) - x.get(0)) - m_left_value));
        } else {
            assert false;
        }
        
        if (m_right == bd_type.second_deriv) {
            // 2 * b[n-1] = f''
            A.set(n - 1, n - 1, 2.0);
            A.set(n - 1, n - 2, 0.0);
            rhs.set(n - 1, m_right_value);
        } else if (m_right == bd_type.first_deriv) {
            // c[n-1] = f', needs to be re-expressed in terms of b:
            // (b[n-2]+2b[n-1])(x[n-1]-x[n-2]) = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
            A.set(n - 1, n - 1, 2.0 * (x.get(n - 1) - x.get(n - 2)));
            A.set(n - 1, n - 2, 1.0 * (x.get(n - 1) - x.get(n - 2)));
            rhs.set(n - 1, 3.0 * (m_right_value - (y.get(n - 1) - y.get(n - 2)) / (x.get(n - 1) - x.get(n - 2))));
        } else {
            assert false;
        }
        
        // solve the equation system to obtain the parameters b[]
        m_b = A.lu_solve(rhs);
        
        // calculate parameters a[] and c[] based on b[]
        m_a = new ArrayList<Double>(n);
        m_c = new ArrayList<Double>(n);
        for (int i = 0; i < n - 1; i++) {
            m_a.add(1.0 / 3.0 * (m_b.get(i + 1) - m_b.get(i)) / (x.get(i + 1) - x.get(i)));
            m_c.add((y.get(i + 1) - y.get(i)) / (x.get(i + 1) - x.get(i)) - 1.0 / 3.0 * (2.0 * m_b.get(i) + m_b.get(i + 1)) * (x.get(i + 1) - x.get(i)));
        }
    } else {
        // linear interpolation
        m_a = new ArrayList<Double>(n);
        m_b = new ArrayList<Double>(n);
        m_c = new ArrayList<Double>(n);
        for (int i = 0; i < n - 1; i++) {
            m_a.add(0.0);
            m_b.add(0.0);
            m_c.add((m_y.get(i + 1) - m_y.get(i)) / (m_x.get(i + 1) - m_x.get(i)));
        }
    }
    
    // for left extrapolation coefficients
    m_b0 = (m_force_linear_extrapolation == false) ? m_b.get(0) : 0.0;
    m_c0 = m_c.get(0);
    
    // for the right extrapolation coefficients
    // f_{n-1}(x) = b * (x - x_{n-1})^2 + c * (x - x_{n-1}) + y_{n-1}
    double h = x.get(n - 1) - x.get(n - 2);
    // m_b[n-1] is determined by the boundary condition
    m_a.set(n - 1, 0.0);
    m_c.set(n - 1, 3.0 * m_a.get(n - 2) * h * h + 2.0 * m_b.get(n - 2) * h + m_c.get(n - 2));   // = f'_{n-2}(x_{n-1})
    if (m_force_linear_extrapolation == true)
        m_b.set(n - 1, 0.0);
}

public double evaluate(double x) {
    int n = m_x.size();
    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
    int idx = 0;
    for (int i = 0; i < m_x.size(); i++) {
        if (m_x.get(i) >= x) {
            break;
        }
        idx = i;
    }
    
    double h = x - m_x.get(idx);
    double interpol;
    if (x < m_x.get(0)) {
        // extrapolation to the left
        interpol = (m_b0 * h + m_c0) * h + m_y.get(0);
    } else if (x > m_x.get(n - 1)) {
        // extrapolation to the right
        interpol = (m_b.get(n - 1) * h + m_c.get(n - 1)) * h + m_y.get(n - 1);
    } else {
        // interpolation
        interpol = ((m_a.get(idx) * h + m_b.get(idx)) * h + m_c.get(idx)) * h + m_y.get(idx);
    }
    return interpol;
}

public double deriv(int order, double x) {
    assert order > 0;
    
    int n = m_x.size();
    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
    int idx = 0;
    for (int i = 0; i < m_x.size(); i++) {
        if (m_x.get(i) >= x) {
            break;
        }
        idx = i;
    }
    
    double h = x - m_x.get(idx);
    double interpol;
    if (x < m_x.get(0)) {
        // extrapolation to the left
        switch(order) {
            case 1:
                interpol = 2.0 * m_b0 * h + m_c0;
                break;
            case 2:
                interpol = 2.0 * m_b0 * h;
                break;
            default:
                interpol = 0.0;
                break;
        }
    } else if (x > m_x.get(n - 1)) {
        // extrapolation to the right
        switch(order) {
            case 1:
                interpol = 2.0 * m_b.get(n - 1) * h + m_c.get(n - 1);
                break;
            case 2:
                interpol = 2.0 * m_b.get(n - 1);
                break;
            default:
                interpol = 0.0;
                break;
        }
    } else {
        // interpolation
        switch(order) {
            case 1:
                interpol = (3.0 * m_a.get(idx) * h + 2.0 * m_b.get(idx)) * h + m_c.get(idx);
                break;
            case 2:
                interpol = 6.0 * m_a.get(idx) * h + 2.0 * m_b.get(idx);
                break;
            case 3:
                interpol = 6.0 * m_a.get(idx);
                break;
            default:
                interpol = 0.0;
                break;
        }
    }
    return interpol;
}

// Add a BandMatrix class for solving the linear system
private class BandMatrix {
    private int n, kl, ku;
    private double[][] data;
    
    public BandMatrix(int n, int kl, int ku) {
        this.n = n;
        this.kl = kl;
        this.ku = ku;
        data = new double[n][kl + ku + 1];
    }
    
    public void set(int i, int j, double value) {
        if (i >= 0 && i < n && j >= i - kl && j <= i + ku) {
            data[i][j - i + kl] = value;
        }
    }
    
    public double get(int i, int j) {
        if (i >= 0 && i < n && j >= i - kl && j <= i + ku) {
            return data[i][j - i + kl];
        }
        return 0.0;
    }
    
    public ArrayList<Double> lu_solve(ArrayList<Double> rhs) {
        // Implementation of LU decomposition and solving the linear system
        // ...
        return new ArrayList<Double>(n); // Replace with actual solution
    }
}
C++ to Java Code Conversion: Spline Implementation

原文地址: http://www.cveoy.top/t/topic/fA7L 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录