【数值分析】拉格朗日插值法Matlab代码+插值回归拟合介绍

---

前言

本文先是对插值、拟合、回归这三种看似相同的方法进行介绍与区分，其次详细介绍插值中的拉格朗日插值法，并采用三种思路方法编写其对应的Matlab代码，供大家思考。

方法一采用多层循环进行编写，码量极小，易于复刻，但并未求出插值函数；

方法二采用符号变量结合矩阵运算，完全按照拉格朗日插值法的思路逐步编写，码量一般，但易于理解；

方法三在方法二的基础上取消了符号变量的使用，码量一般，但并未求出插值函数。

---

一、插值、拟合、回归介绍

插值与拟合、回归均为在数学建模中常用的方法，均可以求得一个函数（一元函数或者多元函数）来刻画两组或者多组数据之间的关系，实际结果看起来似乎很像，但是差别还是不小的。

插值：相较于其他两个方法，插值所得出的函数一定经过全部的已知数据点，即$(x_i,y_i)\in y=f(x)$，如果是利用多项式函数进行插值，则多项式函数的最高次小于等于数据点个数。

拟合：拟合函数不一定会经过全部数据点，如果是多项式拟合，那拟合次数其实是人为设定的，在多项式次数一定的情况下，尽可能的让误差最小，这便是拟合。另外，衡量误差的办法有很多，常见的就是最小二乘法。

提示：数学建模中，如果在寻找最合适的拟合函数次数时，发现某次的误差为0，那需要格外注意此时的拟合函数次数是否就是已知点的个数，如果真是等于，那说明此时已经变成了插值，拟合的函数其实此时大概率不可信，不建议使用。

回归：与拟合类似，回归函数亦不一定经过全部数据点。拟合侧重于最后的拟合函数，而回归更侧重于回归系数，以此来描述单个或者多个自变量与因变量之间的相关性与联系，数学建模中常见的是多元线性回归。

二、拉格朗日插值法

拉格朗日插值法的主要思路是：通过已知点的$x_i$数据，构建拉格朗日基函数，再将基函数与$y_i$进行线性组合，即可得到最终的拉格朗日插值函数。

拉格朗日基函数个数与已知点的个数相同，且其构造仅与自变量$x_i$相关，构造方法如下所示：

（1）计算算子（假设已知点共有$n$个）：
$$f_i(x)=(x-x_1)(x-x_2)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)$$
显然，上式共有$n$个。

（2）构建基函数：
$$l_i(x)=\frac{f_i(x)}{f_i(x_i)}$$
（3）构建拉格朗日插值函数：
$$L(x)=\sum _{i=1}^n{l}_i(x)\cdot y_i$$
由于本文侧重于代码，故拉格朗日插值法的理论部分就到此为止，如果读者对拉格朗日插值法理论详解有兴趣，可以阅读本文第一篇参考文献，或者下方链接：

https://blog.csdn.net/wzk4869/article/details/126589384

三、代码编写

1.方法一

代码如下：

function y=lagrange(X,Y,x)
% X和Y均为n行的列矩阵，存放已知点的X、Y信息
% x为m行列向量，存放目标点的x信息
% y为m行列向量，存放每个x对应的y信息
n=length(X);m=length(x);
l=ones(m,n);
for i=1:m
    for j=1:n
        for k=1:n
            if j~=k
                l(i,j)=l(i,j).*(x(i)-X(k))./(X(j)-X(k));
            end
        end
    end
end
y=l*Y;
end

（该程序来源于网络，侵私删）

2.方法二

代码如下：

function y=lagrange(x_i,y_i,x)

%% 利用已知点的横坐标求解其基函数

% 声明符号变量
syms X

% 符号数组，存储全部的x-x_i
first1=(X-x_i);

% 利用矩阵，删去矩阵主对角线上元素
% 即可得到每一个数据点拉格朗日基函数的分子
all_x_sym=repmat(first1,length(x_i),1);

% 将主对角线上赋值为1
n=size(all_x_sym,1);
idx=logical(eye(n));
all_x_sym(idx)=1;

% 求积得到基函数的分子
frac_1=prod(all_x_sym,2);

% 求解基函数分母
frac_2=zeros(length(x_i),1);
for i=1:length(x_i)
    frac_2(i)=subs(frac_1(i),x_i(i));
end

% 得到基函数
lagrange_funs=frac_1./frac_2;
lagrange_funs=lagrange_funs.';

%% 求解lagrange插值函数
temp=lagrange_funs.*y_i;
lagrange_fun=sum(temp);
disp(fun)

%% 带入实际自变量值
y=subs(fun,x);

end

该方法在代码中逐步求出拉格朗日函数的表达式，最后将需要计算的自变量x带入表达式，最终得出其对应的因变量y。

示例函数为：
$$f(x)=\frac{1}{1+x^2}$$
选取自变量为1到10的十个整数，带入上述函数，得到因变量，以这两组数据作为插值已知数据点，求解拉格朗日插值函数，并随机生成数个自变量，利用插值函数求其因变量值，测试代码如下

clear

close all

x=1:1:10;

fun=@(x)1./(1+x.^2);

y=fun(x);

plot(x,y,'r*')

x_0=sort(rand(1,20)*9+1);

y_0=lagrange(x,y,x_0);

hold on

plot(x_0,y_0,'o','Color','b')

legend("已知数据点","利用插值函数求出的点")

求解 结果如下图所示：

求解结果很不错。

3.方法三

代码如下：

function y=lagrange(x_i,y_i,x)

%% 计算所有基函数分母

% 将每一个x_i-x_j用矩阵表示出来，主对角线为0
first1=repmat(x_i,length(x_i),1);
first2=x_i.'-first1;

% 将主对角线上赋值为1
n=size(first2,1);
idx1=logical(eye(n));
first2(idx1)=1;

% 按hang求积，得到各个基函数的分母
frac_2_part=prod(first2,2);

%% 计算所有基函数分子

% 将x扩充为三维形态，便于后续操作
xx=repmat(x,length(x_i)^2,1);
xxx=reshape(xx,length(x_i),length(x_i),length(x));


% 将每一个x-x_i用矩阵表示出来，此时仍是三维矩阵
second2=xxx-first1;

% 将主对角线上赋值为1
n=size(second2,1);
idx2=logical(eye(n));
idx3=repmat(idx2,1,length(x));
idx3=reshape(idx3,n,n,length(x));
second2(idx3)=1;


% 按hang求积，得到各个基函数的分子，仍是三维矩阵
frac_1_part=prod(second2,2);
% 删除无用的一个维度，此时为二维矩阵，每一列为每个x对应的拉格朗日基函数的pi(x-x_i)
frac_1_x=reshape(frac_1_part,n,length(x));

%% 计算最终结果

lagrange_qx=y_i.'.*frac_1_x./frac_2_part;
y=sum(lagrange_qx);

end

该方法并未求解出拉格朗日的符号表达式，但是有如下考虑：既然可以用subs函数将符号表达式中的x替换为具体的值，那一定可以在求解符号表达式时，就将x替换为实际值。

方法三就是基于上述考虑进行的编写，由于采用了高维矩阵，故所占用的空间复杂度更高。

测试代码与方法二中测试代码相同，故仅展示求解结果：

求解结果亦很不错。

---

四、总结

方法二、三均很好的利用了matlab中矩阵的特性，从某种程度上来说，方法三与方法一其实大差不差，只有方法二中，采用了符号变量，但是一般情况来讲，符号函数的运行效率并没有想象中的高，使用的时候还需细细斟酌。方法三采用了高维矩阵，复刻代码时计算需更谨慎一些。

---

制作不易，如果这篇文章对您有所帮助，还请点个免费的赞吧，感谢！

---

参考文献

[1] 李庆扬,王能超,易大义.数值分析(第五版)[M].北京:清华大学出版社,2008:23-29.