拉格朗日插值原理及实现（Python, Matlab）

拉格朗日插值原理及实现（Python, Matlab）

目录

• 拉格朗日插值原理及实现（Python, Matlab）
  • 一. 前言
  • 二. 3种形式的Lagrange插值函数推导
    • 1. 原始形态的Lagrange插值
    • 2. 第一形式Lagrange插值
    • 3. 第二形式的Lagrange插值（重心插值公式）
  • 三. 利用Python编程实现这三种Lagrange插值
  • 三. Matlab实现第二类和第三类Lagrange插值

一. 前言

Lagrange插值是利用n次多项式来拟合（n+1）个数据点从而得到插值函数的方法。（注意n次多项式的定义是未知数最高次幂为n，但是多项式系数有n+1个，因为还有个常数项）

Lagrange插值和Newton插值本质上相同，都是用(n-1)次多项式来拟合n个数据点。所以这两种插值方法得到的插值函数相同，因为多项式拟合的基本定理：同时通过n个数据点，且最高次幂小于(n-1)的多项式函数唯一。下面顺手证明一下这个重要的定理。

如果已知n+1个数据点\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)，假设\(L_1=k_0+k_1x+k_2x^2+\cdots+k_nx^n\)和\(L_2=k_0’+k_1’x+k_2'x^2+\cdots+k_n’x^n\)都是通过这n个数据点的插值函数。那么应该有\(L_1 - L_2\)通过所有\((x_1,0),(x_2,0),\cdots,(x_n,0)\)

代入这些点得到齐次线性方程组：

\[\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{bmatrix} \begin{bmatrix} k_0-k_0'\\ k_1-k_1'\\ k_2-k_2'\\ \vdots \\ k_n-k_n'\\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \]

它的系数行列式是Vandermonde行列式，所以：

\[\begin{vmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{vmatrix} = \prod_{n\ge i > j \ge 0 }(x_i-x_j) \]

因为每个点都是不同的，所以\(x_i \ne x_j,i\ne j\)，所以齐次线性方程组的系数行列式不等于0，故方程解唯一且为0解：

\[k_i - k_i' = 0\\ k_i = k_i' \\ L_1 = L_2 \]

证毕。

二. 3种形式的Lagrange插值函数推导

1. 原始形态的Lagrange插值

---

为了用n次多项式拟合n+1个数据点：\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)

Lagrange插值函数采用的方法是构建一个这样的函数：

\[L(x) = l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \tag{1} \]

也就是用一组基函数\(\{l_0(x),l_1{x}, \cdots ,l_n(x)\}\)去构建插值函数\(L(x)\)，那么不难想到这样的基函数需要满足这样的条件：

\[l_i(x_j)=\begin{cases} 0 , & j\ne i \\ 1 , & j = i \end{cases} \]

这样对于\(L(x)\)就会有：

\[\begin{split} L(x_i) &= l_0(x_i)y_0 + \cdots l_i(x_i)y_i + \cdots +l_n(x)y_n \\ &=y_i \end{split} \]

这样就实现了\(L(x)\)通过所有的数据点\((x_i,y_i)\)。接下来就构建这样的基函数\(l_i(x)\)。

首先实现\(l_i(x_j) = 0,i\ne j\) ：

\[l_i(x) = (x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1}) \cdots(x-x_n) \]

这个函数实现了\(l_i(x)\)在所有非\((x_i,y_i)\)的点处为0，但是在\(x=x_i\)处：

\[l_i(x_i) = (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \]

为了让\(l_i(x_i)=1\)，我们可以将\(l_i(x)\)除以这个值进行修正：

\[\begin{split} l_i(x_i)&= \frac{(x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}\\ &=\prod_{n\ge j \ge0,j\ne i} \frac{x-x_j}{x_i -x_j} \end{split} \tag{2} \]

将(2)代入(1)就得到了原始形态的Lagrange插值函数.

\[\begin{split} L(x) &= l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \\ &= \prod_{n\ge j \ge0,j\ne 0} \frac{x-x_j}{x_0 -x_j}y_0 +\prod_{n\ge j \ge0,j\ne 1} \frac{x-x_j}{x_1 -x_j}y_1 +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{x-x_j}{x_n -x_j}y_n \end{split} \tag{3} \]

例：已知点(1,1),(2,2),(3,3)，用原始Lagrange插值计算插值函数。

首先计算三个插值基函数：

\[l_0(x) = \frac{(x-2)(x-3)}{(1-2)(1-3)}=\frac{(x-2)(x-3)}{2}\\ l_1(x) = \frac{(x-1)(x-3)}{(2-1)(2-3)}=-{(x-1)(x-3)}\\ l_2(x) = \frac{(x-1)(x-2)}{(3-1)(3-2)}=\frac{(x-1)(x-2)}{2}\\ \]

从而得到插值函数：

\[L(X) = l_0(x)+ 2l_1(x)+3l_2(x)=\frac{(x-2)(x-3)}{2}-2{(x-1)(x-3)}+ \frac{3(x-1)(x-2)}{2} \]

插入点x=4试一下：\(L(4)=4\) 。

原始模式的Lagrange插值函数，每次计算插值点时需要进行n(n-1)次乘法，时间复杂度为\(O(n^2)\)

2. 第一形式Lagrange插值

---

为了降低计算的时间复杂度，我们对原始的Lagrange插值函数进行改进。

为了书写方便，我们令\(w_i\)为

\[\begin{split} w_i &= \prod_{n\ge j \ge0,j\ne i} (x_i -x_j) \\ &= (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \end{split} \tag{4} \]

观察(3)式，我们可以提取一个公因数\((x-x_0)(x-x_1)\cdots (x-x_n) = g(x)\)：

\[\begin{split} L(x) &= (x-x_0)(x-x_1)\cdots (x-x_n)[\prod_{n\ge j \ge0,j\ne 0} \frac{y_0}{(x_0 -x_j)(x-x_0)} +\prod_{n\ge j \ge0,j\ne 1} \frac{y_1}{(x_1 -x_j)(x-x_1)} +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{y_n}{(x_n -x_j)(x-x_n)}]\\ &=g(x)[\frac{y_0}{w_0(x-x_0)}+\frac{y_1}{w_1(x-x_1)}+\cdots+\frac{y_n}{w_n(x-x_n)}]\\ &= g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \end{split} \tag{5} \]

这样就得到了第一形式Lagrange插值。这一形式的Lagrange插值计算的流程如下：

1. 预处理：根据已知数据点，利用公式(4)计算\(w_0,w_1,\cdots,w_n\);
2. 插值：计算\(g(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\)，然后根据公式(5)，先计算中括号内的加式，然后乘以\(g(x)\)得到插值点的值，这样计算一个新的点的时间复杂度就是\(O(n)\)；
3. 补充数据点：如果数据集有更新，只需要更新\(w_i\)即可，加入一个新的点到数据集，只需要将每个\(w_i\)乘以\((x_i-x_{n+1})\)，此外再增加一个\(w_{n+1}\)。

例：已知点(1,1),(2,2),(3,3)，使用Lagrange第一形式计算插值函数。

首先计算\(w_i\)：

\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]

插值函数就是：

\[L(x)=(x-1)(x-2)(x-3)[\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}] \]

插入点x=4试一下：\(L(4)=4\) 。

3. 第二形式的Lagrange插值（重心插值公式）

---

第一形式的Lagrange插值还要计算\(g(x)\)，可以再进行优化。

第一形式：

\[L(x) = g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \tag{6} \]

为了消掉\(g(x)\)，我们取\(y=1\)这条直线上的点进行插值，即取n+1个点：\((x_0,1),(x_1,1),\cdots,(x_n,1)\)那么这n+1个点的插值函数就是：

\[L'(x)=g(x)\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)} \tag{7} \]

而\(L'(x)=1,x=x_0,x_1,\cdots,x_n\) ，所以我们可以用(6)除以(7)消去g(x)：

\[L(x)=\frac{L(x)}{1}=\frac{L(x)}{L'(x)}=\frac{\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)}}{\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)}} \tag{8} \]

这样就得到了第二形式的Lagrange插值，也称为重心插值，通常Lagrange插值采用这种形式。

它的计算过程如下：

1. 预处理：计算\(w_i\)；
2. 插值：对需要插入的点(x,y)，计算（x-x_i)（可以存到一个列表中，计算时直接取用）；
3. 补充数据点：补充新的点到数据集只需要更新\(w_i\) 。

例：已知点(1,1),(2,2),(3,3)，利用重心插值公式计算插值函数。

首先计算\(w_i\)：

\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]

得到\(L(x)\)：

\[L(x)=\frac{\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}}{\frac{1}{2(x-1)}-\frac{1}{(x-2)}+\frac{1}{2(x-3)}} \]

插入点x=4试一下：\(L(4)=4\) 。

三. 利用Python编程实现这三种Lagrange插值

import numpy as np


class LagrangeInterpolation:
    """
    There are three modes of the LagrangeInterpolation.The input data is supposed
    to be two lists.
    ---------------------------------------------------------------------------------------
    self.data : Contain the points known before.
    self.dataLength : Indicate the number of points in the data list.
    self.weight : The self.weight[i] is (xi - x1)(xi - x2)...(xi - x(i-1))(xi - x(i+1))
                  ...(xi - xn)
    self.items : The self.items[i] is (x - x1)(x - x2)...(x - x(i-1))(x - x(i+1))...
                 (x - xn)
    ---------------------------------------------------------------------------------------
    """
    def __init__(self, x, y):
        self.data = {'x': list(x), 'y': list(y)}
        self.dataLength = len(self.data['x'])
        self.weight = []
        self.items = []
        # control is a flag indicating if there is anything wrong with the data.
        self.control = True
        if len(self.data['x']) != len(self.data['y']):
            print("The length of x isn't equal to the length of y!")
            self.control = False
        else:
            self.__preprocess(order=0)

    # Appending function is used to add points to the data list.
    def data_append(self, ap):
        if ap[0] in self.data['x']:
            print("The point already exist.")
            self.control = False
        else:
            self.control = True
            self.data['x'].append(ap[0])
            self.data['y'].append(ap[1])
            self.dataLength = self.dataLength + 1
        self.__preprocess(order=1)

    """
    Preprocessing is used to update the self.weight and self.items.
    order = 0 : Initialize the self.weight.
    order = 1 : Update the self.weight when new point is added to the data list.
    order = 2 : Calculate the self.items for each point waiting for interpolation.
    """
    def __preprocess(self, order=0, x=0):
        if order == 0:
            self.weight = list(np.zeros(self.dataLength))
            for i in range(self.dataLength):
                weight_temp = 1
                for j in range(self.dataLength):
                    if i == j:
                        pass
                    else:
                        weight_temp = weight_temp * (self.data['x'][i] - self.data['x'][j])
                self.weight[i] = weight_temp
        elif order == 1:
            self.weight.append(1)
            for i in range(self.dataLength - 1):
                self.weight[i] = self.weight[i] * (self.data['x'][i] - self.data['x'][-1])
                self.weight[-1] = self.weight[-1] * (self.data['x'][-1] - self.data['x'][i])
        elif order == 2:
            self.items = list(np.zeros(self.dataLength))
            for i in range(self.dataLength):
                self.items[i] = x - self.data['x'][i]

    # The mode1 is the initial mode of Lagrange interpolation.
    def mode1(self, px):
        self.__preprocess(order=2, x=px)
        if self.control:
            dataCheck = False
            for w in self.weight:
                if w != 0:
                    dataCheck = True
                else:
                    dataCheck = False
            if dataCheck:
                py = 0.0
                for i in range(self.dataLength):
                    py_temp = 1
                    for j in range(self.dataLength):
                        if i == j:
                            pass
                        else:
                            py_temp = py_temp * self.items[j]
                    py = py + py_temp * self.data['y'][i] / self.weight[i]
                return py
            else:
                print("There is a same x!")
                return None
        else:
            return None

    def mode2(self, px):
        self.__preprocess(order=2, x=px)
        itemsProd = np.prod(self.items)
        itemsSum = 0
        for i in range(self.dataLength):
            itemsSum = itemsSum + self.data['y'][i] / (self.weight[i] * self.items[i])
        py = itemsProd * itemsSum
        return py

    def mode3(self, px):
        self.__preprocess(order=2, x=px)
        denomSum = 0
        numeSum = 0
        for i in range(self.dataLength):
            dtemp = self.weight[i] * self.items[i]
            numeSum = numeSum + self.data['y'][i] / dtemp
            denomSum = denomSum + 1 / dtemp
        py = numeSum / denomSum
        return py

demo.py :

from lagrange_interpolation import LagrangeInterpolation as lag


x = [1, 2, 3, 4]
y = [1, 2, 3, 4]

inter1 = lag(x, y)
inter1.data_append((5, 5)) #往数据集中追加一个点
z1 = inter1.mode1(6)
z2 = inter1.mode2(6)
z3 = inter1.mode3(6)
print(z1)
print(z2)
print(z3)
"""
Result：
6.0
6.000000000000002
6.0000000000000036
"""

三. Matlab实现第二类和第三类Lagrange插值

第二类lagrange插值，lagrange_type2.m

function f = lagrange_type2(samples)
% Input:
% samples: [n, 2] matrix, n is the number of samples
% every row is a sample like (x, y)
% [x1,x2,...,xn]^T = samples(:,1)
% [y1,y2,...,yn]^T = samples(:,2)
% Return:
% f: the interpolation function f(x), 
% while x is the point to interpolate, x should be a single point.
% Example:
%   samples = [1,1; 2,2; 3,3];
%   inter_f = lagrange_type1(samples);
%   y = inter_f(4) % y=4
[n,~] = size(samples);
    function y = inter_func(x)
        h = x-samples(:, 1); % [x-x1,x-x2,...,x-xn]^T
        g = prod(h, "all");
        w = zeros(n, 1);
        for i=1:n
            w_temp = samples(i, 1)-samples(:,1);
            w_temp(i) = [];
            w(i) = prod(w_temp, "all");
        end
        y = sum(g.*samples(:,2)./(w.*h), "all");
    end
f = @inter_func;
end

第三类lagrange插值，lagrange_type3.m

function f = lagrange_type3(samples)
% Input:
% samples: [n, 2] matrix, n is the number of samples
% every row is a sample like (x, y)
% [x1,x2,...,xn]^T = samples(:,1)
% [y1,y2,...,yn]^T = samples(:,2)
% Return:
% f: the interpolation function f(x), 
% while x is the point to interpolate, x should be a single point.
% Example:
%   samples = [1,1; 2,2; 3,3];
%   inter_f = lagrange_type2(samples);
%   y = inter_f(4) % y=4
[n,~] = size(samples);
    function y = inter_func(x)
        h = x-samples(:, 1); % [x-x1,x-x2,...,x-xn]^T
        w = zeros(n, 1);
        for i=1:n
            w_temp = samples(i, 1)-samples(:,1);
            w_temp(i) = [];
            w(i) = prod(w_temp, "all");
        end
        y = sum(samples(:,2)./(w.*h), "all")/sum(1./(w.*h), "all"); 
    end
f = @inter_func;
end