数值计算程序大放送-插值算法

dormouse_buaa · 发表于 2005-1-19 23:01:13

<

>选自<<徐世良数值计算程序集(C)>>
<

>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。
<

>今天都给贴出来了
//一元全区间不等距插值，
//利用七次(两边为四次)拉格朗日插值多项式计算插值点t的函数值
//x-长度为n的数组，存放n个节点xi,要求增序
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值x
double enlgr(double x[],double y[],int n,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<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-x[j]);
}
z=z+s*y;
 }
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:01:53

//一元全区间等距插值
//利用七次(两边为四次)拉格朗日插值多项式计算插值点t的函数值
//x0-第一个节点x0
//h-等距节点的步长
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值 x
double eelgr(double x0,double h,int n,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-h))/h;
 return(z);
}
if (t>x0)
{
 p=(t-x0)/h;
 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*h;
 for (j=k; j<=m; j++)
 {
if (j!=i)
{
xj=x0+j*h;
s=s*(t-xj)/(xi-xj);
}
 }
 z=z+s*y;
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:02:08

//一元三点不等距插值
//利用抛物线插值公式计算指定插值点t的函数值
//x-长度为n的数组，存放n个不等距节点xi,要求增序
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值x
double enlg3(double x[],double y[],int n,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-x[j]);
}
 }
 z=z+s*y;
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:02:25

//一元三点等距插值
//利用抛物线插值公式计算指定插值点t的函数值
//x0-第一个节点x0
//h-等距节点的步长
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值 x
double eelg3(double x0,double h,int n,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-h))/h;
 return(z);
}
if (t<=x0+h)
{
 k=0;
 m=2;
}
else if (t>=x0+(n-3)*h)
{
 k=n-3;
 m=n-1;
}
else
{
 i=(int)((t-x0)/h)+1;
 if (fabs(t-x0-i*h)>=fabs(t-x0-(i-1)*h))
 {
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*h;
 for (j=k;j<=m;j++)
 {
if (j!=i)
{
xj=x0+j*h;
s=s*(t-xj)/(xi-xj);
}
 }
 z=z+s*y;
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:02:39

//连分式不等距插值
//x-长度为n的数组，存放n个节点xi,要求增序
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值x
double enpqs(double x[],double y[],int n,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;
 for (j=1;j<=i-1;j++)
if (l==0)
{
if (fabs(h-b[j-1])+1.0==1.0)
{
 l=l+1;
}
else
{
 h=(x[i+k-1]-x[j+k-1])/(h-b[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);
}

dormouse_buaa · 发表于 2005-1-19 23:02:54

//连分式等距插值
//x0-第一个节点x0
//h-等距节点的步长
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值 x
double eepqs(double x0,double h,int n,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*h))
{
 k=0;
 m=8;
}
else if (t>(x0+(n-5)*h))
{
 k=n-8;
 m=8;
}
else
{
 k=(int)((t-x0)/h)-3;
 m=8;
}
b[0]=y[k];
for (i=2;i<=m;i++)
{
 hh=y[i+k-1];
 l=0;
 for (j=1;j<=i-1;j++)
if (l==0)
 {
if (fabs(hh-b[j-1])+1.0==1.0)
{
 l=l+1;
}
else
 {
 xi=x0+(i+k-1)*h;
 xj=x0+(j+k-1)*h;
 hh=(xi-xj)/(hh-b[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)*h))/z;
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:03:19

//埃尔米特（Hermite）不等距插值
//x-长度为n的数组，存放n个节点xi,要求增序
//y-长度为n的数组，存放n个节点yi
//dy-长度为n的数组，存放n个节点上的一阶导数值
//n-给定节点的个数
//t-存放指定插值点的值x
double enhmt(double x[],double y[],double dy[],int n,double t)
{
int i,j;
double z,p,q,s;
z=0.0;
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);
}

dormouse_buaa · 发表于 2005-1-19 23:03:35

//埃尔米特（Hermite）等距插值
//x0-第一个节点x0
//h-等距节点的步长
//y-长度为n的数组，存放n个节点yi
//dy-长度为n的数组，存放n个节点上的一阶导数值
//n-给定节点的个数
//t-存放指定插值点的值 x
double eehmt(double x0,double h,int n,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)*h;
 for (j=1;j<=n;j++)
 {
p=x0+(j-1)*h;
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)*h));
}
 }
 q=y[i-1]+(t-q)*(dy[i-1]-2.0*y[i-1]*p);
 z=z+q*s;
}
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:03:55

//埃特金(Aitken)不等距逐步插值
//x-长度为n的数组，存放n个节点xi,要求增序
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值x
//eps-指定的精度要求
double enatk(double x[],double y[],int n,double t,double eps)
{
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;
 for (j=0;j<=i-1;j++)
 {
z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx);
 }
 yy=z;
}
while ((i!=m-1)&&(fabs(yy-yy[i-1])>eps));
return(z);
}

dormouse_buaa · 发表于 2005-1-19 23:04:11

//埃特金(Aitken)等距逐步插值
//x0-第一个节点x0
//h-等距节点的步长
//y-长度为n的数组，存放n个节点yi
//n-给定节点的个数
//t-存放指定插值点的值x
//eps-指定的精度要求
double eeatk(double x0,double h,int n,double y[],double t,double eps)
{
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)*h)
{
 k=n;
}
else
{
 k=1;
 j=n;
 while ((k-j!=1)&&(k-j!=-1))
 {
l=(k+j)/2;
 if (t<x0+(l-1)*h)
{
j=l;
}
 else
{
k=l;
}
 }
 if (fabs(t-(x0+(l-1)*h))>fabs(t-(x0+(j-1)*h)))
 {
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)*h;
 yy[i-1]=y[k-1];
 l=l+1; j=-j;
}
i=0;
do
{
 i=i+1;
 z=yy;
 for (j=0;j<=i-1;j++)
 {
z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx);
 }
 yy=z;
}
while ((i!=m-1)&&(fabs(yy-yy[i-1])>eps));
return(z);
}

		自动登录	找回密码
密码			注-册-帐-号