C#，码海拾贝（04）——拉格朗日（Lagrange）曲线插值算法及其拓展，《C#数值计算算法编程》源代码升级改进版

摘要：本文给出《C#数值计算算法编程》第六章“插值类设计”中“一元全区间不等距插值”和“一元全区间等距插值”，也就是拉格朗日插值算法的源代码，做了比较大的改动。

一、图书介绍

1.1 关于《C#数值计算算法编程》

《C#数值计算——算法编程》是2007年1月1日电子工业出版社出版的图书，由周长发编写。

ISBN:9787121032035

本书囊括了近90个实用经典算法，每一个算法都独立成节。每一节都包括算法原理、算法实现和示例3个部分。算法原理部分分别讨论每一种算法的计算原理；算法实现部分讨论用C#实现算法的技巧，针对不同的算法设计了6个算法类，给出完整的类和算法函数的源程序；示例部分介绍算法类的调用方式，给出调用算法的示例源程序、验证算法的示例数据和运算结果。本书内容丰富，讲解细致，实例丰富，适合涉及科学与工程数值计算工作的科研人员、工程技术人员、管理人员，以及大专院校相关专业的师生参考阅读。

作者虽然出了不少书，写了不少各种语言的代码，但都是实验室产品，无法适用于商业化、工业化场景，要做一些修正才基本能用。

1.2 插值与拟合

插值(Interpolation)是曲线必须经过所有给定数据点。

拟合(Fit)是指定幂次的曲线尽可能符合数据点的趋势（最佳），一般不过给定数据点。

1.3 第六章 插值算法

6.1 插值类设计

6.2 一元全区间不等距插值

6.3 一元全区间等距插值

6.4 一元三点不等距插值

6.5 一元三点等距插值

6.6 连分式不等距插值

6.7 连分式等距插值

6.8 埃尔米特不等距插值

6.9 埃尔米特等距插值

6.10 埃特金不等距逐步插值

6.11 埃特金等距逐步插值

6.12 光滑不等距插值

6.13 光滑等距插佔

6.14 第一种边界条件的三次样条函数插值、微商与积分

6.15 第二种边界条件的三次样条函数插值、微商与积分

6.16 第三种边界条件的三次样条函数插值、微商与积分

6.17 二元三点插们

6.18 二元全区间插值

二、代码升级

2.1 升级的理由

此类代码都是 FORTRAN -> C -> C++ -> C# 等一路改编过来的。

从变量定义到代码逻辑都有很深的 FORTRAN 遗迹。

2.2 升级的代码

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>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        public static double Lagrange(double[] x, double[] y, double t)
        {
            // x,y点数
            int n = x.Length;
            double z = 0.0;

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

            // 开始插值
            int ik = 0;
            while ((x[ik] < t) && (ik < n))
            {
                ik = ik + 1;
            }

            int k = ik - 4;
            if (k < 0)
            {
                k = 0;
            }
            int m = ik + 3;
            if (m > n - 1)
            {
                m = n - 1;
            }
            for (int i = k; i <= m; i++)
            {
                double s = 1.0;
                for (int j = k; j <= m; j++)
                {
                    if (j != i)
                    {
                        // 拉格朗日插值公式
                        s = s * (t - x[j]) / (x[i] - x[j]);
                    }
                }

                z = z + s * y[i];
            }

            return (z);
        }

        /// <summary>
        /// 一元全区间不等距插值
        /// </summary>
        /// <param name="points">点集（含XY坐标）</param>
        /// <param name="t"></param>
        /// <returns></returns>
        public static double Lagrange(PointF[] points, double t)
        {
            double[] x = new double[points.Length];
            double[] y = new double[points.Length];
            for (int i = 0; i < points.Length; i++)
            {
                x[i] = points[i].X;
                y[i] = points[i].Y;
            }
            return Lagrange(x, y, t);
        }

        /// <summary>
        /// 一元全区间不等距插值
        /// </summary>
        /// <param name="points">二元组类型的点集（含XY坐标）</param>
        /// <param name="t"></param>
        /// <returns></returns>
        public static double Lagrange(List<Tuple<double, double>> points, double t)
        {
            double[] x = new double[points.Count];
            double[] y = new double[points.Count];
            for (int i = 0; i < points.Count; i++)
            {
                x[i] = points[i].Item1;
                y[i] = points[i].Item2;
            }
            return Lagrange(x, y, t);
        }

        /// <summary>
        /// 一元全区间不等距插值，获得插值后的曲线（折线拟合）数据
        /// </summary>
        /// <param name="points">点集（含XY坐标）</param>
        /// <param name="segment_count">每数据段的分割数</param>
        /// <returns></returns>
        public static PointF[] Lagrange_Curve(PointF[] points, int segment_count = 10)
        {
            int n = points.Length;
            PointF[] segments = new PointF[n * segment_count + 1];
            for (int i = 0; i < points.Length - 1; i++)
            {
                double dt = (points[i + 1].X - points[i].X) / segment_count;
                double t = points[i].X;
                for (int j = 0; j <= segment_count; j++, t += dt)
                {
                    PointF p = new PointF(0.0F, 0.0F);
                    p.X = (float)t;
                    if (j == 0) p.Y = points[i].Y;
                    else if (j == segment_count) p.Y = points[i + 1].Y;
                    else p.Y = (float)(Lagrange(points, t));
                    segments[i] = p;
                }
            }
            return segments;
        }

        /// <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>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        public static double Lagrange(double x0, double step, double[] y, double t)
        {
            double[] x = new double[y.Length];
            for (int i = 0; i < y.Length; i++, x0 += step)
            {
                x[i] = x0;
            }
            return Lagrange(x, y, t);
        }

#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>
        /// <returns>指定的查指点t的函数近似值y=f(t)</returns>
        public static double Lagrange_Equidistant(double x0, double step, double[] y, double t)
        {
            int n = y.Length;
            double z = 0.0;

            // 特例处理
            if (n < 1)
            {
                return (z);
            }
            else if (n == 1)
            {
                z = y[0];
                return (z);
            }
            else if (n == 2)
            {
                z = (y[1] * (t - x0) - y[0] * (t - x0 - step)) / step;
                return (z);
            }

            // 开始插值
            int ik = 0;
            if (t > x0)
            {
                double p = (t - x0) / step;
                ik = (int)p;
                double q = (float)ik;
                if (p > q)
                {
                    ik = ik + 1;
                }
            }
            else
            {
                ik = 0;
            }
            int k = ik - 4;
            if (k < 0)
            {
                k = 0;
            }
            int m = ik + 3;
            if (m > n - 1)
            {
                m = n - 1;
            }
            for (int i = k; i <= m; i++)
            {
                double s = 1.0;
                double xi = x0 + i * step;

                for (int j = k; j <= m; j++)
                {
                    if (j != i)
                    {
                        double xj = x0 + j * step;
                        // 拉格朗日插值公式
                        s = s * (t - xj) / (xi - xj);
                    }
                }

                z = z + s * y[i];
            }

            return (z);
        }
#endif
    }
}

POWER BY 315SOFT.COM