interpolate.cpp

该程序用于求解科学与工程计算上的数值运算

	// 循环插值
	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);
}

//////////////////////////////////////////////////////////////////////
// 埃尔米特等距插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x0 - 存放等距n个结点中第一个结点的值
// 3. double xStep - 等距结点的步长
// 4. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 4. double dy[] - 一维数组，长度为n，存放给定的n个结点的函数导数值y'(i)，
//                 y'(i) = f'(x(i)), i=0,1,...,n-1
// 5. double t - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::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);
}

//////////////////////////////////////////////////////////////////////
// 埃特金不等距逐步插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x[] - 一维数组，长度为n，存放给定的n个结点的值x(i)
// 3. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 4. double t - 存放指定的插值点的值
// 5. double eps - 控制精度参数，默认值为0.000001
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueAitken(int n, double x[], double y[], double t, double eps /*= 0.000001*/)
{ 
	int i,j,k,m,l;
    double z,xx[10],yy[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 (fabs(t-x[l-1])>fabs(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)&&(fabs(yy[i]-yy[i-1])>eps));
    
	return(z);
}

//////////////////////////////////////////////////////////////////////
// 埃特金等距逐步插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x0 - 存放等距n个结点中第一个结点的值
// 3. double xStep - 等距结点的步长
// 4. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 5. double t - 存放指定的插值点的值
// 6. double eps - 控制精度参数，默认值为0.000001
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueAitken(int n, double x0, double xStep, double y[], double t, double eps /*= 0.000001*/)
{ 
	int i,j,k,m,l;
    double z,xx[10],yy[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 (fabs(t-(x0+(l-1)*xStep))>fabs(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)&&(fabs(yy[i]-yy[i-1])>eps));
    
	return(z);
}

//////////////////////////////////////////////////////////////////////
// 光滑不等距插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x[] - 一维数组，长度为n，存放给定的n个结点的值x(i)
// 3. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 4. double t - 存放指定的插值点的值
// 5. double s[] - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，
//					s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）
// 6. int k - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueAkima(int n, double x[], double y[], double t, double s[], int k /*= -1*/)
{ 
	int kk,m,l;
    double u[5],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])/(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]=fabs(u[3]-u[2]);
    s[1]=fabs(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]=fabs(u[3]-u[4]);
    s[1]=fabs(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];
}

//////////////////////////////////////////////////////////////////////
// 光滑等距插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x0 - 存放等距n个结点中第一个结点的值
// 3. double xStep - 等距结点的步长
// 4. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 5. double t - 存放指定的插值点的值
// 5. double s[] - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，
//					s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）
// 6. int k - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueAkima(int n, double x0, double xStep, double y[], double t, double s[], int k /*= -1*/)
{ 
	int kk,m,l;
    double u[5],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]=fabs(u[3]-u[2]);
    s[1]=fabs(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]=fabs(u[3]-u[4]);