拉格朗日(Lagrange)曲线插值

作者：追风剑情发布于：2017-4-27 23:33 分类：Algorithms

已知若干个离散点的位置值和导数值，求经过这些点的插值多项式，是数学上熟知的Hermite插值问题。这一问题的明确提法是：已知函数f(t)在k+1个点{t_i}处的函数值和导数值 ${f^{(j)} (t_{i})}, i = 0, 1, \dots, k, j = 0, 1, \dots, m_{i} - 1$ 要求确定一个 $N = m_{0} + m_{1} + \dots + m_{k} - 1$ 次的多项式P(t)，满足下面的插值条件： $P^{(j)} (t_{i}) = f^{(j)} (t_{i})$

这样的多项式P(t)就是对于函数f(t)的Hermite插值多项式。数学上已经证明，这样的多项式是存在且唯一的。

一、Lagrange插值公式

当 $j = 0, m_{0} = m_{1} = \dots = m_{k} = 1$ 时， $f^{(j)} (t_{i})$ 为函数f(t)本身，这个问题就是熟知的Lagrange插值问题，即已知f(t)在k+1个点上的函数值 $f (t_{i})$ ，求一个k次多项式使之满足 $P (t_{i}) = f (t_{i}), i = 0, 1, \dots, k$ 。

式中，g_i(t)是混合函数。

使用此混合函数，使 $t = t_{i}$ 时，曲线P(t)通过给定的型值点 $f (t_{i})$ 。

任何类型的曲线都可以通过其控制点的线性组合来表示。

例如，当k=2，m0=m1=m2=1时，已知函数f(t)在三个点t0、t1和t2的函数值f(t0)、f(t1)和f(t2)，则二次多项式P(t)为

$P (t) = f (t_{0}) g_{0} (t) + f (t_{1}) g_{1} (t) + f (t_{2}) g_{2} (t)$

式中，混合函数如下： ${\begin{array}{r} g_{0} (t) = \frac{(t - t_{1}) (t - t_{2})}{(t_{0} - t_{1}) (t_{0} - t_{2})} \\ g_{1} (t) = \frac{(t - t_{0}) (t - t_{2})}{(t_{1} - t_{0}) (t_{1} - t_{2})} \\ g_{2} (t) = \frac{(t - t_{0}) (t - t_{1})}{(t_{2} - t_{0}) (t_{2} - t_{1})} \end{array}$

再如，设表示一条曲线的某个函数f(t)在四点 $t_{0} 、 t_{1} 、 t_{2} 、 t_{3}$ 的函数值 $f (t_{0}) 、 f (t_{1}) 、 f (t_{2}) 、 f (t_{3})$ 。根据Lagrange插值法，则三次多项式P(t)可表示为 $P (t) = f (t_{0}) g_{0} (t) + f (t_{1}) g_{1} (t) + f (t_{2}) g_{2} (t) + f (t_{3}) g_{3} (t)$ 式中 ${\begin{array}{r} g_{0} (t) = \frac{(t - t_{1}) (t - t_{2}) (t - t_{3})}{(t_{0} - t_{1}) (t_{0} - t_{2}) (t_{0} - t_{3})} \\ g_{1} (t) = \frac{(t - t_{0}) (t - t_{2}) (t - t_{3})}{(t_{1} - t_{0}) (t_{1} - t_{2}) (t_{1} - t_{3})} \\ g_{2} (t) = \frac{(t - t_{0}) (t - t_{1}) (t - t_{3})}{(t_{2} - t_{0}) (t_{2} - t_{1}) (t_{2} - t_{3})} \\ g_{3} (t) = \frac{(t - t_{0}) (t - t_{1}) (t - t_{2})}{(t_{3} - t_{0}) (t_{3} - t_{1}) (t_{3} - t_{2})} \end{array}$

一般地，对于k+1个点 $P_{i} ， (i = 0, 1, \dots, k)$ ，若曲线 $P (t) = \sum_{i = 0}^{k} P_{i} g_{i} (t)$ 表达式中 $g_{i} (t)$ 满足
1） $g_{i} (t) ， (i = 0, 1, \dots, k)$ 是连续的；
2） $\sum_{i = 0}^{k} g_{i} (t) = 1$ 。
则 $g_{i} (t) ， (i = 0, 1, \dots, k)$ 称为混合(调和)函数或基函数，k+1个点 $P_{i} ， (i = 0, 1, \dots, k)$ 称为控制点。如果混合函数 $g_{i} (t) ， (i = 0, 1, \dots, k)$ 使曲线经过这k+1个点 $P_{i} ， (i = 0, 1, \dots, k)$ ，则这k+1个点 $P_{i} ， (i = 0, 1, \dots, k)$ 也称为型值点。混合函数定义了如何把控制点的向量值混合在一起，为描述曲线提供了很好的抽象。这样的函数很多，例如，直线段方程 $P (t) = P_{0} (1 - t) + P_{1} t$ 中 $g_{0} (t) = 1 - t ， g_{1} (t) = t$ 就是混合函数，它定义了型值点的线性混合（加权平均）。任何类型的曲线都可以通过其控制点的线性组合来表示，其中权值可以作为关于自由参数的混合函数来计算。

代码实现(C#版)

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
 
namespace TestLine
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
 
        private void Form1_Load(object sender, EventArgs e)
        {
 
        }
 
        protected override void OnPaint(PaintEventArgs e)
        {
            LagrangeCurve curve = new LagrangeCurve();
            curve.input.Add(new Vector2(0, 0));
            curve.input.Add(new Vector2(10, 10));
            curve.input.Add(new Vector2(20, 20));
            curve.input.Add(new Vector2(30, 30));
            curve.input.Add(new Vector2(40, 40));
 
            curve.input.Add(new Vector2(50, 40));
            curve.input.Add(new Vector2(60, 30));
            curve.input.Add(new Vector2(70, 20));
            curve.input.Add(new Vector2(80, 10));
            curve.input.Add(new Vector2(90, 0));
 
            for (float x = 0; x <= 90; x++)
            {
                float y = curve.Lerp(x);
                PaintPoint(e, x, y);
            }
        }
 
        public void PaintPoint(PaintEventArgs e, float x, float y)
        {
            Graphics g = e.Graphics;
            Pen myPen = new Pen(Color.Red, 2);
            g.DrawEllipse(myPen, x, y, 2, 2);
        }
    }
 
    public class Vector2
    {
        public float x;
        public float y;
 
        public Vector2(float x, float y)
        {
            this.x = x;
            this.y = y;
        }
    }
 
    /// <summary>
    /// 拉格朗日曲线插值类
    /// </summary>
    public class LagrangeCurve
    {
        //输入点(即，曲线要经过的点)
        public List<Vector2> keyPoint = new List<Vector2>();
 
        public float Lerp(float x)
        {
            float y = 0;
            //构建input.Count个混合函数
            for (int i = 0; i < keyPoint.Count; i++)
            {
                float yi = keyPoint[i].y;
                float xi = keyPoint[i].x;
                float yj = 1;
                //每个混合函数需要遍历所有输入点(根据混合函数公式)
                for (int j = 0; j < keyPoint.Count; j++)
                {
                    if (i != j)
                    {
                        float xj = keyPoint[j].x;
                        yj *= (x - xj) / (xi - xj);
                    }
                }
                y += yi * yj;
            }
            return y;
        }
    }
}