📄 interpolation.cs
字号:
{
p=t-(x0+kk*xStep);
s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
}
return s[4];
}
/**
* 第一种边界条件的三次样条函数插值、微商与积分
*
* @param n - 结点的个数
* @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)
* @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param dy - 一维数组，长度为n，调用时，dy(0)存放给定区间的左端点处的一阶导数值，
* dy(n-1)存放给定区间的右端点处的一阶导数值。返回时，存放n个给定点处的
* 一阶导数值y'(i)，i=0,1,...,n-1
* @param ddy - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
* i=0,1,...,n-1
* @param m - 指定插值点的个数
* @param t - 一维数组，长度为m，存放m个指定的插值点的值。
* @param z - 一维数组，长度为m，存放m个指定的插值点处的函数值
* @param dz - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
* @param ddz - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
* @return double 型，指定函数的x(0)到x(n-1)的定积分值
*/
public static double 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;
// 初值
double[] 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.0*(y[n-2]-y[n-1])/h1+2.0*(2.0*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];
}
return(g);
}
/**
* 第二种边界条件的三次样条函数插值、微商与积分
*
* @param n - 结点的个数
* @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)
* @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param dy - 一维数组，长度为n，调用时，dy(0)存放给定区间的左端点处的一阶导数值，
* dy(n-1)存放给定区间的右端点处的一阶导数值。返回时，存放n个给定点处的
* 一阶导数值y'(i)，i=0,1,...,n-1
* @param ddy - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
* i=0,1,...,n-1
* @param m - 指定插值点的个数
* @param t - 一维数组，长度为m，存放m个指定的插值点的值。
* @param z - 一维数组，长度为m，存放m个指定的插值点处的函数值
* @param dz - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
* @param ddz - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
* @return double 型，指定函数的x(0)到x(n-1)的定积分值
*/
public static double 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=0,alpha,beta,g;
// 初值
double[] 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.0*(y[n-2]-y[n-1])/h1+2.0*(2.0*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];
}
return(g);
}
/**
* 第三种边界条件的三次样条函数插值、微商与积分
*
* @param n - 结点的个数
* @param x - 一维数组，长度为n，存放给定的n个结点的值x(i)
* @param y - 一维数组，长度为n，存放给定的n个结点的函数值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param dy - 一维数组，长度为n，调用时，dy(0)存放给定区间的左端点处的一阶导数值，
* dy(n-1)存放给定区间的右端点处的一阶导数值。返回时，存放n个给定点处的
* 一阶导数值y'(i)，i=0,1,...,n-1
* @param ddy - 一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，
* i=0,1,...,n-1
* @param m - 指定插值点的个数
* @param t - 一维数组，长度为m，存放m个指定的插值点的值。
* @param z - 一维数组，长度为m，存放m个指定的插值点处的函数值
* @param dz - 一维数组，长度为m，存放m个指定的插值点处的一阶导数值
* @param ddz - 一维数组，长度为m，存放m个指定的插值点处的二阶导数值
* @return double 型，指定函数的x(0)到x(n-1)的定积分值
*/
public static double 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=0,beta=0,u,g;
// 初值
double[] 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.0*(y[n-2]-y[n-1])/h1+2.0*(2.0*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];
}
return(g);
}
/**
* 二元三点插值
*
* @param n - x方向上给定结点的点数
* @param x - 一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)
* @param m - y方向上给定结点的点数
* @param y - 一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)
* @param z - 一维数组，长度为n x m，存放给定的n x m个结点的函数值z(i,j)，
* z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
* @param u - 存放插值点x坐标
* @param v - 存放插值点y坐标
* @return double 型，指定函数值f(u, v)
*/
public static double 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 = new double[3];
double 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 (Math.Abs(u-x[i-1])<Math.Abs(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 (Math.Abs(v-y[i-1])<Math.Abs(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);
}
/**
* 二元全区间插值
*
* @param n - x方向上给定结点的点数
* @param x - 一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)
* @param m - y方向上给定结点的点数
* @param y - 一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)
* @param z - 一维数组，长度为n x m，存放给定的n x m个结点的函数值z(i,j)，
* z(i,j) = f(x(i), y(j)), i=0,1,...,n-1, j=0,1,...,m-1
* @param u - 存放插值点x坐标
* @param v - 存放插值点y坐标
* @return double 型，指定函数值f(u, v)
*/
public static double 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;
double[] b = new double[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);
}
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -