C#,码海拾贝(08)——埃特金(Aitken)逐步曲线插值算法,《C#数值计算算法编程》源代码升级改进版

埃特金逐步线性插值法(Aitken successive linear interpolation method)一种能逐步升阶的插值方法.用拉格朗日插值多项式计算函数近似值时,如果精度不满足要求,需增加插值节点以提高插值多项式次数时,原来算出的结果均不能利用,必须重新计算.用埃特金逐步线性插值法能克服这一缺点。这个算法适用于计算机上计算,且具有自动选取节卢并逐步比较精度的特点,程序也较简单一般以斜结上两个相邻插值多项式的值之差满足所需精度作大计算过程终止标志。

Aitken Interpolation -- from Wolfram MathWorld Applied Mathematics Numerical Methods Approximation Theory Interpolation Aitken Interpolation An algorithm similar to Neville's algorithm for constructing the Lagrange interpolating polynomial

艾特肯插值——来自Wolfram MathWorld应用数学数值方法近似理论插值艾特肯插补一种类似于Neville算法的拉格朗日插值多项式构造算法。

AITKEN'S INTERPOLATION
Aitken's procedure yields systematically and successively better interpolation polynomials corresponding to successively higher order truncation of Newton's divided difference formula.

艾肯插值

Aitken程序系统地、连续地产生更好的插值多项式,对应于牛顿分差公式的连续高阶截断。

using System;
using System.Drawing;
using System.Collections;
using System.Collections.Generic;

namespace Zhou.CSharp.Algorithm
{
    /// <summary>
    /// 插值计算类Interpolation.cs
    /// 作者:周长发
    /// 改编:深度混淆
    /// https://blog.csdn.net/beijinghorn
    /// </summary>
    public static partial class Interpolation
    {
        /// <summary>
        /// 埃特金不等距逐步插值
        /// </summary>
        /// <param name="x">一维数组,长度为n,存放给定的n个结点的值x(i),要求x(0)<x(1)<...<x(n-1)</param>
        /// <param name="y">一维数组,长度为n,存放给定的n个结点的函数值y(i),y(i) = f(x(i)), i=0,1,...,n-1</param>
        /// <param name="t">存放指定的插值点的x值</param>
        /// <param name="eps">控制精度参数</param>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        public static double Aitken(double[] x, double[] y, double t, double eps)
        {
            int n = x.Length;
            double z = 0.0;
            // 特例处理
            if (n < 1)
            {
                return (z);
            }
            if (n == 1)
            {
                z = y[0];
                return (z);
            }

            // 开始插值
            int m = 10;
            if (m > n)
            {
                m = n;
            }
            int k;
            if (t <= x[0])
            {
                k = 1;
            }
            else if (t >= x[n - 1])
            {
                k = n;
            }
            else
            {
                k = 1;
                int j = n;
                int l = 0;
                while ((k - j != 1) && (k - j != -1))
                {
                    l = (k + j) / 2;
                    if (t < x[l - 1])
                    {
                        j = l;
                    }
                    else
                    {
                        k = l;
                    }
                }
                if (Math.Abs(t - x[l - 1]) > Math.Abs(t - x[j - 1]))
                {
                    k = j;
                }
            }

            double[] xx = new double[10];
            double[] yy = new double[10];
            {
                int j = 1;
                int l = 0;
                for (int 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;
                }
            }

            {
                int i = 0;
                do
                {
                    i = i + 1;
                    z = yy[i];

                    for (int 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) && (Math.Abs(yy[i] - yy[i - 1]) > eps));
            }

            return (z);
        }

        /// <summary>
        /// 埃特金等距逐步插值
        /// (使用非等距插值的方法)
        /// </summary>
        /// <param name="x0">存放等距n个结点中第一个结点的值</param>
        /// <param name="step">等距结点的步长</param>
        /// <param name="y">一维数组,长度为n,存放给定的n个结点的函数值y(i),y(i) = f(x(i)), i=0,1,...,n-1</param>
        /// <param name="t">存放指定的插值点的x值</param>
        /// <param name="eps">精算精度</param>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        public static double Aitken(double x0, double step, double[] y, double t, double eps)
        {
            double[] x = new double[y.Length];
            for (int i = 0; i < y.Length; i++, x0 += step)
            {
                x[i] = x0;
            }
            return Aitken(x, y, t, eps);
        }

#if __OLD__
        /// <summary>
        /// 埃特金等距逐步插值
        /// </summary>
        /// <param name="x0">等距n个结点中第一个结点的值</param>
        /// <param name="step">等距结点的步长</param>
        /// <param name="y">一维数组,长度为n,存放给定的n个结点的函数值y(i),y(i) = f(x(i)), i=0,1,...,n-1</param>
        /// <param name="t">存放指定的插值点的x值</param>
        /// <param name="eps">控制精度参数</param>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        /// <returns></returns>
        public static double Aitken(double x0, double step, double[] y, double t, double eps)
        {
            int n = y.Length;
            double z = 0.0;

            // 特例处理
            if (n < 1)
            {
                return (z);
            }
            else if (n == 1)
            {
                z = y[0];
                return (z);
            }

            // 开始插值
            int m = 10;
            if (m > n)
            {
                m = n;
            }
            int k = 0;
            if (t <= x0)
            {
                k = 1;
            }
            else if (t >= x0 + (n - 1) * step)
            {
                k = n;
            }
            else
            {
                k = 1;
                int j = n;
                int w = 0;
                while ((k - j != 1) && (k - j != -1))
                {
                    w = (k + j) / 2;
                    if (t < x0 + (w - 1) * step)
                    {
                        j = w;
                    }
                    else
                    {
                        k = w;
                    }
                }
                if (Math.Abs(t - (x0 + (w - 1) * step)) > Math.Abs(t - (x0 + (j - 1) * step)))
                {
                    k = j;
                }
            }

            double[] xx = new double[10];
            double[] yy = new double[10];
            {
                int j = 1;
                int w = 0;
                for (int i = 1; i <= m; i++)
                {
                    k = k + j * w;
                    if ((k < 1) || (k > n))
                    {
                        w = w + 1;
                        j = -j;
                        k = k + j * w;
                    }

                    xx[i - 1] = x0 + (k - 1) * step;
                    yy[i - 1] = y[k - 1];
                    w = w + 1;
                    j = -j;
                }
            }
            {
                int i = 0;
                do
                {
                    i = i + 1;
                    z = yy[i];
                    for (int 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) && (Math.Abs(yy[i] - yy[i - 1]) > eps));
            }
            return (z);
        }
#endif
    }
}

POWER BY 315SOFT.COM

基于坐标点的计算,从点集计算插值曲线等拓展方法请参阅《拉格朗日插值算法及其拓展

猜你喜欢

转载自blog.csdn.net/beijinghorn/article/details/129829882