interpolation.java

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

/*
 * 进行插值的类Interpolation

 * 周长发编制
 */
package javaalgorithm.algorithm;

/**
 * 进行插值的类Interpolation

 * @author 周长发
 * @version 1.0
 */
public class Interpolation 
{
	/**
	 * 一元全区间不等距插值
	 * 
	 * @param n - 结点的个数
	 * @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)，
	 *            要求x(0)<x(1)<...<x(n-1)
	 * @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
	 *            y(i) = f(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueLagrange(int n, double[] x, double[] y, double t)
	{ 
		int i,j,k,m;
	    double z,s;

		// 初值
	    z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0];
			return(z);
		}
	    
		if (n==2)
	    { 
			z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
	        return(z);
	    }
	    
		// 开始插值
		i=0;
	    while ((x[i]<t)&&(i<n)) 
			i=i+1;
	    
		k=i-4;
	    if (k<0) 
			k=0;
	    
		m=i+3;
	    if (m>n-1) 
			m=n-1;
	    for (i=k;i<=m;i++)
	    { 
			s=1.0;
	        for (j=k;j<=m;j++)
			{
				if (j!=i) 
					// 拉格朗日插值公式
					s=s*(t-x[j])/(x[i]-x[j]);
			}

	        z=z+s*y[i];
	    }
	    
		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 - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueLagrange(int n, double x0, double xStep, double[] y, double t)
	{ 
		int i,j,k,m;
	    double z,s,xi,xj;
	    double p,q;
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    if (n==2)
	    { 
			z=(y[1]*(t-x0)-y[0]*(t-x0-xStep))/xStep;
	        return(z);
	    }
	    
		// 开始插值
		if (t>x0)
	    { 
			p=(t-x0)/xStep; 
			i=(int)p; 
			q=(float)i;
	        
			if (p>q) 
				i=i+1;
	    }
	    else 
			i=0;
	    
		k=i-4;
	    if (k<0) 
			k=0;
	    
		m=i+3;
	    if (m>n-1) 
			m=n-1;
	    
		for (i=k;i<=m;i++)
	    { 
			s=1.0; 
			xi=x0+i*xStep;
	        
			for (j=k; j<=m; j++)
	        {
				if (j!=i)
	            { 
					xj=x0+j*xStep;
	                // 拉格朗日插值公式
					s=s*(t-xj)/(xi-xj);
	            }
			}

	        z=z+s*y[i];
	    }
	    
		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 - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueLagrange3(int n, double[] x, double[] y, double t)
	{ 
		int i,j,k,m;
	    double z,s;
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    if (n==2)
	    { 
			z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
	        return(z);
	    }
	    
		// 开始插值
		if (t<=x[1]) 
		{ 
			k=0; 
			m=2;
		}
	    else if (t>=x[n-2]) 
		{ 
			k=n-3; 
			m=n-1;
		}
	    else
	    { 
			k=1; 
			m=n;
	        while (m-k!=1)
	        { 
				i=(k+m)/2;
	            
				if (t<x[i-1]) 
					m=i;
	            else 
					k=i;
	        }
	        
			k=k-1; 
			m=m-1;
	        
			if (Math.abs(t-x[k])<Math.abs(t-x[m])) 
				k=k-1;
	        else 
				m=m+1;
	    }
	    
		z=0.0;
	    for (i=k;i<=m;i++)
	    { 
			s=1.0;
	        for (j=k;j<=m;j++)
	        {
				if (j!=i) 
	                // 抛物线插值公式
					s=s*(t-x[j])/(x[i]-x[j]);
			}

	        z=z+s*y[i];
	    }
	    
		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 - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueLagrange3(int n, double x0, double xStep, double[] y, double t)
	{ 
		int i,j,k,m;
	    double z,s,xi,xj;
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    if (n==2)
	    { 
			z=(y[1]*(t-x0)-y[0]*(t-x0-xStep))/xStep;
	        return(z);
	    }
	    
		// 开始插值
		if (t<=x0+xStep) 
		{ 
			k=0; 
			m=2;
		}
	    else if (t>=x0+(n-3)*xStep) 
		{ 
			k=n-3; 
			m=n-1;
		}
	    else
	    { 
			i=(int)((t-x0)/xStep)+1;

	        if (Math.abs(t-x0-i*xStep)>=Math.abs(t-x0-(i-1)*xStep))
	        { 
				k=i-2; 
				m=i;
			}
	        else 
			{
				k=i-1; 
				m=i+1;
			}
	    }
	    
		z=0.0;
	    for (i=k;i<=m;i++)
	    { 
			s=1.0; 
			xi=x0+i*xStep;

	        for (j=k;j<=m;j++)
	        {
				if (j!=i)
	            { 
					xj=x0+j*xStep; 
	                // 抛物线插值公式
					s=s*(t-xj)/(xi-xj);
				}
			}

	        z=z+s*y[i];
	    }
	    
		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 - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValuePqs(int n, double[] x, double[] y, double t)
	{ 
		int i,j,k,m,l;
	    double z,h;
	    double[] b = new double[8];
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    
		// 连分式插值
		if (n<=8) 
		{ 
			k=0; 
			m=n;
		}
	    else if (t<x[4]) 
		{ 
			k=0; 
			m=8;
		}
	    else if (t>x[n-5]) 
		{ 
			k=n-8; 
			m=8;
		}
	    else
	    { 
			k=1; 
			j=n;
	        
			while (j-k!=1)
	        { 
				i=(k+j)/2;
	            if (t<x[i-1]) 
					j=i;
	            else 
					k=i;
	        }
	        
			k=k-4; 
			m=8;
	    }
	    
		b[0]=y[k];
	    for (i=2;i<=m;i++)
	    { 
			h=y[i+k-1]; 
			l=0; 
			j=1;
	        
			while ((l==0)&&(j<=i-1))
	        { 
				if (Math.abs(h-b[j-1])+1.0==1.0) 
					l=1;
	            else 
					h=(x[i+k-1]-x[j+k-1])/(h-b[j-1]);
	              
				j=j+1;
	        }
	        
			b[i-1]=h;
	        
			if (l!=0) 
				b[i-1]=1.0e+35;
	    }
	    
		z=b[m-1];
	    for (i=m-1;i>=1;i--) 
			z=b[i-1]+(t-x[i+k-1])/z;
	    
		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 - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValuePqs(int n, double x0, double xStep, double[] y, double t)
	{ 
		int i,j,k,m,l;
	    double z,hh,xi,xj;
	    double[] b = new double[8];
	    
		// 初值
		z=0.0;
	    
		// 特例处理
		if (n<1) 
			return(z);
	    if (n==1) 
		{ 
			z=y[0]; 
			return(z);
		}
	    
		// 连分式插值
		if (n<=8) 
		{ 
			k=0; 
			m=n;
		}
	    else if (t<(x0+4.0*xStep)) 
		{ 
			k=0; 
			m=8;
		}
	    else if (t>(x0+(n-5)*xStep)) 
		{ 
			k=n-8; 
			m=8;
		}
	    else 
		{ 
			k=(int)((t-x0)/xStep)-3; 
			m=8;
		}
	    
		b[0]=y[k];
	    for (i=2;i<=m;i++)
	    { 
			hh=y[i+k-1]; 
			l=0; 
			j=1;
	        
			while ((l==0)&&(j<=i-1))
	        { 
				if (Math.abs(hh-b[j-1])+1.0==1.0) 
					l=1;
	            else
	            { 
					xi=x0+(i+k-1)*xStep;
	                xj=x0+(j+k-1)*xStep;
	                hh=(xi-xj)/(hh-b[j-1]);
	            }
	        
				j=j+1;
	        }

	        b[i-1]=hh;
	        if (l!=0) 
				b[i-1]=1.0e+35;
	    }
	    
		z=b[m-1];
	    for (i=m-1;i>=1;i--)
			z=b[i-1]+(t-(x0+(i+k-1)*xStep))/z;
	    
		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 dy - 一维数组，长度为n，存放给定的n个结点的函数导数值y'(i)，
	 *             y'(i) = f'(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */
	public static double getValueHermite(int n, double[] x, double[] y, double[] dy, double t)
	{ 
		int i,j;
	    double z,p,q,s;
	    
		// 初值
		z=0.0;
	    
		// 循环插值
		for (i=1;i<=n;i++)
	    { 
			s=1.0;
	        
			for (j=1;j<=n;j++)
			{
				if (j!=i) 
					s=s*(t-x[j-1])/(x[i-1]-x[j-1]);
			}

	        s=s*s;
	        p=0.0;
	        
			for (j=1;j<=n;j++)
	        {
				if (j!=i) 
					p=p+1.0/(x[i-1]-x[j-1]);
			}

	        q=y[i-1]+(t-x[i-1])*(dy[i-1]-2.0*y[i-1]*p);
	        z=z+q*s;
	    }
	    
		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 dy - 一维数组，长度为n，存放给定的n个结点的函数导数值y'(i)，
	 *             y'(i) = f'(x(i)), i=0,1,...,n-1
	 * @param t - 存放指定的插值点的值
	 * @return double 型，指定的查指点t的函数近似值f(t)
	 */