?? interpolation.cs
字號:
/*
* 進行插值的類Interpolation
* 周長發編制
*/
using System;
namespace 土地適宜評價系統.Algorithm
{
/**
* 進行插值的類Interpolation
* @author 周長發
* @version 1.0
*/
public class Interpolation
{
/**
* 一元全區間不等距插值
*
* @param n - 結點的個數
* @param x - 一維數組，長度為n，存放給定的n個結點的值x(i)，
* 要求x(0)<x(1)<...<x(n-1)
* @param y - 一維數組，長度為n，存放給定的n個結點的函數值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double 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);
}
/**
* 一元全區間等距插值
*
* @param n - 結點的個數
* @param x0 - 存放等距n個結點中第一個結點的值
* @param xStep - 等距結點的步長
* @param y - 一維數組，長度為n，存放給定的n個結點的函數值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double GetValueLagrange(int n, double x0, double xStep, double[] y, double t)
{
int i,j,k,m;
double z,s,xi,xj;
double 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);
}
/**
* 一元三點不等距插值
*
* @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 t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double 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 (Math.Abs(t-x[k])<Math.Abs(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);
}
/**
* 一元三點等距插值
*
* @param n - 結點的個數
* @param x0 - 存放等距n個結點中第一個結點的值
* @param xStep - 等距結點的步長
* @param y - 一維數組，長度為n，存放給定的n個結點的函數值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double 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 (Math.Abs(t-x0-i*xStep)>=Math.Abs(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);
}
/**
* 連分式不等距插值
*
* @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 t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double GetValuePqs(int n, double[] x, double[] y, double t)
{
int i,j,k,m,l;
double z,h;
double[] b = new double[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 (Math.Abs(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);
}
/**
* 連分式等距插值
*
* @param n - 結點的個數
* @param x0 - 存放等距n個結點中第一個結點的值
* @param xStep - 等距結點的步長
* @param y - 一維數組，長度為n，存放給定的n個結點的函數值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double GetValuePqs(int n, double x0, double xStep, double[] y, double t)
{
int i,j,k,m,l;
double z,hh,xi,xj;
double[] b = new double[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 (Math.Abs(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);
}
/**
* 埃爾米特不等距插值
*
* @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，存放給定的n個結點的函數導數值y'(i)，
* y'(i) = f'(x(i)), i=0,1,...,n-1
* @param t - 存放指定的插值點的值
* @return double 型，指定的查指點t的函數近似值f(t)
*/
public static double GetValueHermite(int n, double[] x, double[] y, double[] dy, 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);
}
/**
* 埃爾米特等距插值
*
* @param n - 結點的個數
* @param x0 - 存放等距n個結點中第一個結點的值
* @param xStep - 等距結點的步長
* @param y - 一維數組，長度為n，存放給定的n個結點的函數值y(i)，
* y(i) = f(x(i)), i=0,1,...,n-1
* @param dy - 一維數組，長度為n，存放給定的n個結點的函數導數值y'(i)，
* y'(i) = f'(x(i)), i=0,1,...,n-1
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -