spline.c

这是一个三次样条插值函数实现的C代码。如果插值点需要改变

/*
* 1．给定插值条件如下：
	 * i     0      1        2        3       4        5        6        7
	 * Xi  8.125    8.4      9.0      9.485    9.6      9.959    10.166     10.2
	 * Yi  0.0774   0.099    0.280    0.60     0.708    1.200     1.800     2.177
	 * 作三次样条函数插值，取第一类边界条件 Y0’=0.01087  Y7’=100
	 */
              double xi=() ;    //所要计算点的自变量值
			double []X={8.125,8.4,9.0,9.485,9.6,9.959,10.166,10.2};      //初始数据
			double []Y={0.0774,0.099,0.280,0.60,0.708,1.200,1.800,2.177};
			double Yd0=0.01087;
			double Yd7=100;
			int n=7;
double []h=new double[n+1];
			double []f=new double[n+1];
			double []lambda=new double[n+1];
			double []mu=new double[n+1];
			double []g=new double[n+1];
			double []m=new double[n+1];
			double x=0;
			double Sx=0;
          Interpolation.Get_h(h,X);            //计算 ,f[ , ]
			Interpolation.Get_f(f,Y,h);          //计算f[ , ]
			Interpolation.Get_lambda(lambda,h);  //计算 
			Interpolation.Get_mu(mu,h);          //计算 
			Interpolation.Get_g(g,lambda,mu,f);  //计算 
			double []b={2,2,2,2,2,2,2,2};
			g[1]=g[1]-lambda[1]*Yd0;
			g[n-1]=g[n-1]-mu[n-1]*Yd7;
			m[0]=Yd0;
			m[n]=Yd7;
			Interpolation.FEBS(m,lambda,b,mu,g);   //用追赶法解方程求出 
			for(int i=0;i<n+1;i++)		//查找自变量所在区间
				{
					if((xi==X[i]||x>X[i])&&x<X[i+1]) 
					{
						j=i;
						break;
					}
				}
                 
Sx=(x-X[j+1])*(x-X[j+1])*(h[j]+2*(x-X[j]))*Y[j]/(h[j]*h[j]*h[j])+(x-X[j])*(x-X[j])*(h[j]+2*(X[j+1]-x))*Y[j+1]/(h[j]*h[j]*h[j])+(x-X[j+1])*(x-X[j+1])*(x-X[j])*m[j]/(h[j]*h[j])+(x-X[j])*(x-X[j])*(x-X[j+1])*m[j+1]/(h[j]*h[j]);
		    	
string s1=" ";              //显示 ， f[ , ]， ， ， ， 
			string s2=" ";
			string s3=" ";
			string s4=" ";
			string s5=" ";
			string s6=" ";
			for(int i=0;i<n+1;i++)
			{		
				s1=s1+" "+h[i].ToString();
				s2=s2+" "+m[i].ToString();
				s3=s3+" "+f[i].ToString();
				s4=s4+" "+lambda[i].ToString();
				s5=s5+" "+mu[i].ToString();
				s6=s6+" "+g[i].ToString();
			}
			this.textBox1.Text=s1;
			this.textBox2.Text=s2;
			this.textBox3.Text=s3;
			this.textBox4.Text=s4;
			this.textBox5.Text=s5;
			this.textBox6.Text=s6;
			
	public class Interpolation
	{
		public static int n=7;
		public static void Get_h(double[]h,double[]X)            //计算 
		{
			for(int j=0;j<n;j++)
			{
				h[j]=X[j+1]-X[j];
			}
		}
		public static void Get_f(double[]f,double[]Y,double[]h)   //计算f[ , ]
		{
			for(int j=0;j<n;j++)
			{
				f[j]=(Y[j+1]-Y[j])/h[j];
			}
		}
		public static void Get_lambda(double[]lambda,double[]h) //计算 
		{
			for(int j=1;j<n;j++)
			{
				lambda[j]=h[j]/(h[j-1]+h[j]);
			}
		}
		public static void Get_mu(double[]mu,double[]h)        //计算 
		{
			for(int j=1;j<n;j++)
			{
				mu[j]=h[j-1]/(h[j-1]+h[j]);
			}
		}
		public static void Get_g(double[]g,double[]lambda,double[]mu,double[]f)// //计算 
		{
			for(int j=1;j<n;j++)
			{
				g[j]=3*(lambda[j]*f[j-1]+mu[j]*f[j]);
			}
		}
public static void FEBS(double[]x,double[]a,double[]b,double[]c,double[]f)      //forward elimination and backward substitution 追赶法计算 
		{
			int n=f.Length-2;
			double[]beta=new double[n+1];
			double[]y=new double[n+1];
			beta[1]=c[1]/b[1];
			for(int i=2;i<n;i++)
			{
				beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);
			}
			y[1]=f[1]/b[1];
			for(int i=2;i<n+1;i++)
			{
				y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);
			}
			x[n]=y[n];
			for(int i=n-1;i>0;i--)
			{
				x[i]=y[i]-beta[i]*x[i+1];
			}
		}
	}