www.pudn.com > test.rar > test.c

 
#include "stdafx.h" 
#include  
#include  
#include  
#include "InsertValue.h"  
 
//第二种边界条件的三次样条函数插值、微商与积分 
double espl2(double x[], double y[], int n, double dy[], double ddy[], double t[], int m, double z[], double dz[], double ddz[]) 
{ 
 int i,j; 
 double h0,h1,alpha,beta,g,*s; 
 s= (double *)malloc(n*sizeof(double)); 
 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]; 
 } 
 free(s); 
 return(g); 
}  
 
//光滑不等距插值 
void enspl(double x[], double y[], int n, int k, double t, double s[5]) 
{ 
 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; 
 if (n==1) 
 { 
  s[0]=y[0]; 
  s[4]=y[0]; 
  return; 
 } 
 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; 
 } 
 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=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; 
}  
 
//埃尔米特不等距插值 
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); 
}  
 
//一元全区间等距插值 
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= float((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[i]; 
 } 
 return(z); 
}  
 
//第一种边界条件的三次样条函数插值、微商与积分 
double espl1(double x[], double y[], int n, double dy[], double ddy[], double t[], int m, double z[], double dz[], double ddz[]) 
{ 
 int i,j; 
 double h0,h1,alpha,beta,g,*s; 
 s= (double *)malloc(n*sizeof(double)); 
 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]; 
 } 
 free(s); 
 return(g); 
}  
 
//连分式不等距插值 
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 (tx[n-5]) 
 { 
  k=n-8; 
  m=8; 
 } 
 else 
 { 
  k=1; 
  j=n; 
  while (j-k!=1) 
  { 
   i=(k+j)/2; 
   if (t=1;i--) 
  z=b[i-1]+(t-x[i+k-1])/z; 
 return(z); 
}  
 
//埃待金不等距逐步插值 
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 (tfabs(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); 
}  
 
//一元三点等距插值 
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[i]; 
 } 
 return(z); 
}  
 
//二元三点插值 
double eslq3(double x[], double y[], double z[], int n, int m, 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=y[m-2]) 
  iq=m-3; 
 else 
 { 
  i=1; 
  j=m; 
  while (((i-j)!=1)&&((i-j)!=-1)) 
  { 
   l=(i+j)/2; 
   if (vn-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); 
}  
 
//光滑等距插值 
void eespl(double x0, double h, int n, double y[], int k, double t, double s[5]) 
{ 
 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; 
 if (n==1) 
 { 
  s[0]=y[0]; 
  s[4]=y[0]; 
  return; 
 } 
 if (n==2) 
 { 
  s[0]=y[0]; 
  s[1]=(y[1]-y[0])/h; 
  if (k<0) 
   s[4]=(y[1]*(t-x0)-y[0]*(t-x0-h))/h; 
  return; 
 } 
 if (k<0) 
 { 
  if (t<=x0+h) 
   kk=0; 
  else if (t>=x0+(n-1)*h) 
   kk=n-2; 
  else 
  { 
   kk=1; 
   m=n; 
   while (((kk-m)!=1)&&((kk-m)!=-1)) 
   { 
    l=(kk+m)/2; 
    if (t=n-1) 
  kk=n-2; 
 u[2]=(y[kk+1]-y[kk])/h; 
 if (n==3) 
 { 
  if (kk==0) 
  { 
   u[3]=(y[2]-y[1])/h; 
   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])/h; 
   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])/h; 
   if (kk==1) 
   { 
    u[1]=(y[1]-y[0])/h; 
    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])/h; 
   } 
   else 
   { 
    u[1]=2.0*u[2]-u[3]; 
    u[0]=2.0*u[1]-u[2]; 
    u[4]=(y[3]-y[2])/h; 
   } 
  } 
  else if (kk>=(n-3)) 
  { 
   u[1]=(y[kk]-y[kk-1])/h; 
   if (kk==(n-3)) 
   { 
    u[3]=(y[n-1]-y[n-2])/h; 
    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])/h; 
   } 
   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])/h; 
   } 
  } 
  else 
  { 
   u[1]=(y[kk]-y[kk-1])/h; 
   u[0]=(y[kk-1]-y[kk-2])/h; 
   u[3]=(y[kk+2]-y[kk+1])/h; 
   u[4]=(y[kk+3]-y[kk+2])/h; 
  } 
 } 
 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]=h; 
 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*h); 
  s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p; 
 } 
 return; 
}  
 
//埃尔米特等距插值 
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); 
}  
 
//第三种边界条件的三次样条函数插值、微商与积分 
double espl3(double x[], double y[], int n, double dy[], double ddy[], double t[], int m, double z[], double dz[], double ddz[]) 
{ 
 int i,j; 
 double h0,y0,h1,y1,alpha,beta,u,g,*s; 
 s = (double *)malloc(n*sizeof(double)); 
 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=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 (h0x[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]; 
 } 
 free(s); 
 return(g); 
}  
 
//二元全区间插值 
double eslgq(double x[], double y[], double z[], int n, int m, 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 (un) 
  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 (vm) 
  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); 
}  
 
//一元三点不等距插值 
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(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; 
  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)*h; 
    xj=x0+(j+k-1)*h; 
    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)*h))/z; 
 return(z); 
}  
 
//埃特金等距逐步插值 
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 (tfabs(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[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); 
}