spline.cs

C#编写的三次样条插值函数

using System;
class MathMethod
{
	/// <summary>
	/// TSS法解三对角线方程组
	/// </summary>
	/// <param name="a">a[n]</param>
	/// <param name="b">b[n]</param>
	/// <param name="c">c[n]</param>
	/// <param name="d">d[n]</param>
	/// <param name="n">n</param>
	static void TSS(double[] a,double[] b,double[] c,double[] d,int n)
	{
		d[0] /= b[0];
		int i = 0;
		for(i=1;i<n;i++)
		{
			c[i-1] /= b[i-1];
			b[i] -= a[i]*c[i-1];
			d[i] = (d[i]-a[i]*d[i-1])/b[i];
		}
		for(i=n-2;i>=0;i--)                                                                                                                                                                              
			d[i] -= c[i]*d[i+1];
	}
	/// <summary>
	/// 三次样条曲线插值
	/// </summary>
	/// <param name="x">x序列</param>
	/// <param name="y">y序列</param>
	/// <param name="M">用于存放插值系数，从C++转过来的，懒得改了，保留</param>
	/// <param name="n">序列长度</param>
	/// <param name="x0">样条插值的端点处理，详细见数值计算数</param>
	/// <param name="d0"></param>
	/// <param name="un"></param>
	/// <param name="dn"></param>
	static void SplineM(double[] x,double[] y,double[] M,int n,double x0,double d0,double un,double dn)
	{
		int i = 0;
		for(i = 0;i<=n;++i)
			M[i] = y[i];
		for(int k=1;k<=2;++k)
		{
			for(i = n;i>=k;--i)
			{
				M[i] = (M[i] -M[i-1])/(x[i] -x[i-k]);
			}
		}
		double[] h = new double[n+1];
		double[] a = new double[n+1];
		double[] b = new double[n+1];
		double[] c = new double[n+1];
		h[1] = x[1]-x[0];
		for(i = 1;i<n;++i)
		{
			h[i+1] = x[i+1]-x[i];
			c[i]=h[i+1]/(h[i]+h[i+1]);
			a[i]=1-c[i];
			b[i]=2;
			M[i]=6*M[i+1];
		}
		M[0]=d0;
		M[n]=dn;
		c[0]=x0;
		b[0]=2;
		a[n]=un;
		b[n]=2;
		TSS(a,b,c,M,n);
	}
	/// <summary>
	/// 计算插值
	/// </summary>
	/// <param name="xv">x序列</param>
	/// <param name="y">y序列</param>
	/// <param name="M">插值系数序列</param>
	/// <param name="n">序列长度</param>
	/// <param name="xIn">插值点坐标</param>
	/// <returns>插值点结果</returns>
	static double Evaspline(double[] xv,double[] y,double[] M,int n,double xIn)
	{
		double r =0.0;
		int i;
		double difa,difb,h;
		i = -1;
		do
		{
			i++;
		}while(!(((xIn>=xv[i])&&(xIn<xv[i+1]))||(i==n-1)));
		difa = xIn-xv[i];
		difb = xv[i+1]-xIn;
		h = xv[i+1]-xv[i];
		r = M[i]*Math.Pow(difb,3)/6+M[i+1]*Math.Pow(difa,3)/6+(y[i]-M[i]*h*h/6)*difb+(y[i+1]-M[i+1]*h*h/6)*difa;
		r /= h;
		return r;
	}
}