三次样条插值--自然边界（Java）

通过Java math3中矩阵的运算实现算法中相关的矩阵操作，可以省下不少事，先附上代码，最近要考试了，写的时间不多，后面再完善吧。

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;

import java.text.DecimalFormat;

public class Sanci {
    final int max = 10;
    double[] h = new double[max];
    double[] b = new double[max-1];
    double[] m = new double[max-1];
    double[] n = new double[max-1];//步长，常量2，μ，λ
    double[] D = new double[max+1];

    double s; //函数结果s（x)

    //求出步长h, μ值，λ值；
    public void step(double[] x){
        for (int i = 0; i < max - 1; i++)
            b[i] = 2;
        for (int i = 0; i < max; i++)
        {
            h[i] = x[i + 1] - x[i];
        }
        for (int k = 1; k < max - 1; k++)
        {
            m[k] = h[k - 1] / (h[k - 1] + h[k]);//μ值
            n[k] = 1 - m[k];//λ值
        }

       /* System.out.println("h的值：");
        for (int i =0;i<max;i++)
            System.out.println(h[i]);

        System.out.println("m的值(μ值)：");
        for (int i =1;i<max-1;i++)
            System.out.println(m[i]);

        System.out.println("n的值(λ值)：");
        for (int i =1;i<max-1;i++)
            System.out.println(n[i]);*/
    }

    //计算D
    public RealMatrix conect(double[] y){
        D[0] = 0;
        D[max] = 0;
        double[] D1 = new double[max-1];
        RealMatrix MD = null;
        for (int i = 1; i < max; i++)
        {
            D[i] = (6 / (h[i - 1] + h[i]) * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]));
            D1[i-1] =D[i];  //去掉D0和Dmax
        }
        MD = new Array2DRowRealMatrix(D1);   //将数组转成矩阵
        //System.out.println("D的值："+MD);
        return MD;
    }

    //初始化矩阵A；
   public RealMatrix  det(){
        //初始化
        double[][] A =new double[max-1][max-1]; //追赶法矩阵

        for (int i = 0; i < max -1; i++){
            for (int j = 0; j < max-1; j++){
                if (i == j){//对角线为2
                    A[i][i] = b[i];
                }
                else if (j - i == 1)
                {
                    A[i][j] = n[i + 1];
                }
                else if (i - j == 1)
                {
                    A[i][j] = m[j + 1];
                }
                else {              //其余为0
                    A[i][j] = 0;
                }
            }
        }
        RealMatrix matrix = new Array2DRowRealMatrix(A);
        return matrix;
    }

    //求逆函数
    public static RealMatrix inverseMatrix(RealMatrix A) {
        RealMatrix result = new LUDecomposition(A).getSolver().getInverse();
        return result;
    }

    //根据s（x）公式，分别计算出每个段的插值
    public String get(double[] h, double[] M, double[] x, double[] y, double l){
        double s = 0;
        DecimalFormat df = new DecimalFormat("0.0000");
        String s1 = null;
        for (int i = 0; i < max; i++)
        {
            if (l >=x[i] && l <= x[i + 1])
            {
                s = (M[i] * (Math.pow((x[i + 1] - l), 3))) / (6 * h[i]) +
                        (M[i + 1] * (Math.pow((l - x[i]), 3))) / (6 * h[i]) +
                        (y[i] - (M[i] * h[i] * h[i] / 6)) * (x[i + 1] - l) / h[i] +
                        (y[i + 1] - (M[i + 1] * h[i] * h[i] / 6)) * (l - x[i]) / h[i];
            }

        }
        s1 = df.format(s);
        return s1;
    }
    //输出矩阵
    public boolean print1(double[] array){
      int count = 0;
      for (int n=0;n<array.length;n++){
          /*   for (int q=0;q<array.length;q++){*/
          System.out.print(array[n]+"  ");
          count++;
          if (count%9==0)
              System.out.println("\n");
          /* }*/
      }
      return true;
  }
    //将二维数组转成一维
    public double[]  change(double[][] M){
        //将二维数组转成一维
        double[] m;
        int len = 0;
        // 计算一维数组长度
        for (double[] element : M) {
            len += element.length;
        }
        // 复制元素
        m = new double[len];
        int index = 0;
        for (double[] element : M) {
            for (double element2 : element) {
                m[index++] = element2;
            }
        }
        return m;
    }
    //重新加回M0，M1
    public double[] xxx(double[] m){
        double[] m1 = new double[m.length+2];
        m1[0] =0;
        m1[m.length+1]=0;
        for (int i=0;i<m.length;i++){
          m1[i+1] = m[i];
        }
        return m1;
    }



    public static void main(String[] args) {
        Sanci sc = new Sanci();
        double[] h = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
        double[] x = { -5.0000,-4.0000,-3.0000,-2.0000,-1.0000, 0 ,1.0000,2.0000,3.0000,4.0000,5.0000};//自变量
        double[] y = { -0.0000,-0.0002,-0.0004,-0.0366,-0.3679, 0 ,0.3679,0.0366,0.0004,0.0002,0.0000};//因变量
        double[] x0 ={ -5.0000,-4.5000,-4.0000,-3.5000,-3.0000,-2.5000,-2.0000,
                -1.5000,-1.0000,-0.5000,0,0.5000,1.0000,1.5000,
                2.0000,2.5000,3.0000,3.5000,4.0000,4.5000,5.0000};

        sc.step(x);
        //求M
        RealMatrix MD = sc.conect(y); //D矩阵
        RealMatrix  A = sc.det();//A矩阵
        //矩阵求逆
        RealMatrix inversetestMatrix = inverseMatrix(A);
        //A矩阵的逆乘以D
        RealMatrix MM =  inversetestMatrix.multiply(MD);
        //System.out.println("两个矩阵相乘后的结果为：\t"+ MM);
        //将M矩阵转换为M数组
        double[][] M = MM.getData();
        //将二维数组转成一维
        double[] m = sc.change(M);
        double[] m1 =  sc.xxx(m);
        String s;
        for(int p=0;p<x0.length;p++){
           s = sc.get(h,m1,x,y,x0[p]);
            System.out.println("S:"+s);
        }

    }
}