拉格朗日(Lagrange)曲线插值

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

已知若干个离散点的位置值和导数值，求经过这些点的插值多项式，是数学上熟知的Hermite插值问题。这一问题的明确提法是：已知函数f(t)在k+1个点{t_i}处的函数值和导数值$\{f^{(j)}(t_i)\},i=0,1,\cdots,k,\quad j=0,1,\cdots,m_i-1$要求确定一个$N=m_0+m_1+\cdots+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= \cdots =m_k=1$时，$f^{(j)}(t_i)$为函数f(t)本身，这个问题就是熟知的Lagrange插值问题，即已知f(t)在k+1个点上的函数值$f(t_i)$，求一个k次多项式使之满足$P(t_i)=f(t_i),\quad i=0,1,\cdots,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{equation} \left\{ \begin{aligned} 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{aligned} \right. \end{equation} $$

再如，设表示一条曲线的某个函数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{equation} \left\{ \begin{aligned} 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{aligned} \right. \end{equation} $$

一般地，对于k+1个点$P_i，(i=0,1,\cdots,k)$，若曲线 $$ P(t)=\displaystyle\sum_{i=0}^{k}P_ig_i(t) $$ 表达式中$g_i(t)$满足
1）$g_i(t)，(i=0,1,\cdots,k)$是连续的；
2）$\displaystyle\sum_{i=0}^{k}g_i(t)=1$。
则$g_i(t)，(i=0,1,\cdots,k)$称为混合(调和)函数或基函数，k+1个点$P_i，(i=0,1,\cdots,k)$称为控制点。如果混合函数$g_i(t)，(i=0,1,\cdots,k)$使曲线经过这k+1个点$P_i，(i=0,1,\cdots,k)$，则这k+1个点$P_i，(i=0,1,\cdots,k)$也称为型值点。混合函数定义了如何把控制点的向量值混合在一起，为描述曲线提供了很好的抽象。这样的函数很多，例如，直线段方程$P(t)=P_0(1-t)+P_1t$中$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;
        }
    }
}

运行测试