interpolation.inl

本程序是在CCS环境下运行的

// Interpolation.inl		插值函数定义(实现)头文件
// Ver 1.0.0.0
// 版权所有(C)  何渝(HE Yu) 2002

// 最后修改: 2002.5.31.

#ifndef _INTERPOLATION_INL	//避免多次编译
#define _INTERPOLATION_INL

#include <valarray>			//数组模板类标准头文件
#include "Matrix.h"    	    //矩阵类头文件
#include "comm.h"			//公共头文件

//一元全区间不等距插值
template <class _Ty>
_Ty Interpolation1VariableNotIsometry(valarray<_Ty>& x, 
											valarray<_Ty>& y, _Ty t)
{
	int i,j,k,m;
    _Ty z(0), s;
    
	int n = x.size();	//插值点个数(x数组元素个数)

	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++;
    
	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);
}


//一元全区间等距插值
template <class _Ty>
_Ty Interpolation1VariableIsometry(_Ty x0, _Ty h, valarray<_Ty>& y, _Ty t)
{
	int i, j, k, m;
    _Ty z(0), s, xi, xj, p, q;
    
	int n =  y.size();				//等距结点的个数

    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 - h)) / h;
		return(z);
	}
    
	if(t > x0)
	{ 
		p = (t - x0) / h;
		i = (int)p; 
		q = (float)i;
        if(p > q) i++;
	}

    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 * h;

        for(j = k; j<=m; j++)
			if(j != i)
            {
				xj = x0 + j * h;
				s = s * (t - xj) / (xi - xj);
            }
			z = z + s * y[i];
	}
    return(z);
}

//一元三点不等距插值
template <class _Ty>
_Ty Interpolation1Variable3PointsNotIsometry(valarray<_Ty>& x, 
											valarray<_Ty>& y, _Ty t)
{
	int i, j, k, m;
	_Ty z(0.0), s;

	int n =  x.size();	//给定不等距结点的个数
	
	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(Abs(t - x[k]) < 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);
}


//一元三点等距插值
template <class _Ty>
_Ty Interpolation1Variable3PointsIsometry(_Ty x0, _Ty h, 
											valarray<_Ty>& y, _Ty t)
{
	int i, j, k, m;
    _Ty z(0.0), s, xi, xj;

	int n =  y.size();		//给定等距结点的个数

    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 - h)) / h;
		return(z);
	}
    
	if(t <= (x0 + h))
	{ 
		k = 0;
		m = 2;
	}
    
	else if(t >= (x0+(n-3)*h))
	{
		k = n -3 ; 
		m = n - 1;
	}
    
	else
	{
		i = (int)((t - x0) / h) + 1;
        
		if(Abs(t - x0 - i * h) >= Abs(t - x0 - (i - 1) * h))
		{
			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 * h;
        for(j = k; j <= m; j++)
			if(j != i)
            {
				xj = x0 + j * h;
				s = s * (t - xj) / (xi - xj);
			}
        z = z + s * y[i];
	}
 
	return(z);
}


//连分式不等距插值
template <class _Ty>
_Ty InterpolationFractionNotIsometry(valarray<_Ty>& x, 
												valarray<_Ty>& y, _Ty t)
{
	int i,j,k,m,l;
	_Ty z(0), h, b[8];

	int n = x.size();	//给定不等距结点的个数
		
	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((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);
}


//连分式等距插值
template <class _Ty>
_Ty InterpolationFractionIsometry(_Ty x0, _Ty h, valarray<_Ty>& y, _Ty t)
{
	int i,j,k,m,l;
	_Ty z(0.0), hh, xi, xj, b[8];
	
	int n = y.size();		//给定等距结点的个数
	
	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 * h)) 
	{
		k = 0; 
		m = 8;
	}
	
	else if(t > (x0 + (n - 5) * h)) 
	{
		k = n - 8;
		m = 8;
	}
	
	else
	{
		k = (int)((t - x0) / h) - 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((Abs(hh - b[j - 1]) + 1.0) == 1.0 )
			{
				l=1;
			}
			
			else
			{
				xi = x0 + (i + k - 1) * h;
				xj = x0 + (j + k - 1) * h;
				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) * h)) / z;
	}

	return(z);
}

//埃尔米特不等距插值
template <class _Ty>
_Ty InterpolationHermiteNotIsometry(valarray<_Ty>& x, 
							valarray<_Ty>& y, valarray<_Ty>& dy, _Ty t)
{
	int i,j;
	_Ty z(0.0), p, q, s;
	
	int n =  y.size();		//给定不等距结点的个数
	
	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);
}

//埃尔米特等距插值
template <class _Ty>
_Ty InterpolationHermiteIsometry(_Ty x0, _Ty h, 
						valarray<_Ty>& y, valarray<_Ty>& dy, _Ty t)
{ 
	int i, j;
	_Ty z(0.0), s, p, q;
	
	int n =  y.size();		//给定等距结点的个数
	
	for(i = 1; i <= n; i ++)
	{
		s = 1.0;
		q = x0 + (i - 1) * h;
		
		for(j = 1; j <= n; j ++)
		{
			p = x0 + (j - 1) * h;
			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) * h));
			}
		}
		
		q = y[i - 1] + (t - q) * (dy[i - 1] - 2.0 * y[i - 1] * p);
		z = z + q * s;
		
	}
	
	return(z);
}

//埃特金不等距逐步插值
template <class _Ty>
_Ty InterpolationAitkenNotIsometry(valarray<_Ty>& x, 
									valarray<_Ty>& y, _Ty t, _Ty eps)
{
	int i,j,k,m,l;
	_Ty z(0), xx[10], yy[10];
	
	int n =  y.size();		//给定不等距结点的个数
	
	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(Abs(t - x[l - 1]) > 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)) && (Abs(yy[i] - yy[i - 1]) > eps));

	return(z);
}

//埃特金等距逐步插值
template <class _Ty>
_Ty InterpolationAitkenIsometry(_Ty x, _Ty h, 
									valarray<_Ty>& y, _Ty t, _Ty eps)
{
	int i, j, k, m, l;
	_Ty z(0), xx[10], yy[10];

	int n =  y.size();	//给定等距结点的个数

	if (n < 1)
		return z;
	if (n == 1)
	{ 
		z = y[0];
		return z;
	}

	m = 10;
	if (m > n)
		m = n;
	if (t <= x)
		k = 1;
	else
	{
		if (t >= (x + (n - 1) * h))
			k = n;
		else
		{
			k = 1;
			j = n;
			while ((k - j != 1) && (k - j != -1))
			{
				l = (k + j) / 2;
				if (t < (x + (l - 1) * h))
					j = l;
				else
					k = l;
			}
			if (Abs(t - (x + (l - 1) * h)) > Abs(t - (x + (j - 1) * h)))
				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) * h;
		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) && (Abs(yy[i] - yy[i - 1]) > eps));

	return z;
}

//光滑不等距插值
template <class _Ty>
void InterpolationSmoothNotIsometry(valarray<_Ty>& x, valarray<_Ty>& y, 
										int k, _Ty t, valarray<_Ty>& s)
{
	int kk, m, l;
	_Ty u[5], p, q;

	int n =  y.size();	//给定不等距结点的个数
	
	for(m=0; m<5; m++) s[m] = 0.0;	

	if(n < 1) goto END;
	
	if(n == 1)
	{ 
		s[0] = y[0];
		s[4] = y[0];
		
		goto END;
	}
	
	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]);
		}
		
		goto END;
	}
	
	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))