数值分析（三） Lagrange（拉格朗日）插值法及Matlab代码实现

前言

  本篇为插值法专栏第三篇内容讲述，此章主要讲述 Lagrange（拉格朗日）插值法及matlab代码，其中也给出详细的例子让大家更好的理解Lagrange插值法
提示 之前已经介绍牛顿插值法和三次样条插值，如果没看过前两篇的可以点击以下链接阅读

1. 数值分析（一）牛顿插值法
2. 数值分析（二）三次样条插值法
3. 数值分析（二续） 三次样条插值二类边界完整matlab代码
4. 数值分析（三） Lagrange（拉格朗日）插值法及Matlab代码实现
5. 数值分析（四） Hermite（埃尔米特）插值法及matlab代码

一、Lagrange（拉格朗日）插值

  为了构造满足插值条件$p(x_i)=f(x_i),(i=0,1,2,\dots ,n)$的便于使用的插值多项式 $P(x)$ ，在介绍Lagrange插值前先补充一下 线性插值 和 抛物插值 的知识点。

1. 线性插值

   线性插值是代数插值的最简单形式。假设给定了函数 $f(x)$ 在两个互异的点的值，$x_0,x_1,y_0=f(x_0),y_1=f(x_1)$ ，现要求用线性函数$p(x)=ax+b$ 近似的代替$f(x)$。选择参数 $a$ 和 $b$, 使 $p(x_i)=f(x_i),(i=0,1)$。称这样的线性函数 $P(x)$ 为 $f(x)$ 的线性插值函数。
   线性插值的几何意义：用通过点 $A(x_0,f(x_0))$ 和 $B(x_1,f(x_1))$ 的直线近似地代替曲线 $y=f(x)$ 由解析几何知道，这条直线用点斜式表示为（如下图所示）: $$p(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)\to p(x)=\frac{x-x_1}{x_0-x_1}{y}_0+\frac{x-x_0}{x_1-x_0}{y}_1$$

为了便于推广，那么记$$l_0(x)=\frac{x-x_1}{x_0-x_1},l_1(x)=\frac{x-x_0}{x_1-x_0}$$这是一次函数，具有如下性质$$\begin{gathered} \begin{array}{cc}{l}_0(x_0)=1,& l_0(x_1)=0,\\ l_1(x_0)=0& l_1(x_1)=1\end{array} \\ {l}_0(x)+l_1(x)=1 \end{gathered}$$
那么可以记为：$$l_k(x_i)={\delta }_{ki}=\{\begin{array}{ccc}1& (i=k)& \\ 0& (i\ne k)& \end{array}$$那么 $l_0(x)$ 和 $l_1(x)$ 称为线性插值基函数，且有 $$l_k(x)=\prod _{j=0,j\ne k}^1\frac{x-x_j}{x_k-x_j},(k=0,1)$$ 则上述线性插值函数可以表示为与基函数的线性组合 $$p(x)=l_0(x)y_0+l_1(x)y_1$$

2. 抛物插值

  抛物插值又称二次插值，它也是常用的代数插值之一。设已知 $f(x)$ 在三个互异点$x_0,x_1,x_2$的函数值$y_0,y_1,y_2$，要构造次数不超过二次的多项式
$$P(x)=a_2{x}^2+a_{1}x+a_0$$ 使满足二次插值条件：$$P(x_i)=y_i(i=0,1,2)$$ 这就是二次插值问题。其几何意义是用经过3个点 $(x{}_0,y_0),(x{}_1,y_1),(x{}_2,y_2)$ 的抛物线 $y=P(x)$ 近似代替曲线 $y=f(x)$ , 如下图所示。因此也称之为抛物插值。

线性插值 和 抛物插值 两个示意图使用matlab绘制，并借助使用了 slandarer 大佬的一个绘制带箭头的坐标轴工具箱，代码将在下方给出

其中 $P(x)$ 的参数 $a_0,a_1,a_2$ 直接由插值条件决定，即 $a_0,a_1,a_2$ 满足下面的代数方程组：$$\{\begin{array}{l}{a}_0+a_1{x}_0+a_2{x}_0^2=y_0\\ a_0+a_1{x}_1+a_2{x}_1^2=y_1\\ a_0+a_1{x}_2+a_2{x}_2^2=y_2\end{array}$$ 那么其系数矩阵为：$$\left[\begin{array}{ccc}1& x_0& x_0^2\\ 1& x_1& x_1^2\\ 1& x_2& x_2^2\end{array}\right]$$ 此行列式是范德蒙行列式，当 $x_0\ne x_1\ne x_2$ 时，方程组的解唯一，为了与下一节的Lagrange插值公式比较,仿线性插值,用基函数的方法求解方程组，求二次式 $l_0(x)$，使其满足条件：$$l_0(x_0)=1,l_0(x_1)=0,l_0(x_2)=0$$，于是可以知道 $x_1,x_2$是$l_0(x)$ 的两个零点，则$$l_0(x)=c(x-x_1)(x-x_2)$$再由另一个条件$l_0(x_0)=1$解得系数为$$c=\frac{1}{(x_0-x_1)(x_0-x_2)}$$从而导出$$l_0(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}$$类似地可以构造出满足条件：$l_0(x_0)=0,l_0(x_1)=1,l_0(x_2)=0$的插值多项式$$l_1(x)=\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}$$以及满足条件$l_0(x_0)=0,l_0(x_1)=0,l_0(x_2)=1$ 的插值多项式 $$l_2(x)=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$$这样构造出来的$l_0(x),l_0(x),l_0(x)$称为抛物插值的基函数，取已知数据$y_0,y_1,y_2$ 作为线性组合系数，将基函数$l_0(x),l_1(x),l_2(x)$线性组合可得$$P(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}{y}_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}{y}_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}{y}_2$$容易看出，$P(x)$满足条件 $P(x_i)=y_i(i=0,1,2)$

3. 拉格朗日插值多项式

  通过上述描述可以发现，两个插值点可求出一次插值多项式，而三个插值点可求出二次插值多项式。插值点增加到 $n+1$ 个时，也就是通过 $n+1$ 个不同的已知点$(x_i,y_i)(i=0,1,\cdots ,n)$，来构造一个次数为n的代数多项式 $P(x)$。与推导抛物插值的基函数类似，先构造一个特殊 $n$ 次多项式 $l_i(x)$ 的插值问题，使其在各节点 $x_i$ 上满足 $$l_k(x_0)=0,\cdots ,l_k(x_{k-1})=0,l_k(x_k)=1,l_k(x_{k+1})=0,\cdots ,l_k(x_n)=0$$即$$l_k(x_i)={\delta }_{ki}=\{\begin{array}{cc}1& (i=k)\\ 0& (i\ne k)\end{array}$$由条件$l_k(x_i)=0,(i\ne k)$可知，$x_0,x_1,\cdots ,x_{k-1},x_{k+1},\cdots ,x_n$ 都是 $n$ 次 $l_k(x)$ 的零点，故可设$$l_k(x)=A_k(x-x_0)(x-x_1)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n)$$其中 $A_k$为待定常数。有条件$l_k(x_k)=1$，可求得 $A_k$为：$$A_k\prod _{j=0,j\ne k}^n(x_k-x_j)=1$$于是$$A_k=\frac{1}{{\prod }_{j=0,j\ne k}^n(x_k-x_j)}$$代入上式，得$$l_k(x)=\frac{{\prod }_{j=0,j\ne k}^n(x-x_j)}{{\prod }_{j=0,j\ne k}^n(x_k-x_j)}=\prod _{j=0,j\ne k}^n\frac{x-x_j}{x_k-x_j}$$称$l_k(x)$为关于基点 $x_i$ 的 $n$ 次插值基函数 $(i=0,1,\cdots ,n)$，那么以$n+1$个$n$次基本插值多项式$l_k(x)(k=0,1,\cdots ,n)$为基础,就能直接写出满足插值条件$P(x_i)=f(x_i)(i=0,1,2,\cdots ,n)$的$n$次代数插值多项式$$P(x)=l_0(x)y_0+l_1(x)y_1+\cdots +l_n(x)y_n$$其中每个插值基函数$l_k(x)(k=0,1,\cdots ,n)$都是$n$次值多项式,所以他们的线性组合$$P(x)=\sum _{k=0}^n{l}_k(x)y_k$$
即为$n$次拉格朗日插值多项式，并记作$L_n(x)$，引入${\omega }_{n+1}(x)$为$${\omega }_{n+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_n)$$那么得到$${{\omega }^{\prime }}_{n+1}(x_k)=(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)$$ 那么$L_n(x)$可以转化为$$L_n(x)=\sum _{k=0}^n{y}_k\frac{{\omega }_{n+1}(x)}{(x-x_k){{\omega }^{\prime }}_{n+1}(x_k)}$$

二、Lagrange插值算法及matlab代码

Matlab 版本号 2022b
根据Lagrange 插值算法构建 lagrange.m文件
为了让代码具有可输出性，让大家通俗易懂，编写了一个lab函数，并构建了lab.m文件

1. Lagrange 插值算法matlab实现

下面展示拉格朗日插值函数代码: lagrange.m（复制代码，保存此名称）

function varargout = lagrange(x0,y0,varargin)
% x0 为节点
% y0 为节点对应的函数值
% x 为插值点
m = length(x0);
L1 = zeros(m,m);
l = [];

% step 1 构造基函数lk
for k = 1 : m
    V = 1;
    for i = 1 : m
        if k ~= i
            V = conv(V,poly(x0(i))) / (x0(k) - x0(i));
        end
    end
    L1(k,:) = V;
    l = cat(1,l,poly2sym(V));
end
fprintf('基函数为：\n');

for k = 1:m
    fprintf(['l',lab(k),'(x)=%s\n'],l(k));fprintf('\n');
end

% 通过矩阵乘积求L(x)
L = y0 * l;
L1 = str2sym(['P(x)=',char(L)]);
disp(L1)
% L = simplify(L1);
S = ['$',latex(L1),'$'];
fprintf('拉格朗日多项式为:\nP(x)=%s\n',L);fprintf('\n');
% symdisp(L)
L = matlabFunction(L);
if nargin == 3
    varargout{1} = feval(L,varargin{1});
    plot(x0,y0,'bo')
    hold on
    plot(varargin{1},varargout{1},'rp')
    plot(min([varargin{1},x0]):0.01:max([varargin{1},x0]),feval(L,min([varargin{1},x0]):0.01:max([varargin{1},x0])),'k-','LineWidth',1)
    subtitle(S,'Interpreter','latex','FontSize',14)
    title('Lagrange插值','FontSize',14)
    legend('原始点','插值点','插值函数P(x)')
else
plot(x0,y0,'bo')
hold on
plot(min(x0):0.01:max(x0),feval(L,min(x0):0.01:max(x0)),'k-','LineWidth',1)
subtitle(S,'Interpreter','latex','FontSize',14)
title('Lagrange插值','FontSize',14)
legend('原始点','插值函数P(x)')
end
end

下面展示lab代码: lab.m（复制代码，保存此名称）

function label = lab(k)
name = {'₀','₁','₂','₃','₄','₅','₆','₇','₈','₉'};
label = '';
k=k-1;
if k == 0
    label = name{k+1};
else
    while 1
        if mod(k,10)==0 && floor(k/10)==0
            break
        end
        j = mod(k,10);
        k = floor(k/10);
        label = cat(2,name{j+1},label);
    end
end
end

2 实例

  话不多说，直接讲述lagrange函数使用说明

  1. lagrange(x0,y0)： 输入参数只有2个时，x0 为节点，y0 为节点对应的函数值；自动输出：基函数，拉格朗日多项式，插值函数图像
  2. y = lagrange(x0,y0,x)： 输入参数为3个时，x0 为节点，y0 为节点对应的函数值，x 为插值点；输出：y为插值点对应函数值，自动输出：基函数，拉格朗日多项式，插值函数图像

示例展示一下：
  已知 $f(x)$ 的观测数据为如下表所示，并构造Lagrange插值多项式

$x$	0	1	2	4
$f(x)$	1	9	23	3

手动计算一下，四个点可构造三次Lagrange插值多项式，其基函数为：
$$l_0(x)=\frac{(x-1)(x-2)(x-4)}{(0-1)(0-2)(0-4)}=-\frac{1}{8}{x}^3+\frac{7}{8}{x}^2-\frac{7}{4}x+1$$$$l_1(x)=\frac{(x-0)(x-2)(x-4)}{(1-0)(1-2)(1-4)}=\frac{1}{3}{x}^3-2x^2+\frac{8}{3}x$$$$l_2(x)=\frac{(x-0)(x-1)(x-4)}{(2-0)(2-1)(2-4)}=-\frac{1}{4}{x}^3+\frac{5}{4}{x}^2-x$$$$l_3(x)=\frac{(x-0)(x-1)(x-2)}{(4-0)(4-1)(4-2)}=\frac{1}{24}{x}^3-\frac{1}{8}{x}^2+\frac{1}{12}x$$那么Lagrange插值多项式为：$$\begin{array}{c}{L}_3(x)=\sum _{k=0}^3{y}_k{l}_k(x)=l_0(x)+9l_1(x)+23l_2(x)+3l_3(x)=-\frac{11}{4}{x}^3+\frac{45}{4}{x}^2-\frac{1}{2}x+1\end{array}$$

接下来利用上述lagrange(x0,y0)，检验一下，先构建一个叫demo.m文件，接下来复制此代码放入其中，切记要将demo.m`文件、lagrange.m文件和lab.m文件放在一个文件内（新手注意！！！）

x0 = [0 1 2 4];
y0 = [1 9 23 3];
lagrange(x0,y0);

直接启动运行，输出结果如下，发现和手动计算结果一摸一样

这里代码中还给出了此插值函数的图像

在上述的基础上需要求出当$x=1.5,2.5,3,3.5$时，对应的函数值是什么？复制如下代码直接运行看看结果。

x0 = [0 1 2 4];
y0 = [1 9 23 3];
x = [1.5 2.5 3 3.5];
y = lagrange(x0,y0,x)

输出函数值为：

3. 线性插值示意图代码

注意：一定需要有坐标轴工具箱函数否则会报错，请前往上述链接下载。

% 使用了CSDN: Slandarer大佬开源的带箭头的坐标函数包
% 剩余代码由本人撰写
% 此代码为CSDN文章中线性插值示意图代码
% @Author
% Copyright© 2024.3.12
% CSDN name: cugautozp
% GitHub: https://github.com/cug-auto-zp
clc
clear
close all
x=0.2:0.01:1.8;
y=x-exp(x-1);
plot(x,y,'LineWidth',1.5)
% 修改坐标轴属性
ax=gca;
ax.XLim = [0,2];
ax.YLim = [-0.6,0.3];
ax.YTick = [];
ax.XTick = [];
arrowAxes(ax)


x0 = [0.5,1.6];
y0 = x0-exp(x0-1);
hold on
for i=1:length(x0)
    plot(x0(i)*[1,1],[-0.6,y0(i)],'LineStyle','--','LineWidth',1.5,'Color',[0,0,0])
    text(x0(i)+0.05,(-0.6+y0(i))/2,['$y_',num2str(i-1),'$'],'Interpreter','latex','FontSize',15,'Color',[0,0,0])
    text(x0(i),-0.65,['$x_',num2str(i-1),'$'],'Interpreter','latex','FontSize',15,'Color',[0,0,0])
end
plot(x0,y0,'LineStyle','none','Marker','.','MarkerSize',20,'Color',[1,0,0])
t=polyfit(x0,y0,1);
plot(x,polyval(t,x),'-k','LineWidth',1.5)
text(1,0.1,'$y = f(x)$','Interpreter','latex','FontSize',15,'Color',[0.00,0.45,0.74])
text(0.1,0,'$p(x)=ax+b$','Interpreter','latex','FontSize',15,'Color',[0,0,0])
text(x0(1),y0(1)-0.05,'$A(x_0,f(x_0))$','Interpreter','latex','FontSize',15,'Color',[1,0,0])
text(x0(2),y0(2)-0.05,'$B(x_0,f(x_0))$','Interpreter','latex','FontSize',15,'Color',[1,0,0])
title('线性插值示意图','FontSize',14)
text(-0.1,-0.65,'O','FontSize',15,'Color',[0,0,0])

4. 抛物插值示意图代码

% 使用了CSDN: Slandarer大佬开源的带箭头的坐标函数包
% 剩余代码由本人撰写
% 此代码为CSDN文章中抛物插值示意图代码
% @Author
% Copyright© 2024.3.12
% CSDN name: cugautozp
% GitHub: https://github.com/cug-auto-zp

clc
clear
close all
x=0:0.1:20;
y=-0.5*(x-4).^3+15*(x-4)+6*(x-4).^2+5;
plot(x,y,'LineWidth',1.5)


ax=gca;
ax.XLim = [-1,21];
ax.YLim = [-300,300];
ax.YTick = [];
ax.XTick = [];
arrowAxes(ax)


x0 = [3,9,18];
y0 = -0.5*(x0-4).^3+15*(x0-4)+6*(x0-4).^2+5;
hold on
plot(x0,y0,'LineStyle','none','Marker','.','MarkerSize',20,'Color',[1,0,0])
text(-2,-310,'O','FontSize',15,'Color',[0,0,0])
for i=1:length(x0)
    plot(x0(i)*[1,1],[-300,y0(i)],'LineStyle','--','LineWidth',1.5,'Color',[0,0,0])
    text(x0(i)+0.3,(-300+y0(i))/2,['$y_',num2str(i-1),'$'],'Interpreter','latex','FontSize',15,'Color',[0,0,0])
    text(x0(i),-320,['$x_',num2str(i-1),'$'],'Interpreter','latex','FontSize',15,'Color',[0,0,0])
end
t=polyfit(x0,y0,2);
plot(x,polyval(t,x),'-k','LineWidth',1.5)
text(15,200,'$y = f(x)$','Interpreter','latex','FontSize',15,'Color',[0.00,0.45,0.74])
text(3,150,'$y = L_2(x)$','Interpreter','latex','FontSize',15,'Color',[0,0,0])
title('抛物插值示意图')

三、总结

  此次内容主要讲的是拉格朗日插值的原理，及根据原理利用matlab编写一个通用计算公式函数，然后举例来验证代码的正确性。为了方便，这里直接将函数图像，以及中间过程的基函数和最后的插值函数都自动输出，非常方便。如果对拉格朗日插值有兴趣的话可以更加深入的做一做一下实验，还蛮有意思的。

四、插值法专栏

专栏链接：插值法专栏，如果对你有帮助的话可以点个赞，点个订阅，我将完善此专栏

1. 数值分析（一） 牛顿插值法及matlab代码
2. 数值分析（二) 三次样条插值法matlab程序
3. 数值分析（二续） 三次样条插值二类边界完整matlab代码
4. 数值分析（三） Lagrange（拉格朗日）插值法及Matlab代码实现
5. 数值分析（四） Hermite（埃尔米特）插值法及matlab代码