interpolation.cs

csharp版常见数值计算源码

/*
 * 进行插值的类Interpolation

 * 周长发编制
 */

using System;

namespace 土地适宜评价系统.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