/****************函数说明*********************/
//px py为已知的数据点，xs为要插值的x坐标，最终会得到xs坐标下的y值
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace spline
{
    class Program
    {
writeline函数
        static void Main(string[] args)
        {
            point[] points = new point[13];
            double[] px = { 64, 304, 544, 1035, 1502, 2061, 2540, 2939, 3498, 3897, 4456, 4696, 4936 };
            double []py={4.663,4.476,4.969,4.853,4.766,5.2415,4.795,4.7055,5.4565,5.5515,5.411,5.8385,6.7085};
            for (int i = 0; i < 13; i++)
            {
                points[i] = new point();
                points[i].x =px[i];
                points[i].y = py[i];
            }
            point.DeSortX(points);
            double[] xs = {64, 304, 544, 1023, 1502, 1981, 2460, 2939, 3418, 3897, 4376, 4616, 4856};
            splineInsertPoint(points, xs, 1);
            Console.ReadLine();
        }
        static double[] splineInsertPoint(point[] points, double[] xs, int chf)
        {
          int plength = points.Length;
            double[] h = new double[plength];
            double[] f = new double[plength];
            double[] l = new double[plength];
            double[] v = new double[plength];
            double[] g = new double[plength];
            //三转角法的待定一阶系数法。基本方程为l(i)m(i-1)+2m(i)+u(i)m(i+1)=g(i),i=1,2,...,n-1
            //目前有n-1个方程，n+1个未知量（从m0到mn），需要添加两个边界条件
            for (int i = 0; i < plength - 1; i++)  //共有plength(n+1)个点，h总有plength-1(n)个值
            {
                h[i] = points[i + 1].x - points[i].x;
            }
            for (int i = 1; i < plength - 1; i++)
            {
                l[i] = h[i] / (h[i - 1] + h[i]);   
                v[i] = h[i - 1] / (h[i - 1] + h[i]);
                g[i] = 3 * (l[i] * f[i - 1] + v[i] * f[i]);  //i=1,2,3,...,n-1
            }
            //求解矩阵，用追赶法
            double[] U = new double[plength];  //上三角矩阵U ,plength=n+1阶矩阵
            double[] L = new double[plength];  //下三角矩阵L
            double[] m = new double[plength];  //用于保存m0-mn的一阶系数数值
            double[] Y = new double[plength];  //中间矩阵Y
            /***********************矩阵分界，求解待定的一阶系数，m1到m(n-1)*************************/
            U[0] = 2;
            v[0] = 1;
            for (int i = 1; i < plength;i++ )
            {
                L[i] = l[i] / U[i-1];
                U[i] = 2 - v[i - 1] * L[i];
            }
            Y[0] = g[0];
            for(int i = 1; i < plength;i++)
            {
                Y[i] = g[i] - L[i] * Y[i - 1];
            }
            m[plength - 1] = Y[plength - 1] / U[plength - 1];
            for(int i = plength - 2; i>=0; i--)
            {
                m[i] = (Y[i] - v[i] * m[i + 1]) / U[i];
            }
            /*****************************************************************************************/
            int xlength = xs.Length;  //记录需要插值的点的个数
            double[] insertRes = new double[xlength];  //用于返回插值点的Y坐标
            for (int i = 0; i < xlength; i++)
            {
                int j = 0;
                for (j = 0; j < plength; j++)
                {
                    if (xs[i] < points[j].x)
                        break;
                }
                j = j - 1;
                //Console.WriteLine(j);  //表示要插值的点数组的坐标在第j+1组区间
                if (j == -1 || j == plength - 1)
                {
                    if (j == -1)
                        throw new Exception("插值下边界超出");
                    if (j == plength - 1 && xs[i] == points[j].x) 
                        //这种情况是，需要插值的点正好等于一系列点x坐标最大的那个点
                        insertRes[i] = points[j].y;
                    else
                        throw new Exception("插值下边界超出");
                }
                else
                {
                    //根据每个分段的三次多项式求值
                    insertRes[i] = (h[j]+2*(xs[i]-points[j].x))*Math.Pow(xs[i]-points[j+1].x,2)*points[j].y/Math.Pow(h[j],3) + (h[j]+2*(points[j+1].x-xs[i]))*Math.Pow(xs[i]-points[j].x,2)*points[j+1].y/Math.Pow(h[j],3) + (xs[i]-points[j].x)*Math.Pow(xs[i]-points[j+1].x,2)*m[j]/Math.Pow(h[j],2) + (xs[i]-points[j+1].x)*Math.Pow(xs[i]-points[j].x,2)*m[j+1]/Math.Pow(h[j],2);

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。