?? 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]);
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -