?? interpolate.cpp
字號(hào):
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];
}
//////////////////////////////////////////////////////////////////////
// 第一種邊界條件的三次樣條函數(shù)插值、微商與積分
//
// 參數(shù)：
// 1. int n - 結(jié)點(diǎn)的個(gè)數(shù)
// 2. double x[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的值x(i)
// 3. double y[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的函數(shù)值y(i)，
// y(i) = f(x(i)), i=0,1,...,n-1
// 4. double dy[] - 一維數(shù)組，長(zhǎng)度為n，調(diào)用時(shí)，dy(0)存放給定區(qū)間的左端點(diǎn)處的一階導(dǎo)數(shù)值，
// dy(n-1)存放給定區(qū)間的右端點(diǎn)處的一階導(dǎo)數(shù)值。返回時(shí)，存放n個(gè)給定點(diǎn)處的
// 一階導(dǎo)數(shù)值y'(i)，i=0,1,...,n-1
// 5. double ddy[] - 一維數(shù)組，長(zhǎng)度為n，返回時(shí)，存放n個(gè)給定點(diǎn)處的二階導(dǎo)數(shù)值y''(i)，
// i=0,1,...,n-1
// 6. int m - 指定插值點(diǎn)的個(gè)數(shù)
// 7. double t[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)的值。
// 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的函數(shù)值
// 9. double dz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的一階導(dǎo)數(shù)值
// 10. double ddz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的二階導(dǎo)數(shù)值
//
// 返回值：double 型，指定函數(shù)的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);
}
//////////////////////////////////////////////////////////////////////
// 第二種邊界條件的三次樣條函數(shù)插值、微商與積分
//
// 參數(shù)：
// 1. int n - 結(jié)點(diǎn)的個(gè)數(shù)
// 2. double x[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的值x(i)
// 3. double y[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的函數(shù)值y(i)，
// y(i) = f(x(i)), i=0,1,...,n-1
// 4. double dy[] - 一維數(shù)組，長(zhǎng)度為n，返回時(shí)，存放n個(gè)給定點(diǎn)處的一階導(dǎo)數(shù)值y'(i)，
// i=0,1,...,n-1
// 5. double ddy[] - 一維數(shù)組，長(zhǎng)度為n，返回時(shí)，存放n個(gè)給定點(diǎn)處的二階導(dǎo)數(shù)值y''(i)，
// i=0,1,...,n-1，調(diào)用時(shí)，ddy(0)存放給定區(qū)間的左端點(diǎn)處的二階導(dǎo)數(shù)值，
// ddy(n－1)存放給定區(qū)間的右端點(diǎn)處的二階導(dǎo)數(shù)值
// 6. int m - 指定插值點(diǎn)的個(gè)數(shù)
// 7. double t[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)的值。
// 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的函數(shù)值
// 9. double dz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的一階導(dǎo)數(shù)值
// 10. double ddz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的二階導(dǎo)數(shù)值
//
// 返回值：double 型，指定函數(shù)的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);
}
//////////////////////////////////////////////////////////////////////
// 第三種邊界條件的三次樣條函數(shù)插值、微商與積分
//
// 參數(shù)：
// 1. int n - 結(jié)點(diǎn)的個(gè)數(shù)
// 2. double x[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的值x(i)
// 3. double y[] - 一維數(shù)組，長(zhǎng)度為n，存放給定的n個(gè)結(jié)點(diǎn)的函數(shù)值y(i)，
// y(i) = f(x(i)), i=0,1,...,n-1
// 4. double dy[] - 一維數(shù)組，長(zhǎng)度為n，返回時(shí)，存放n個(gè)給定點(diǎn)處的一階導(dǎo)數(shù)值y'(i)，
// i=0,1,...,n-1
// 5. double ddy[] - 一維數(shù)組，長(zhǎng)度為n，返回時(shí)，存放n個(gè)給定點(diǎn)處的二階導(dǎo)數(shù)值y''(i)，
// i=0,1,...,n-1
// 6. int m - 指定插值點(diǎn)的個(gè)數(shù)
// 7. double t[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)的值。
// 要求x(0)<t(j)<x(n-1), j=0,1,…,m-1
// 8. double z[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的函數(shù)值
// 9. double dz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的一階導(dǎo)數(shù)值
// 10. double ddz[] - 一維數(shù)組，長(zhǎng)度為m，存放m個(gè)指定的插值點(diǎn)處的二階導(dǎo)數(shù)值
//
// 返回值：double 型，指定函數(shù)的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);
}
//////////////////////////////////////////////////////////////////////
// 二元三點(diǎn)插值
//
// 參數(shù)：
// 1. int n - x方向上給定結(jié)點(diǎn)的點(diǎn)數(shù)
// 2. double x[] - 一維數(shù)組，長(zhǎng)度為n，存放給定n x m 個(gè)結(jié)點(diǎn)x方向上的n個(gè)值x(i)
// 3. int m - y方向上給定結(jié)點(diǎn)的點(diǎn)數(shù)
// 4. double y[] - 一維數(shù)組，長(zhǎng)度為m，存放給定n x m 個(gè)結(jié)點(diǎn)y方向上的m個(gè)值y(i)
// 5. double z[] - 一維數(shù)組，長(zhǎng)度為n x m，存放給定的n x m個(gè)結(jié)點(diǎn)的函數(shù)值z(mì)(i,j)，
// z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
// 6. double u - 存放插值點(diǎn)x坐標(biāo)
// 7. double v - 存放插值點(diǎn)y坐標(biāo)
//
// 返回值：double 型，指定函數(shù)值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);
}
//////////////////////////////////////////////////////////////////////
// 二元全區(qū)間插值
//
// 參數(shù)：
// 1. int n - x方向上給定結(jié)點(diǎn)的點(diǎn)數(shù)
// 2. double x[] - 一維數(shù)組，長(zhǎng)度為n，存放給定n x m 個(gè)結(jié)點(diǎn)x方向上的n個(gè)值x(i)
// 3. int m - y方向上給定結(jié)點(diǎn)的點(diǎn)數(shù)
// 4. double y[] - 一維數(shù)組，長(zhǎng)度為m，存放給定n x m 個(gè)結(jié)點(diǎn)y方向上的m個(gè)值y(i)
// 5. double z[] - 一維數(shù)組，長(zhǎng)度為n x m，存放給定的n x m個(gè)結(jié)點(diǎn)的函數(shù)值z(mì)(i,j)，
// z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
// 6. double u - 存放插值點(diǎn)x坐標(biāo)
// 7. double v - 存放插值點(diǎn)y坐標(biāo)
//
// 返回值：double 型，指定函數(shù)值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);
}
?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -