interpolate.cpp

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

//////////////////////////////////////////////////////////////////////
// Interpolate.cpp
//
// 插值类 CInterpolate 的实现代码
//
// 周长发编制, 2002/8
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Interpolate.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
// 基本构造函数
//////////////////////////////////////////////////////////////////////
CInterpolate::CInterpolate()
{
}

//////////////////////////////////////////////////////////////////////
// 析构函数
//////////////////////////////////////////////////////////////////////
CInterpolate::~CInterpolate()
{
}

//////////////////////////////////////////////////////////////////////
// 将字符串转化为结点的值
//
// 参数：
// 1. CString s - 数字和分隔符构成的字符串
// 2. int n - 结点的个数
// 3. double dblNodes[] - 一维数组，长度为n，返回结点的值
// 4. const CString& sDelim - 数字之间的分隔符，默认为空格
//
// 返回值：int 型，转换成功的结点的个数
//////////////////////////////////////////////////////////////////////
int CInterpolate::GetNodesFromString(CString s, int n, double dblNodes[], const CString& sDelim /*= " "*/)
{
	// 将结点值初始化为0
	memset(dblNodes, 0, n*sizeof(double));

	// 构造根据指定的分界符将字符串分解为不同的子串的对象
	CTokenizer tk(s, sDelim);

	CString sElement;
	
	// 分解字符串，给结点赋值
	int i = 0;
	while (tk.Next(sElement))
	{
		sElement.TrimLeft();
		sElement.TrimRight();
		double v = atof(sElement);
		dblNodes[i++] = v;

		if (i == n)
			break;
	}

	return i;
}

//////////////////////////////////////////////////////////////////////
// 一元全区间不等距插值
//
// 参数：
// 1. int n - 结点的个数
// 2. double x[] - 一维数组，长度为n，存放给定的n个结点的值x(i)，
//                 要求x(0)<x(1)<...<x(n-1)
// 3. double y[] - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
//                 y(i) = f(x(i)), i=0,1,...,n-1
// 4. double t - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::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);
}

//////////////////////////////////////////////////////////////////////
// 一元全区间等距插值
//
// 参数：
// 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 - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueLagrange(int n, double x0, double xStep, double y[], double t)
{ 
	int i,j,k,m;
    double z,s,xi,xj;
    float 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);
}

//////////////////////////////////////////////////////////////////////
// 一元三点不等距插值
//
// 参数：
// 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 - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::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 (fabs(t-x[k])<fabs(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);
}

//////////////////////////////////////////////////////////////////////
// 一元三点等距插值
//
// 参数：
// 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 - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::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 (fabs(t-x0-i*xStep)>=fabs(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);
}

//////////////////////////////////////////////////////////////////////
// 连分式不等距插值
//
// 参数：
// 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 - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValuePqs(int n, double x[], double y[], double t)
{ 
	int i,j,k,m,l;
    double z,h,b[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 (fabs(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);
}

//////////////////////////////////////////////////////////////////////
// 连分式等距插值
//
// 参数：
// 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 - 存放指定的插值点的值
//
// 返回值：double 型，指定的查指点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValuePqs(int n, double x0, double xStep, double y[], double t)
{ 
	int i,j,k,m,l;
    double z,hh,xi,xj,b[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 (fabs(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);
}

//////////////////////////////////////////////////////////////////////
// 埃尔米特不等距插值
//
// 参数：
// 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 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 x[], double y[], double dy[], double t)
{ 
	int i,j;
    double z,p,q,s;
    
	// 初值
	z=0.0;