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