实现采用三弯矩法的三次样条插值

/**
 * @(#)work2.java
 *
 * work2 application 分段二次多项式插值法
 *
 * @author 503所 马灵霞
 * @version 1.00 2008/12/4
 */
import java.awt.Graphics2D;
import java.awt.geom.Line2D;
import java.awt.Graphics;
import java.awt.Container;
import javax.swing.JFrame;
import java.awt.Toolkit;
import java.awt.Image;
import java.awt.Dimension;
import javax.swing.JPanel;

class Method{
	static double sanWanJu(double x){
		double vectorX[]={0,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00};
   		double vectorY[]={0.00,5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00};
   		double dcY0=1.12;
   		double dcYn=-0.22;
   		double vectorH[]=new double[vectorX.length];
   		double vectorA[]=new double[vectorX.length-1];
   		double vectorR[]=new double[vectorX.length-1];
   		double vectorB[]=new double[vectorX.length];
   		double matrixA[][]=new double[vectorX.length][vectorX.length];
   		double vectorM[]=new double[vectorX.length];
   		vectorH[0]=0;
		for(int i=1;i<=vectorH.length-1;i++){
			vectorH[i]=vectorX[i]-vectorX[i-1];
		}//求 vectorH
		
   		vectorA[0]=1;
   		for(int i=1;i<=vectorA.length-1;i++){
   			vectorA[i]=vectorH[i+1]/(vectorH[i]+vectorH[i+1]);
   			vectorR[i-1]=1-vectorA[i];
   		}
   		vectorR[vectorR.length-1]=1;//vectorA,vectorR
   		vectorB[0]=-6/vectorH[1]*((vectorY[1]-vectorY[0])/vectorH[1]-dcY0);
   		for(int i=1;i<=vectorB.length-2;i++){
   			vectorB[i]=6/(vectorH[i]+vectorH[i+1])*((vectorY[i+1]-vectorY[i])/vectorH[i+1]-(vectorY[i]-vectorY[i-1])/vectorH[i]);
   		}
   		vectorB[12]=-6/vectorH[12]*(dcYn-(vectorY[12]-vectorY[11])/vectorH[12]);//vectorB

   		for(int i=0;i<=matrixA.length-1;i++){
   			matrixA[i][i]=2;
   		}
   		for(int i=0;i<=matrixA.length-2;i++){
   			matrixA[i][i+1]=vectorA[i];
   			matrixA[i+1][i]=vectorR[i];
   		}//matrixA

   		//用追赶 法 求解方程 matrixA*vectorM=vectorB;
   		double matrixL[][]=new double[matrixA.length][matrixA.length];
   		double matrixU[][]=new double[matrixA.length][matrixA.length];
   		double y[]=new double[matrixA.length];
   		matrixL[0][0]=matrixA[0][0];
   		for(int i=0;i<=matrixA.length-2;i++){
   			matrixU[i][i+1]=matrixA[i][i+1]/matrixL[i][i];
   			matrixL[i+1][i+1]=matrixA[i+1][i+1]-matrixA[i+1][i]*matrixU[i][i+1];
   		}
   		y[0]=vectorB[0]/matrixL[0][0];
   		for(int i=1;i<=matrixA.length-1;i++){
   			y[i]=(vectorB[i]-matrixA[i][i-1]*y[i-1])/matrixL[i][i];
   		}
   		vectorM[matrixA.length-1]=y[matrixA.length-1];
   		for(int i=matrixA.length-2;i>=0;i--){
   			vectorM[i]=y[i]-matrixU[i][i+1]*vectorM[i+1];
   		}
   		
   		//for(int i=0;i<=x.length-1;i++){
   			if(x>vectorX[0]&&x<vectorX[1]){
   				return vectorM[0]/6/vectorH[1]*Math.pow((vectorX[1]-x),3)+vectorM[1]/6/vectorH[1]*Math.pow((x-vectorX[0]),3)+(vectorY[0]/vectorH[1]-vectorM[0]*vectorH[1]/6)*(vectorX[1]-x)+(vectorY[1]/vectorH[1]-vectorM[1]*vectorH[1]/6)*(x-vectorX[0]);
   			}
   			if(x>=vectorX[1]&&x<vectorX[2]){
   				return vectorM[1]/6/vectorH[2]*Math.pow((vectorX[2]-x),3)+vectorM[2]/6/vectorH[2]*Math.pow((x-vectorX[1]),3)+(vectorY[1]/vectorH[2]-vectorM[1]*vectorH[2]/6)*(vectorX[2]-x)+(vectorY[2]/vectorH[2]-vectorM[2]*vectorH[2]/6)*(x-vectorX[1]);
   			}
   			if(x>=vectorX[2]&&x<vectorX[3]){
   				return vectorM[2]/6/vectorH[3]*Math.pow((vectorX[3]-x),3)+vectorM[3]/6/vectorH[3]*Math.pow((x-vectorX[2]),3)+(vectorY[2]/vectorH[3]-vectorM[2]*vectorH[3]/6)*(vectorX[3]-x)+(vectorY[3]/vectorH[3]-vectorM[3]*vectorH[3]/6)*(x-vectorX[2]);
   			}
   			if(x>=vectorX[3]&&x<vectorX[4]){
   				return vectorM[3]/6/vectorH[4]*Math.pow((vectorX[4]-x),3)+vectorM[4]/6/vectorH[4]*Math.pow((x-vectorX[3]),3)+(vectorY[3]/vectorH[4]-vectorM[3]*vectorH[4]/6)*(vectorX[4]-x)+(vectorY[4]/vectorH[4]-vectorM[4]*vectorH[4]/6)*(x-vectorX[3]);
   			}
   			if(x>=vectorX[4]&&x<vectorX[5]){
   				return vectorM[4]/6/vectorH[5]*Math.pow((vectorX[5]-x),3)+vectorM[5]/6/vectorH[5]*Math.pow((x-vectorX[4]),3)+(vectorY[4]/vectorH[5]-vectorM[4]*vectorH[5]/6)*(vectorX[5]-x)+(vectorY[5]/vectorH[5]-vectorM[5]*vectorH[5]/6)*(x-vectorX[4]);
   			}
   			if(x>=vectorX[5]&&x<vectorX[6]){
   				return vectorM[5]/6/vectorH[6]*Math.pow((vectorX[6]-x),3)+vectorM[6]/6/vectorH[6]*Math.pow((x-vectorX[5]),3)+(vectorY[5]/vectorH[6]-vectorM[5]*vectorH[6]/6)*(vectorX[6]-x)+(vectorY[6]/vectorH[6]-vectorM[6]*vectorH[6]/6)*(x-vectorX[5]);
   			}
   			if(x>=vectorX[6]&&x<vectorX[7]){
   				return vectorM[6]/6/vectorH[7]*Math.pow((vectorX[7]-x),3)+vectorM[7]/6/vectorH[7]*Math.pow((x-vectorX[6]),3)+(vectorY[6]/vectorH[7]-vectorM[6]*vectorH[7]/6)*(vectorX[7]-x)+(vectorY[7]/vectorH[7]-vectorM[7]*vectorH[7]/6)*(x-vectorX[6]);
   			}
   			if(x>=vectorX[7]&&x<vectorX[8]){
   				return vectorM[7]/6/vectorH[8]*Math.pow((vectorX[8]-x),3)+vectorM[8]/6/vectorH[8]*Math.pow((x-vectorX[7]),3)+(vectorY[7]/vectorH[8]-vectorM[7]*vectorH[8]/6)*(vectorX[8]-x)+(vectorY[8]/vectorH[8]-vectorM[8]*vectorH[8]/6)*(x-vectorX[7]);
   			}
   			if(x>=vectorX[8]&&x<vectorX[9]){
   				return vectorM[8]/6/vectorH[9]*Math.pow((vectorX[9]-x),3)+vectorM[9]/6/vectorH[9]*Math.pow((x-vectorX[8]),3)+(vectorY[8]/vectorH[9]-vectorM[8]*vectorH[9]/6)*(vectorX[9]-x)+(vectorY[9]/vectorH[9]-vectorM[9]*vectorH[9]/6)*(x-vectorX[8]);
   			}
   			if(x>=vectorX[9]&&x<vectorX[10]){
   				return vectorM[9]/6/vectorH[10]*Math.pow((vectorX[10]-x),3)+vectorM[10]/6/vectorH[10]*Math.pow((x-vectorX[9]),3)+(vectorY[9]/vectorH[10]-vectorM[9]*vectorH[10]/6)*(vectorX[10]-x)+(vectorY[10]/vectorH[10]-vectorM[10]*vectorH[10]/6)*(x-vectorX[9]);
   			}
   			if(x>=vectorX[10]&&x<vectorX[11]){
   				return vectorM[10]/6/vectorH[11]*Math.pow((vectorX[11]-x),3)+vectorM[11]/6/vectorH[11]*Math.pow((x-vectorX[10]),3)+(vectorY[10]/vectorH[11]-vectorM[10]*vectorH[11]/6)*(vectorX[11]-x)+(vectorY[11]/vectorH[11]-vectorM[11]*vectorH[11]/6)*(x-vectorX[10]);
   			}
   			if(x>=vectorX[11]&&x<vectorX[12]){
   				return vectorM[11]/6/vectorH[12]*Math.pow((vectorX[12]-x),3)+vectorM[12]/6/vectorH[12]*Math.pow((x-vectorX[11]),3)+(vectorY[11]/vectorH[12]-vectorM[11]*vectorH[12]/6)*(vectorX[12]-x)+(vectorY[12]/vectorH[12]-vectorM[12]*vectorH[12]/6)*(x-vectorX[11]);
   			}
   			else{
   				return 0;
   			}
   		//}
	}
} 
class BlankFrame extends JFrame{
	private static final int WIDTH=1200;
	private static final int HEIGHT=200;
	public BlankFrame(){
		Container con=getContentPane();
		setSize(WIDTH,HEIGHT);
		setTitle("三次样条插值");
		setLocation(10,10);
		StringPanel panel=new StringPanel();
		con.add(panel);
	}
}
class StringPanel extends JPanel{
	public static final double WIDTH=1200;
	public static final double HEIGHT=200;
	public void paintComponent(Graphics g){
		super.paintComponent(g);
		Graphics2D g2=(Graphics2D)g;
		double vectorX[]={0,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00};
    	double vectorF[]={0.00,5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00};
		Line2D.Double lineX=new Line2D.Double(0,100,1000.00,100.00);
		g2.draw(lineX);
		g2.drawString("x",1000,90);
		Line2D.Double lineY=new Line2D.Double(0,100,0,0);
		g2.draw(lineY);
		g2.drawString("y",10,10);
    	for(int i=0;i<=11;i++){
    		g2.drawString("。",(float)vectorX[i]*5,100-(float)vectorF[i]*5);
    	}
    	double i=0;
    	do{
    		g2.drawString(".",(float)i*5,100-(float)Method.sanWanJu(i)*5);
    		i=i+0.2;
    	}while(i<=190);
	}
}

public class work4 {
    
    public static void main(String[] args) {
    	System.out.println("用三次样条插值法细化后的全部数据为：");
    	for(int i=0;i<=190;i++){
    		System.out.println("坐标x为	"+i+"	时，y为：	"+Method.sanWanJu(i));
    	}
    	BlankFrame frame=new BlankFrame();
    	frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
    	frame.show();

    }
}