MATLAB之拉格朗日插值

一、算法原理1、插值问题的定义当精确函数 y = f(x) 非常复杂或未知时，在区间【a，b】上一系列节点x0,x1,x2,......xn处测得函数值f(x0)、f(x1)、f(x2)......f(xn)；由此构造一个简单的简单易算的近似函数g(x)≈f(x)，满足条件g(xi)=f(xi) (i=0,1,2...,n)。这个问题就被称为插值问题。2、上一篇文章，我们介绍了多项...

天涯铭

19474人浏览 · 2020-04-13 09:16:32

天涯铭 · 2020-04-13 09:16:32 发布

一、算法原理

1、插值问题的定义

当精确函数 y = f(x) 非常复杂或未知时，在区间【a，b】上一系列节点x0,x1,x2,......xn处测得函数值f(x0)、f(x1)、f(x2)......f(xn)；

由此构造一个简单的简单易算的近似函数g(x)≈f(x)，满足条件g(xi)=f(xi) (i=0,1,2...,n)。这个问题就被称为插值问题。

其中，g(x)被称为插值函数；

x0,x1,x2,......xn位插值节点；

g(x)≈f(x)为插值条件；

区间【a,b】为插值区间；

2、上一篇文章，我们介绍了多项式插值：

多项式插值及其matlab代码

问题引入：

a、要构造插值函数，首先要确定插值多项式的次数，对于n+1个点，一般可以找到一个次数不超过n的插值多项式；

b、在n次多项式空间Pn中找一组合适的基函数 $\Phi (x)$ ,使插值函数 $g(x)=a_{0}\Phi _{0}(x)+a_{1}\Phi _{1}(x)+...+a_{n}\Phi _{n}(x)$ ；

c、不同基函数的选取方法，导致了不同的插值方法(多项式插值，牛顿插值，埃米特插值)；

例如，多项式插值，选取的基函数为自然基函数{1，x,x^2,...,x^n}

3、拉格朗日插值推导

现在我们来介绍一下拉格朗日插值。

次数n=1,一次拉格朗日插值（取两个点，插值多项式最高次数n=2-1=1）

如图所示

已知点x0、x1及其函数值y0、y1，则根据多项式插值原理，可以确定一个次数不超过n-1（2-1）的代数多项式（n为插值点的个数）g(x)=a0+a1*x，使得g(x0)=y0,g(x1)=y1。其推导过程如下；

$g(x)=y_{0}+(y_{1}-y_{0})*(x-x_{0})/(x_{1}-x_{0})$

$=(x-x_{1})*y_{0}/(x_{0}-x_{1})+(x-x_{0})*y_{1}/(x_{1}-x_{0})$

$=L_{0}(x)y_{0}+L_{1}(x)y_{1}=\sum_{i=0}^{1}L_{i}y_{i}$

式中，L0（x）、L0（x）为拉个朗日插值基函数。

表达式最后的形式即为拉格朗日插值多项式的形式。

n=2,二次拉格朗日插值（取三个点，插值多项式最高次数n=3-1=2）

已知点x0、x1、x2及其函数值y0、y1、y2，则根据多项式插值原理，可以确定一个次数不超过n-1（3-1）的代数多项式（n为插值点的个数）g(x)=a0+a1*x+a2*x^2，使得g(x0)=y0,g(x1)=y1,g(x2)=y2。写成拉个朗日擦插值的形式为：

$g(x_{i})=\sum_{k=0}^{2}L_{k}(x_{i})*y_{k}=y_{i},(i=0,1,2)$

要使g(xi)=yi，其中Lk应满足Lk（xi）=1,k=i ； Lk（xi）=0 ,k≠i；根据Lk(x)的定义，xk以外所有的结点都是Lk(x)的根，因此

$L_{0}(x_{0})=\lambda *(x-x_{1})(x-x_{2})$

$L_{1}(x_{1})=\lambda *(x-x_{0})(x-x_{2})$

$L_{2}(x_{2})=\lambda *(x-x_{0})(x-x_{1})$

又由Lk(xk) = 1，得

$\lambda _{0}=\frac{1}{(x_{0}-x_{1})*(x_{0}-x_{2})}$

$\lambda _{1}=\frac{1}{(x_{1}-x_{0})*(x_{1}-x_{2})}$

$\lambda _{2}=\frac{1}{(x_{2}-x_{0})*(x_{2}-x_{1})}$

将求得的系数带入即可求得二次插值的拉格朗日多项式

n次拉格朗日插值（取n+1个点，插值多项式最高次数n）

其拉个朗日多项式的形式为：

$g(x_{i})=\sum_{k=0}^{n}L_{k}(x_{i})*y_{k}=y_{i},(i=0,1,2,...,n)$

要使g(xi)=yi，其中Lk应满足Lk（xi）=1,k=i ； Lk（xi）=0 ,k≠i；根据Lk(x)的定义，xk以外所有的结点都是Lk(x)的根，因此

$L_{0}(x_{0})=\lambda *(x-x_{0})(x-x_{1})...(x-x_{k-1})(x-x_{k+1})...(x-x_{n})$

又由Lk(xk) = 1，得

$\lambda _{1}=\frac{1}{(x_{k}-x_{0})*(x_{k}-x_{1})...(x_{k}-x_{k-1})(x_{k}-x_{k+1})...(x_{k}-x_{n})}$

综上所述

综上所述，我们有1次，2次，推导出来n次拉个朗日插值多项式的数学公式。

二、matlab代码

%% 拉格朗日插值 (创建函数脚本方可运行)
close
clc
clear
X=[2 3 5 7 8 9];
Y=[5 9 12 16 23 35];
function D=ji(X,n)  %求拉格朗日多项式的系数子函数
syms x
%法一
% switch n
%     case 1
%         B=x-X(2:end);
%         C=X(1)-X(2:end);
%     case length(X)
%         B=x-X(1:n-1);
%         C=X(n)-X(1:n-1);
%     otherwise
%         B=x-X([1:n-1,n+1:end]);
%         C=X(n)-X([1:n-1,n+1:end]);
% end
%法二
Xn=X(n);
X(n)=[];
B=x-X;
C=Xn-X;
D=prod(B)/prod(C);%prod数组元素连乘积
end

m=length(X);
for i=1:m
    L1(i)=Y(i)*ji(X,i);
end
L=sym2poly(expand(sum(L1))) %将符号变量变为数据
x=X(1):0.1:X(end);
y=polyval(L,x); %y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值

plot(x,y)  %画拟合曲线
hold on
plot(X,Y,'o') %画散点图

开放原子开发者工作坊

开放原子开发者工作坊旨在鼓励更多人参与开源活动，与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动，如meetup、训练营等，主打技术交流，干货满满，真诚地邀请各位开发者共同参与！

更多推荐

2019 CCF CSP-J第一轮

开放原子开发者工作坊

宏观经济学马工程教材个人笔记整理

宏观经济学from Miracle_Zero文章目录宏观经济学from Miracle_Zero第二章总产出及其衡量一、国民收入核算中的六个总量1、国内生产总值GDP2、国民生产总值GNP3、国民生产净值NNP与国内生产净值NDP4、国民收入NI与个人收入PI5、个人可支配收入DPI二、国内生产总值的核算方法1、生产法2、收入法3、支出法4、贡献率的计算三、总产出核算的校正1、消费价格平减指数（