C#实现分段三次Hermite插值
目录
一、Hermite插值介绍
1、功能说明
2、数学方法
二、代码实现
1、CubicHermiteInterpolator类封装
2、应用示例
三、导数值的获取方式
1、数学方法介绍
2、代码应用示例
四、其它封装的分段三次Hermite插值类
1、方式一
（1）封装代码
（2）应用示例
2、方法二
（1）封装代码
（2）应用示例
一、Hermite插值介绍
1、功能说明
功能：实现 三次埃尔米特插值（Cubic Hermite Interpolation），用于在两个已知点之间构造一条平滑曲线，保证插值函数在端点处匹配原函数的值 和一阶导数。
核心方法 ：根据输入的坐标点和导数，计算插值结果。
2、数学方法
教材《数值分析》-第五版 作者：李庆扬
二、代码实现
1、CubicHermiteInterpolator类封装
public class CubicHermiteInterpolator{/// <summary>/// x数值数组/// </summary>private readonly double[] _x;/// <summary>/// y轴数值数组/// </summary>private readonly double[] _y;/// <summary>/// y'导数值数组/// </summary>private readonly double[] _yPrime;/// <summary>/// 构造函数/// </summary>/// <param name="x">x数值数组</param>/// <param name="y">y数值数组</param>/// <param name="yPrime">y'导数值数组</param>public CubicHermiteInterpolator(double[] x, double[] y, double[] yPrime){ValidateInput(x, y, yPrime);_x = (double[])x.Clone();_y = (double[])y.Clone();_yPrime = (double[])yPrime.Clone();}/// <summary>/// 数据检验/// </summary>/// <param name="x">x数值数组</param>/// <param name="y">y数值数组</param>/// <param name="yPrime">y'导数值数组</param>private void ValidateInput(double[] x, double[] y, double[] yPrime){if (x == null || y == null || yPrime == null)throw new ArgumentNullException("Input arrays cannot be null.");if (x.Length != y.Length || x.Length != yPrime.Length)throw new ArgumentException("All input arrays must have the same length.");if (x.Length < 2)throw new ArgumentException("At least two points are required for interpolation.");for (int i = 0; i < x.Length - 1; i++){if (x[i] >= x[i + 1])throw new ArgumentException("x must be strictly increasing.");}}/// <summary>/// 经过插值计算后，x值对应输出y值/// </summary>/// <param name="x">输入x值</param>/// <returns>输出插值后的y值</returns>public double Evaluate(double x){int segmentIndex = FindSegmentIndex(x);double x0 = _x[segmentIndex];double x1 = _x[segmentIndex + 1];double y0 = _y[segmentIndex];double y1 = _y[segmentIndex + 1];double y0Prime = _yPrime[segmentIndex];double y1Prime = _yPrime[segmentIndex + 1];double h = x1 - x0;double t = (x - x0) / h;// Calculate basis functionsdouble t2 = t * t;double t3 = t2 * t;double h00 = 2 * t3 - 3 * t2 + 1;double h10 = -2 * t3 + 3 * t2;double h01 = t3 - 2 * t2 + t;double h11 = t3 - t2;// Compute interpolated valuereturn y0 * h00 + y1 * h10 + y0Prime * h01 * h + y1Prime * h11 * h;}/// <summary>/// 检查输入的x值，有没有超出数组范围/// </summary>/// <param name="x"></param>/// <returns></returns>private int FindSegmentIndex(double x){if (x < _x[0] || x > _x[_x.Length - 1])throw new ArgumentOutOfRangeException("x", "x is outside the interpolation range.");int low = 0;int high = _x.Length - 1;int segmentIndex = 0;while (low <= high){int mid = (low + high) / 2;if (_x[mid] <= x){segmentIndex = mid;low = mid + 1;}else{high = mid - 1;}}// Handle case where x is exactly the last pointif (segmentIndex == _x.Length - 1)segmentIndex--;return segmentIndex;}}2、应用示例
private void button1_Click(object sender, EventArgs e){// 示例数据double[] x = { 0, 1, 2, 3 };double[] y = { 0, 1, 0, -1 };double[] yPrime = { 1, 0, -1, 0 };var interpolator = new CubicHermiteInterpolator(x, y, yPrime);// 测试插值double y0 = interpolator.Evaluate(0.0); //输出0.0double y1 = interpolator.Evaluate(0.5); //输出0.625double y2 = interpolator.Evaluate(1.0); //输出1.0double y3 = interpolator.Evaluate(1.5); //输出0.625double y4 = interpolator.Evaluate(2.0); //输出0.0double y5 = interpolator.Evaluate(2.5); //输出-0.625double y6 = interpolator.Evaluate(3.0); //输出-1}三、导数值的获取方式
1、数学方法介绍
代码
// 计算 yPrime[1]（x=1 处的导数）
double h = x[2] - x[1]; // h = 1
yPrime[1] = (y[2] - y[0]) / (2 * h); // (0 - 0)/(2*1) = 0// 计算 yPrime[2]（x=2 处的导数）
yPrime[2] = (y[3] - y[1]) / (2 * h); // (-1 - 1)/(2*1) = -1// 边界点 x=0 用向前差分
yPrime[0] = (y[1] - y[0]) / h; // (1 - 0)/1 = 1// 边界点 x=3 用向后差分
yPrime[3] = (y[3] - y[2]) / h; // (-1 - 0)/1 = -1
这种设定可以生成一个平滑的插值曲线
2、代码应用示例
代码
private void button2_Click(object sender, EventArgs e){double[] x = { 0, 1, 2, 3 };double[] y = { 0, 1, 0, -1 };double[] yPrime = new double[4];// 计算 yPrime[1]（x=1 处的导数）double h = x[2] - x[1]; // h = 1yPrime[1] = (y[2] - y[0]) / (2 * h); // (0 - 0)/(2*1) = 0// 计算 yPrime[2]（x=2 处的导数）yPrime[2] = (y[3] - y[1]) / (2 * h); // (-1 - 1)/(2*1) = -1// 边界点 x=0 用向前差分yPrime[0] = (y[1] - y[0]) / h; // (1 - 0)/1 = 1// 边界点 x=3 用向后差分yPrime[3] = (y[3] - y[2]) / h; // (-1 - 0)/1 = -1var interpolator = new CubicHermiteInterpolator(x, y, yPrime);// 测试插值double y0 = interpolator.Evaluate(0.0); //输出0.0double y1 = interpolator.Evaluate(0.5); //输出0.625double y2 = interpolator.Evaluate(1.0); //输出1.0double y3 = interpolator.Evaluate(1.5); //输出0.625double y4 = interpolator.Evaluate(2.0); //输出0.0double y5 = interpolator.Evaluate(2.5); //输出-0.5double y6 = interpolator.Evaluate(3.0); //输出-1}四、其它封装的分段三次Hermite插值类
1、方式一
（1）封装代码
public class CubicHermiteInterpolator{private double[] xPoints; // 已知的x点private double[] yPoints; // 对应的y点（函数值）private double[] yDerivatives; // 对应的y'点（导数值）public CubicHermiteInterpolator(double[] xPoints, double[] yPoints, double[] yDerivatives){if (xPoints.Length != yPoints.Length || xPoints.Length != yDerivatives.Length)throw new ArgumentException("xPoints, yPoints, and yDerivatives must have the same length.");this.xPoints = xPoints;this.yPoints = yPoints;this.yDerivatives = yDerivatives;}public double Interpolate(double x){int n = xPoints.Length;if (n == 0) throw new InvalidOperationException("No data points available.");if (x <= xPoints[0]) return yPoints[0]; // 边界情况处理if (x >= xPoints[n - 1]) return yPoints[n - 1]; // 边界情况处理int i = 0;while (x > xPoints[i + 1]) i++; // 找到正确的区间索引idouble h = xPoints[i + 1] - xPoints[i];double t = (x - xPoints[i]) / h;double t2 = t * t;double t3 = t2 * t;double a = 2 * t3 - 3 * t2 + 1; // a(t) = 2t^3 - 3t^2 + 1double b = -2 * t3 + 3 * t2; // b(t) = -2t^3 + 3t^2double c = t3 - 2 * t2 + t; // c(t) = t^3 - 2t^2 + tdouble d = t3 - t2; // d(t) = t^3 - t^2return a * yPoints[i] + b * yPoints[i + 1] + (c * h * yDerivatives[i] + d * h * yDerivatives[i + 1]);}}（2）应用示例
private void Form1_Load(object sender, EventArgs e){double[] xPoints = { 0, 1, 2, 3 }; // x坐标点double[] yPoints = { 0, 1, 4, 9 }; // y坐标点（函数值）double[] yDerivatives = { 1, 0, -3, 0 }; // 导数值（例如，在x=1处导数为0，在x=2处导数为-3）var interpolator = new CubicHermiteInterpolator(xPoints, yPoints, yDerivatives);double result = interpolator.Interpolate(1.5); // 在x=1.5处计算插值结果Console.WriteLine("Interpolated value at x=1.5: " + result); // 应该输出接近2.875的值，因为这是从y=4到y=9的过渡点附近的值。}2、方法二
（1）封装代码
// 定义数据点结构public struct HermitePoint{public double X { get; set; }public double Y { get; set; }public double M { get; set; } // 导数值}public class CubicHermiteInterpolator{private List<HermitePoint> points;public CubicHermiteInterpolator(IEnumerable<HermitePoint> dataPoints){points = dataPoints.OrderBy(p => p.X).ToList();}// 核心插值方法public double Interpolate(double x){// 1. 边界检查if (points.Count < 2) throw new ArgumentException("至少需要两个点");if (x <= points[0].X) return points[0].Y;if (x >= points.Last().X) return points.Last().Y;// 2. 查找所在区间int index = FindSegmentIndex(x);HermitePoint p0 = points[index];HermitePoint p1 = points[index + 1];// 3. 计算归一化参数tdouble deltaX = p1.X - p0.X;double t = (x - p0.X) / deltaX;// 4. 计算Hermite基函数double t2 = t * t;double t3 = t2 * t;double h00 = 2 * t3 - 3 * t2 + 1;double h10 = t3 - 2 * t2 + t;double h01 = -2 * t3 + 3 * t2;double h11 = t3 - t2;// 5. 组合计算结果（注意导数项需要乘以deltaX）return h00 * p0.Y +h10 * deltaX * p0.M +h01 * p1.Y +h11 * deltaX * p1.M;}// 二分查找所在区间private int FindSegmentIndex(double x){int low = 0;int high = points.Count - 2;while (low <= high){int mid = (low + high) / 2;if (points[mid].X <= x && x < points[mid + 1].X)return mid;if (x < points[mid].X)high = mid - 1;elselow = mid + 1;}return points.Count - 2;}}（2）应用示例
private void Form1_Load(object sender, EventArgs e){// 创建插值点（包含导数值）//var points = new List<HermitePoint>//{// new HermitePoint { X = 0, Y = 0, M = 1 }, // 导数为1// new HermitePoint { X = 2, Y = 4, M = 0 }, // 导数为0// new HermitePoint { X = 5, Y = 1, M = -1 } // 数为为-1 //};var points = new List<HermitePoint>{new HermitePoint { X = 0, Y = 0, M = 1 }, //new HermitePoint { X = 1, Y = 1, M = 0 }, //new HermitePoint { X = 2, Y = 4, M = -3 }, //new HermitePoint { X = 3, Y = 9, M = 0 }};var interpolator = new CubicHermiteInterpolator(points);// 在x=1.5处插值double result = interpolator.Interpolate(1.5); //输出2.875// 在x=3处插值result = interpolator.Interpolate(3); //输出9.0}