C++ to Java Code Conversion: Spline Implementation
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
}
}
原文地址: http://www.cveoy.top/t/topic/fA7L 著作权归作者所有。请勿转载和采集!