interpolation.java

实现用于工程计算的的数学方法

	public static double getValueHermite(int n, double x0, double xStep, double[] y, double[] dy, double t)
	{ 
		int i,j;
	    double z,s,p,q;
	    
		// 初值
		z=0.0;
	    
		// 循环插值
		for (i=1;i<=n;i++)
	    { 
			s=1.0; 
			q=x0+(i-1)*xStep;
	        
			for (j=1;j<=n;j++)
	        { 
				p=x0+(j-1)*xStep;
	            if (j!=i) 
					s=s*(t-p)/(q-p);
	        }
	        
			s=s*s;
	        p=0.0;
	        
			for (j=1;j<=n;j++)
	        {
				if (j!=i) 
					p=p+1.0/(q-(x0+(j-1)*xStep));
			}

	        q=y[i-1]+(t-q)*(dy[i-1]-2.0*y[i-1]*p);
	        z=z+q*s;
	    }
	    
		return(z);
	}

	/**
	 * 埃特金不等距逐步插值
	 * 
	 * @param n - 结点的个数
	 * @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)
	 * @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
	 *            y(i) = f(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @param eps - 控制精度参数
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueAitken(int n, double[] x, double[] y, double t, double eps)
	{ 
		int i,j,k,m,l=0;
	    double z;
	    double[] xx = new double[10];
	    double[] yy = new double[10];
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    
		// 开始插值
		m=10;
	    if (m>n) 
			m=n;
	    
		if (t<=x[0]) 
			k=1;
	    else if (t>=x[n-1]) 
			k=n;
	    else
	    { 
			k=1; 
			j=n;
	        
			while ((k-j!=1)&&(k-j!=-1))
	        { 
				l=(k+j)/2;
	            if (t<x[l-1]) 
					j=l;
	            else 
					k=l;
	        }
	        
			if (Math.abs(t-x[l-1])>Math.abs(t-x[j-1])) 
				k=j;
	    }
	    
		j=1; 
		l=0;
	    
		for (i=1;i<=m;i++)
	    { 
			k=k+j*l;
	        if ((k<1)||(k>n))
	        { 
				l=l+1; 
				j=-j; 
				k=k+j*l;
			}
	        
			xx[i-1]=x[k-1]; 
			yy[i-1]=y[k-1];
	        l=l+1; 
			j=-j;
	    }
	    
		i=0;
	    
		do
	    { 
			i=i+1; 
			z=yy[i];
	        
			for (j=0;j<=i-1;j++)
				z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx[i]);
	        
			yy[i]=z;
	    } while ((i!=m-1)&&(Math.abs(yy[i]-yy[i-1])>eps));
	    
		return(z);
	}

	/**
	 * 埃特金等距逐步插值
	 * 
	 * @param n - 结点的个数
	 * @param x0 - 存放等距n个结点中第一个结点的值
	 * @param xStep - 等距结点的步长
	 * @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
	 *            y(i) = f(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @param eps - 控制精度参数
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueAitken(int n, double x0, double xStep, double[] y, double t, double eps)
	{ 
		int i,j,k,m,l=0;
	    double z;
	    double[] xx = new double[10];
	    double[] yy = new double[10];
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    
		// 开始插值
		m=10;
	    if (m>n) 
			m=n;
	    
		if (t<=x0) 
			k=1;
	    else if (t>=x0+(n-1)*xStep) 
			k=n;
	    else
	    { 
			k=1; 
			j=n;
	        
			while ((k-j!=1)&&(k-j!=-1))
	        { 
				l=(k+j)/2;
	            
				if (t<x0+(l-1)*xStep) 
					j=l;
	            else 
					k=l;
	        }
	        
			if (Math.abs(t-(x0+(l-1)*xStep))>Math.abs(t-(x0+(j-1)*xStep))) 
				k=j;
	    }
	    
		j=1; 
		l=0;
	    for (i=1;i<=m;i++)
	    { 
			k=k+j*l;
	        if ((k<1)||(k>n))
	        { 
				l=l+1; 
				j=-j; 
				k=k+j*l;
			}
	        
			xx[i-1]=x0+(k-1)*xStep; 
			yy[i-1]=y[k-1];
	        l=l+1; 
			j=-j;
	    }
	    
		i=0;
	    do
	    { 
			i=i+1; 
			z=yy[i];
	        
			for (j=0;j<=i-1;j++)
				z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx[i]);
	        
			yy[i]=z;
	    } while ((i!=m-1)&&(Math.abs(yy[i]-yy[i-1])>eps));
	    
		return(z);
	}

	/**
	 * 光滑不等距插值
	 * 
	 * @param n - 结点的个数
	 * @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)
	 * @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
	 *            y(i) = f(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @param s - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，
	 *  		  s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）
	 * @param k - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueAkima(int n, double[] x, double[] y, double t, double[] s, int k)
	{ 
		int kk,m,l;
	    double p,q;
	    double[] u = new double[5];
	    
		// 初值
		s[4]=0.0; 
		s[0]=0.0; 
		s[1]=0.0; 
		s[2]=0.0; 
		s[3]=0.0;
	    
		// 特例处理
		if (n<1) 
			return s[4];
	    if (n==1) 
		{ 
			s[0]=y[0]; 
			s[4]=y[0]; 
			return s[4];
		}
	    if (n==2)
	    { 
			s[0]=y[0]; 
			s[1]=(y[1]-y[0])/(x[1]-x[0]);
	        if (k<0)
				s[4]=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
	        return s[4];
	    }
	    
		// 插值
		if (k<0)
	    { 
			if (t<=x[1]) 
				kk=0;
	        else if (t>=x[n-1]) 
				kk=n-2;
	        else
	        { 
				kk=1; 
				m=n;
	            while (((kk-m)!=1)&&((kk-m)!=-1))
	            { 
					l=(kk+m)/2;
	                if (t<x[l-1]) 
						m=l;
	                else 
						kk=l;
	            }
	            
				kk=kk-1;
	        }
	    }
	    else 
			kk=k;
	    
		if (kk>=n-1) 
			kk=n-2;
	    
		u[2]=(y[kk+1]-y[kk])/(x[kk+1]-x[kk]);
	    if (n==3)
	    { 
			if (kk==0)
	        { 
				u[3]=(y[2]-y[1])/(x[2]-x[1]);
	            u[4]=2.0*u[3]-u[2];
	            u[1]=2.0*u[2]-u[3];
	            u[0]=2.0*u[1]-u[2];
	        }
	        else
	        { 
				u[1]=(y[1]-y[0])/(x[1]-x[0]);
	            u[0]=2.0*u[1]-u[2];
	            u[3]=2.0*u[2]-u[1];
	            u[4]=2.0*u[3]-u[2];
	        }
	    }
	    else
	    { 
			if (kk<=1)
	        { 
				u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
	            if (kk==1)
	            { 
					u[1]=(y[1]-y[0])/(x[1]-x[0]);
	                u[0]=2.0*u[1]-u[2];
	                
					if (n==4) 
						u[4]=2.0*u[3]-u[2];
	                else 
						u[4]=(y[4]-y[3])/(x[4]-x[3]);
	            }
	            else
	            { 
					u[1]=2.0*u[2]-u[3];
	                u[0]=2.0*u[1]-u[2];
	                u[4]=(y[3]-y[2])/(x[3]-x[2]);
	            }
	        }
	        else if (kk>=(n-3))
	        { 
				u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
	            if (kk==(n-3))
	            { 
					u[3]=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
	                u[4]=2.0*u[3]-u[2];
	                if (n==4) 
						u[0]=2.0*u[1]-u[2];
	                else 
						u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
	            }
	            else
	            { 
					u[3]=2.0*u[2]-u[1];
	                u[4]=2.0*u[3]-u[2];
	                u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
	            }
	        }
	        else
	        { 
				u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
	            u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
	            u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
	            u[4]=(y[kk+3]-y[kk+2])/(x[kk+3]-x[kk+2]);
	        }
	    }
	    
		s[0]=Math.abs(u[3]-u[2]);
	    s[1]=Math.abs(u[0]-u[1]);
	    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
	         p=(u[1]+u[2])/2.0;
	    else 
			p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
	    
		s[0]=Math.abs(u[3]-u[4]);
	    s[1]=Math.abs(u[2]-u[1]);
	    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
	        q=(u[2]+u[3])/2.0;
	    else 
			q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
	    
		s[0]=y[kk];
	    s[1]=p;
	    s[3]=x[kk+1]-x[kk];
	    s[2]=(3.0*u[2]-2.0*p-q)/s[3];
	    s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
	    if (k<0)
	    { 
			p=t-x[kk];
	        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
	    }
	    
		return s[4];
	}

	/**
	 * 光滑等距插值
	 * 
	 * @param n - 结点的个数
	 * @param x0 - 存放等距n个结点中第一个结点的值
	 * @param xStep - 等距结点的步长
	 * @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
	 *            y(i) = f(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @param s - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，
	 *  		  s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）
	 * @param k - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueAkima(int n, double x0, double xStep, double[] y, double t, double[] s, int k)
	{ 
		int kk,m,l;
	    double[] u = new double[5];
	    double p,q;
	    
		// 初值
		s[4]=0.0; 
		s[0]=0.0; 
		s[1]=0.0; 
		s[2]=0.0; 
		s[3]=0.0;
	    
		// 特例处理
		if (n<1) 
			return s[4];
	    if (n==1) 
		{ 
			s[0]=y[0]; 
			s[4]=y[0]; 
			return s[4];
		}
	    if (n==2)
	    { 
			s[0]=y[0]; 
			s[1]=(y[1]-y[0])/xStep;
	        if (k<0)
				s[4]=(y[1]*(t-x0)-y[0]*(t-x0-xStep))/xStep;
	        return s[4];
	    }
	    
		// 插值
		if (k<0)
	    { 
			if (t<=x0+xStep) 
				kk=0;
	        else if (t>=x0+(n-1)*xStep) 
				kk=n-2;
	        else
	        { 
				kk=1; 
				m=n;
	            while (((kk-m)!=1)&&((kk-m)!=-1))
	            { 
					l=(kk+m)/2;
	                if (t<x0+(l-1)*xStep) 
						m=l;
	                else 
						kk=l;
	            }
	            
				kk=kk-1;
	        }
	    }
	    else 
			kk=k;
	    
		if (kk>=n-1) 
			kk=n-2;
	    
		u[2]=(y[kk+1]-y[kk])/xStep;
	    if (n==3)
	    { 
			if (kk==0)
	        { 
				u[3]=(y[2]-y[1])/xStep;
	            u[4]=2.0*u[3]-u[2];
	            u[1]=2.0*u[2]-u[3];
	            u[0]=2.0*u[1]-u[2];
	        }
	        else
	        { 
				u[1]=(y[1]-y[0])/xStep;
	            u[0]=2.0*u[1]-u[2];
	            u[3]=2.0*u[2]-u[1];
	            u[4]=2.0*u[3]-u[2];
	        }
	    }
	    else
	    { 
			if (kk<=1)
	        { 
				u[3]=(y[kk+2]-y[kk+1])/xStep;
	            if (kk==1)
	            { 
					u[1]=(y[1]-y[0])/xStep;
	                u[0]=2.0*u[1]-u[2];
	                if (n==4) 
						u[4]=2.0*u[3]-u[2];
	                else 
						u[4]=(y[4]-y[3])/xStep;
	            }
	            else
	            { 
					u[1]=2.0*u[2]-u[3];
	                u[0]=2.0*u[1]-u[2];
	                u[4]=(y[3]-y[2])/xStep;
	            }
	        }
	        else if (kk>=(n-3))
	        { 
				u[1]=(y[kk]-y[kk-1])/xStep;
	            if (kk==(n-3))
	            { 
					u[3]=(y[n-1]-y[n-2])/xStep;
	                u[4]=2.0*u[3]-u[2];
	                if (n==4) 
						u[0]=2.0*u[1]-u[2];
	                else 
						u[0]=(y[kk-1]-y[kk-2])/xStep;
	            }
	            else
	            { 
					u[3]=2.0*u[2]-u[1];
	                u[4]=2.0*u[3]-u[2];
	                u[0]=(y[kk-1]-y[kk-2])/xStep;
	            }
	        }
	        else
	        { 
				u[1]=(y[kk]-y[kk-1])/xStep;
	            u[0]=(y[kk-1]-y[kk-2])/xStep;
	            u[3]=(y[kk+2]-y[kk+1])/xStep;
	            u[4]=(y[kk+3]-y[kk+2])/xStep;
	        }
	    }
	    
		s[0]=Math.abs(u[3]-u[2]);
	    s[1]=Math.abs(u[0]-u[1]);
	    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
			p=(u[1]+u[2])/2.0;
	    else 
			p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
	    
		s[0]=Math.abs(u[3]-u[4]);
	    s[1]=Math.abs(u[2]-u[1]);
	    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
			q=(u[2]+u[3])/2.0;
	    else 
			q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
	    
		s[0]=y[kk];
	    s[1]=p;
	    s[3]=xStep;
	    s[2]=(3.0*u[2]-2.0*p-q)/s[3];
	    s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
	    
		if (k<0)
	    { 
			p=t-(x0+kk*xStep);