interpolate.cpp

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

    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]=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)
    { 
		p=t-(x0+kk*xStep);
        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
    }
    
	return s[4];
}

//////////////////////////////////////////////////////////////////////
// 第一种边界条件的三次样条函数插值、微商与积分
//
// 参数：
// 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，调用时，dy(0)存放给定区间的左端点处的一阶导数值，
//                  dy(n-1)存放给定区间的右端点处的一阶导数值。返回时，存放n个给定点处的
//                  一阶导数值y'(i)，i=0,1,...,n-1
// 5. double ddy[] - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
//                  i=0,1,...,n-1
// 6. int m - 指定插值点的个数
// 7. double t[] - 一维数组，长度为m，存放m个指定的插值点的值。
//                 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一维数组，长度为m，存放m个指定的插值点处的函数值
// 9. double dz[] - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
// 10. double ddz[] - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
//
// 返回值：double 型，指定函数的x(0)到x(n-1)的定积分值
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueSpline1(int n, double x[], double y[], double dy[], double ddy[], 
					  int m, double t[], double z[], double dz[], double ddz[])
{ 
	int i,j;
    double h0,h1,alpha,beta,g,*s;
    
	// 初值
	s=new double[n];
    s[0]=dy[0]; 
	dy[0]=0.0;
    h0=x[1]-x[0];
    
	for (j=1;j<=n-2;j++)
    { 
		h1=x[j+1]-x[j];
        alpha=h0/(h0+h1);
        beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
        beta=3.0*(beta+alpha*(y[j+1]-y[j])/h1);
        dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
        s[j]=(beta-(1.0-alpha)*s[j-1]);
        s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
        h0=h1;
    }
    
	for (j=n-2;j>=0;j--)
		dy[j]=dy[j]*dy[j+1]+s[j];
    
	for (j=0;j<=n-2;j++) 
		s[j]=x[j+1]-x[j];
    
	for (j=0;j<=n-2;j++)
    { 
		h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
    }
    
	h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    
	for (i=0;i<=n-2;i++)
    { 
		h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
    }
    
	for (j=0;j<=m-1;j++)
    { 
		if (t[j]>=x[n-1]) 
			i=n-2;
        else
        { 
			i=0;
            while (t[j]>x[i+1]) 
				i=i+1;
        }
        
		h1=(x[i+1]-t[j])/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(t[j]-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
    }
    
	delete[] s;

    return(g);
}

//////////////////////////////////////////////////////////////////////
// 第二种边界条件的三次样条函数插值、微商与积分
//
// 参数：
// 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)，
//                 i=0,1,...,n-1
// 5. double ddy[] - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
//                  i=0,1,...,n-1，调用时，ddy(0)存放给定区间的左端点处的二阶导数值，
//                  ddy(n－1)存放给定区间的右端点处的二阶导数值
// 6. int m - 指定插值点的个数
// 7. double t[] - 一维数组，长度为m，存放m个指定的插值点的值。
//                 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一维数组，长度为m，存放m个指定的插值点处的函数值
// 9. double dz[] - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
// 10. double ddz[] - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
//
// 返回值：double 型，指定函数的x(0)到x(n-1)的定积分值
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueSpline2(int n, double x[], double y[], double dy[], double ddy[], 
					  int m, double t[], double z[], double dz[], double ddz[])
{ 
	int i,j;
    double h0,h1,alpha,beta,g,*s;
    
	// 初值
	s=new double[n];
    dy[0]=-0.5;
    h0=x[1]-x[0];
    s[0]=3.0*(y[1]-y[0])/(2.0*h0)-ddy[0]*h0/4.0;
    
	for (j=1;j<=n-2;j++)
    { 
		h1=x[j+1]-x[j];
        alpha=h0/(h0+h1);
        beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
        beta=3.0*(beta+alpha*(y[j+1]-y[j])/h1);
        dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
        s[j]=(beta-(1.0-alpha)*s[j-1]);
        s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
        h0=h1;
    }
    
	dy[n-1]=(3.0*(y[n-1]-y[n-2])/h1+ddy[n-1]*h1/2.0-s[n-2])/(2.0+dy[n-2]);
    for (j=n-2;j>=0;j--)
		dy[j]=dy[j]*dy[j+1]+s[j];
    
	for (j=0;j<=n-2;j++) 
		s[j]=x[j+1]-x[j];
    
	for (j=0;j<=n-2;j++)
    { 
		h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
    }
    
	h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    
	for (i=0;i<=n-2;i++)
    { 
		h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
    }
    
	for (j=0;j<=m-1;j++)
    { 
		if (t[j]>=x[n-1]) 
			i=n-2;
        else
        { 
			i=0;
            while (t[j]>x[i+1]) 
				i=i+1;
        }
        
		h1=(x[i+1]-t[j])/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(t[j]-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
    }
    
	delete[] s;

    return(g);
}

//////////////////////////////////////////////////////////////////////
// 第三种边界条件的三次样条函数插值、微商与积分
//
// 参数：
// 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)，
//                 i=0,1,...,n-1
// 5. double ddy[] - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
//                  i=0,1,...,n-1
// 6. int m - 指定插值点的个数
// 7. double t[] - 一维数组，长度为m，存放m个指定的插值点的值。
//                 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一维数组，长度为m，存放m个指定的插值点处的函数值
// 9. double dz[] - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
// 10. double ddz[] - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
//
// 返回值：double 型，指定函数的x(0)到x(n-1)的定积分值
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueSpline3(int n, double x[], double y[], double dy[], double ddy[], 
					  int m, double t[], double z[], double dz[], double ddz[])
{ 
	int i,j;
    double h0,y0,h1,y1,alpha,beta,u,g,*s;
    
	// 初值
	s=new double[n];
    h0=x[n-1]-x[n-2];
    y0=y[n-1]-y[n-2];
    dy[0]=0.0; ddy[0]=0.0; ddy[n-1]=0.0;
    s[0]=1.0; s[n-1]=1.0;

    for (j=1;j<=n-1;j++)
    { 
		h1=h0; y1=y0;
        h0=x[j]-x[j-1];
        y0=y[j]-y[j-1];
        alpha=h1/(h1+h0);
        beta=3.0*((1.0-alpha)*y1/h1+alpha*y0/h0);
        
		if (j<n-1)
        { 
			u=2.0+(1.0-alpha)*dy[j-1];
            dy[j]=-alpha/u;
            s[j]=(alpha-1.0)*s[j-1]/u;
            ddy[j]=(beta-(1.0-alpha)*ddy[j-1])/u;
        }
    }
    
	for (j=n-2;j>=1;j--)
    { 
		s[j]=dy[j]*s[j+1]+s[j];
        ddy[j]=dy[j]*ddy[j+1]+ddy[j];
    }
    
	dy[n-2]=(beta-alpha*ddy[1]-(1.0-alpha)*ddy[n-2])/
            (alpha*s[1]+(1.0-alpha)*s[n-2]+2.0);
    
	for (j=2;j<=n-1;j++)
        dy[j-2]=s[j-1]*dy[n-2]+ddy[j-1];
    
	dy[n-1]=dy[0];
    for (j=0;j<=n-2;j++) 
		s[j]=x[j+1]-x[j];
    
	for (j=0;j<=n-2;j++)
    { 
		h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
    }
    
	h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    
	for (i=0;i<=n-2;i++)
    { 
		h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
    }
    
	for (j=0;j<=m-1;j++)
    { 
		h0=t[j];
        while (h0>=x[n-1]) 
			h0=h0-(x[n-1]-x[0]);
        
		while (h0<x[0]) 
			h0=h0+(x[n-1]-x[0]);
        
		i=0;
        while (h0>x[i+1]) 
			i=i+1;
        
		u=h0;
        h1=(x[i+1]-u)/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(u-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
    }
 
	delete[] s;

    return(g);
}

//////////////////////////////////////////////////////////////////////
// 二元三点插值
//
// 参数：
// 1. int n - x方向上给定结点的点数
// 2. double x[] - 一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)
// 3. int m - y方向上给定结点的点数
// 4. double y[] - 一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)
// 5. double z[] - 一维数组，长度为n x m，存放给定的n x m个结点的函数值z(i,j)，
//                 z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
// 6. double u - 存放插值点x坐标
// 7. double v - 存放插值点y坐标
//
// 返回值：double 型，指定函数值f(u, v)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueTqip(int n, double x[], int m, double y[], double z[], double u, double v)
{ 
	int nn,mm,ip,iq,i,j,k,l;
    double b[3],h,w;
    
	// 初值
	nn=3;

	// 特例
    if (n<=3) 
	{ 
		ip=0;  
		nn=n;
	}
    else if (u<=x[1]) 
		ip=0;
    else if (u>=x[n-2]) 
		ip=n-3;
    else					
    { 
		i=1; j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
        { 
			l=(i+j)/2;
            if (u<x[l-1]) 
				j=l;
            else 
				i=l;
        }
        
		if (fabs(u-x[i-1])<fabs(u-x[j-1])) 
			ip=i-2;
        else 
			ip=i-1;
    }
    
	mm=3;
    
	if (m<=3) 
	{ 
		iq=0; 
		mm=m;
	}
    else if (v<=y[1]) 
		iq=0;
    else if (v>=y[m-2]) 
		iq=m-3;
    else
    { 
		i=1; 
		j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
        { 
			l=(i+j)/2;
            if (v<y[l-1]) 
				j=l;
            else 
				i=l;
        }
        
		if (fabs(v-y[i-1])<fabs(v-y[j-1])) 
			iq=i-2;
        else 
			iq=i-1;
    }
    
	for (i=0;i<=nn-1;i++)
    { 
		b[i]=0.0;
        for (j=0;j<=mm-1;j++)
        { 
			k=m*(ip+i)+(iq+j);
            h=z[k];
            for (k=0;k<=mm-1;k++)
            {
				if (k!=j)
					h=h*(v-y[iq+k])/(y[iq+j]-y[iq+k]);
			}

            b[i]=b[i]+h;
        }
    }
    
	w=0.0;
    for (i=0;i<=nn-1;i++)
    { 
		h=b[i];
        for (j=0;j<=nn-1;j++)
        {
			if (j!=i)
				h=h*(u-x[ip+j])/(x[ip+i]-x[ip+j]);
		}

        w=w+h;
    }
    
	return(w);
}

//////////////////////////////////////////////////////////////////////
// 二元全区间插值
//
// 参数：
// 1. int n - x方向上给定结点的点数
// 2. double x[] - 一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)
// 3. int m - y方向上给定结点的点数
// 4. double y[] - 一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)
// 5. double z[] - 一维数组，长度为n x m，存放给定的n x m个结点的函数值z(i,j)，
//                 z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
// 6. double u - 存放插值点x坐标
// 7. double v - 存放插值点y坐标
//
// 返回值：double 型，指定函数值f(u, v)
//////////////////////////////////////////////////////////////////////
double CInterpolate::GetValueLagrange2(int n, double x[], int m, double y[], double z[], double u, double v)
{ 
	int ip,ipp,i,j,l,iq,iqq,k;
    double h,w,b[10];
    
	// 特例
	if (u<=x[0]) 
	{ 
		ip=1; 
		ipp=4;
	}
    else if (u>=x[n-1]) 
	{ 
		ip=n-3; 
		ipp=n;
	}
    else
    { 
		i=1; 
		j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
        { 
			l=(i+j)/2;
            if (u<x[l-1]) 
				j=l;
            else 
				i=l;
        }
        
		ip=i-3; 
		ipp=i+4;
    }
    
	if (ip<1) 
		ip=1;

    if (ipp>n) 
		ipp=n;

    if (v<=y[0]) 
	{ 
		iq=1; 
		iqq=4;
	}
    else if (v>=y[m-1]) 
	{ 
		iq=m-3; 
		iqq=m;
	}
    else
    { 
		i=1; 
		j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
        { 
			l=(i+j)/2;
            if (v<y[l-1]) 
				j=l;
            else 
				i=l;
        }
        
		iq=i-3; 
		iqq=i+4;
    }
    
	if (iq<1) 
		iq=1;

    if (iqq>m) 
		iqq=m;

    for (i=ip-1;i<=ipp-1;i++)
    { 
		b[i-ip+1]=0.0;
        for (j=iq-1;j<=iqq-1;j++)
        { 
			h=z[m*i+j];
            for (k=iq-1;k<=iqq-1;k++)
            {
				if (k!=j) 
					h=h*(v-y[k])/(y[j]-y[k]);
			}

            b[i-ip+1]=b[i-ip+1]+h;
        }
    }
    
	w=0.0;
    for (i=ip-1;i<=ipp-1;i++)
    { 
		h=b[i-ip+1];
        for (j=ip-1;j<=ipp-1;j++)
        {
			if (j!=i) 
				h=h*(u-x[j])/(x[i]-x[j]);
		}

        w=w+h;
    }
    
	return(w);
}