interpolation.cs

csharp版常见数值计算源码

		 * @param t - 存放指定的插值点的值
		 * @return double 型，指定的查指点t的函数近似值f(t)
		 */
		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)