C#常用数学插值法

目录

1、分段线性插值

2、三次样条插值

3、拉格朗日插值

（1）一元全区间不等距插值

（2）一元全区间等距插值

4、埃尔米特插值

（1）埃尔米特不等距插值

（2）埃尔米特等距插值

---

1、分段线性插值

        /// <summary>  /// 分段线性插值，将一组数插值为所需点数  /// </summary>  /// <param name="dataIn">待插值的数据数组</param>  /// <param name="n">插值点数</param>  /// <returns>插值后的数据数组</returns>  public static float[] Interpolation(float[] dataIn, int n){float[] dataOut = new float[n];int lenIn = dataIn.Length;float[] a = new float[lenIn];float[] divIn = new float[lenIn];float[] divOut = new float[n];divIn[0] = 0;for (int i = 1; i < lenIn; i++){divIn[i] = divIn[i - 1] + 1;}divOut[0] = 0;for (int i = 1; i < n; i++){divOut[i] = divOut[i - 1] + lenIn / (float)n;}int k = 0;for (int i = k; i < n; i++){for (int j = 0; j < lenIn - 1; j++){if (divOut[i] >= divIn[j] && divOut[i] < divIn[j + 1]){dataOut[i] = (dataIn[j + 1] - dataIn[j]) * (divOut[i] - divIn[j]) / (divIn[j + 1] - divIn[j]) + dataIn[j];k = i;}}}return dataOut;}

2、三次样条插值

三次样条插值 C#代码实现_c# 三次样条插值_Big_潘大师的博客-CSDN博客

        /// <summary>/// 三次样条插值/// </summary>/// <param name="points">排序好的x、y点集合</param>/// <param name="xs">输入x轴数据，插值计算出对应的y轴点</param>/// <param name="chf">写1</param>/// <returns>返回计算好的Y轴数值</returns>public static double[] SplineInsertPoint(PointClass[] 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];for (int i = 0; i < plength - 1; i++){h[i] = points[i + 1].x - points[i].x;f[i] = (points[i + 1].y - points[i].y) / h[i];}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]);}double[] b = new double[plength];double[] tem = new double[plength];double[] m = new double[plength];double f0 = (points[0].y - points[1].y) / (points[0].x - points[1].x);double fn = (points[plength - 1].y - points[plength - 2].y) / (points[plength - 1].x - points[plength - 2].x);b[1] = v[1] / 2;for (int i = 2; i < plength - 2; i++){// Console.Write(" " + i);b[i] = v[i] / (2 - b[i - 1] * l[i]);}tem[1] = g[1] / 2;for (int i = 2; i < plength - 1; i++){//Console.Write(" " + i);tem[i] = (g[i] - l[i] * tem[i - 1]) / (2 - l[i] * b[i - 1]);}m[plength - 2] = tem[plength - 2];for (int i = plength - 3; i > 0; i--){//Console.Write(" " + i);m[i] = tem[i] - b[i] * m[i + 1];}m[0] = 3 * f[0] / 2.0;m[plength - 1] = fn;int xlength = xs.Length;double[] insertRes = new double[xlength];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);if (j == -1 || j == points.Length - 1){if (j == -1)throw new Exception("插值下边界超出");if (j == points.Length - 1 && xs[i] == points[j].x)insertRes[i] = points[j].y;elsethrow new Exception("插值下边界超出");}else{double p1;p1 = (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x);p1 = p1 * p1;double p2; p2 = (xs[i] - points[j].x) / (points[j + 1].x - points[j].x);p2 = p2 * p2;double p3; p3 = p1 * (1 + 2 * (xs[i] - points[j].x) / (points[j + 1].x - points[j].x)) * points[j].y + p2 * (1 + 2 * (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x)) * points[j + 1].y;double p4; p4 = p1 * (xs[i] - points[j].x) * m[j] + p2 * (xs[i] - points[j + 1].x) * m[j + 1];//         Console.WriteLine(m[j] + " " + m[j + 1] + " " + j);p4 = p4 + p3;insertRes[i] = p4;//Console.WriteLine("f(" + xs[i] + ")= " + p4);}}//Console.ReadLine();return insertRes;}

排序计算

     public class PointClass{public double x = 0;public double y = 0;public PointClass(){x = 0; y = 0;}//-------写一个排序函数，使得输入的点按顺序排列，是因为插值算法的要求是，x轴递增有序的---------public static PointClass[] DeSortX(PointClass[] points){int length = points.Length;double temx, temy;for (int i = 0; i < length - 1; i++){for (int j = 0; j < length - i - 1; j++)if (points[j].x > points[j + 1].x){temx = points[j + 1].x;points[j + 1].x = points[j].x;points[j].x = temx;temy = points[j + 1].y;points[j + 1].y = points[j].y;points[j].y = temy;}}return points;}}

3、拉格朗日插值

（1）一元全区间不等距插值

        /// <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);}

（2）一元全区间等距插值

        /// <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);}

4、埃尔米特插值

（1）埃尔米特不等距插值

        /// <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="dy">一维数组，长度为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 Hermite(double[] x, double[] y, double[] dy, double t){int n = x.Length;double z = 0.0;// 循环插值for (int i = 1; i <= n; i++){double s = 1.0;for (int j = 1; j <= n; j++){if (j != i){s = s * (t - x[j - 1]) / (x[i - 1] - x[j - 1]);}}s = s * s;double p = 0.0;for (int j = 1; j <= n; j++){if (j != i){p = p + 1.0 / (x[i - 1] - x[j - 1]);}}double q = y[i - 1] + (t - x[i - 1]) * (dy[i - 1] - 2.0 * y[i - 1] * p);z = z + q * s;}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="dy">一维数组，长度为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 Hermite(double x0, double step, double[] y, double[] dy, double t){double[] x = new double[y.Length];for (int i = 0; i < y.Length; i++, x0 += step){x[i] = x0;}return Hermite(x, y, dy, t);}

（2）埃尔米特等距插值

        /// <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="dy">一维数组，长度为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 Hermite(double x0, double step, double[] y, double[] dy, double t){int n = y.Length;double z = 0.0;// 循环插值for (int i = 1; i <= n; i++){double s = 1.0;double q = x0 + (i - 1) * step;double p;for (int j = 1; j <= n; j++){p = x0 + (j - 1) * step;if (j != i){s = s * (t - p) / (q - p);}}s = s * s;p = 0.0;for (int j = 1; j <= n; j++){if (j != i){p = p + 1.0 / (q - (x0 + (j - 1) * step));}}q = y[i - 1] + (t - q) * (dy[i - 1] - 2.0 * y[i - 1] * p);z = z + q * s;}return (z);}