插值与逼近（一）：Lagrange插值

Lagrange插值

对多个数据点 $f(x_i)=y_i, \forall i=1,2,...n$ $L_n(x)=\sum_{i=0} ^n y_i*l_i(x)$
基函数$l_i(x)=\prod_{j=0,j\neq i} ^n \frac{x-x_{j}}{x_i-x_j}$

Lagrange方法推导过程较简单，但是由于其没有承继性，即每次添加一个插值点都要重复计算, 不是很实用。

下面用三种方法实现该算法

数值方法

%lagrange.m
% 利用数值结果
function y=lagrange(x0,y0,x)
    %实现Lagrange插值
    %(x0(k),y0(k))为插值结点，x为拟合点，y为拟合值
    n=length(x0); %n个结点
    m=length(x);
    if(length(x0)==n)  %正确输入对应的n组节点值
        y=zeros(1,m); %存放拟合结果
        for i=1:n  %n个插值点的循环
            temp=ones(1,m);
            for j=1:n  %计算li(x)
                if(i~=j)
                    temp=temp.*(x-x0(j))./(x0(i)-x0(j));
                end
            end
            y=y+temp*y0(i);
        end
    end
end

RMK:
1.该段程序是数值计算，可以同时求得多个点的插值，并且利用向量的运算，充分发挥MATLAB矩阵计算的优势。
2.无法得到插值多项式$L_n(x)$ 的具体表达式，使得程序没有通用性。

符号计算

%lagrange2.m
function y=lagrange2(x0,y0,x)
    syms z 
    fun=@(x)0;  %初始化
    %实现Lagrange插值，并输出多项式的系数
    %(x0(k),y0(k))为插值结点，x为拟合点，y为拟合值
    n=length(x0); %n个插值节点
    m=length(x); % m个需要计算的节点
    if(length(x0)==n)  %正确输入对应的n组节点值
        for i=1:n  %n个插值点的循环
            temp=@(x)1;
            for j=1:n  %计算li(x)
                if(i~=j)
                    temp=temp*(z-x0(j))/(x0(i)-x0(j));
                end
            end
            fun=fun+temp*double(y0(i));    %存放拟合结果
            fun=vpa(fun,4)
        end
    end
    disp('该插值多项式为')
    fun=vpa(simplify(fun),3)
    for j=1:m  %计算所求点的函数值
    y(i)=double(subs(fun,z,x(i)));
    end
end

RMK:
1.该函数在内部设置了一个符号变量，并采用符号计算求得插值函数的精确表达式，每次较复杂的函数计算后都用$vpa()$控制精度，牺牲一定的精度来提高函数运行速度（因为MATLAB的符号计算每步都是精确的分数解，会严重影响函数计算速度，但实际操作中，我们并不需要这么高的精度要求），最后用$simplify()$ 函数来化简该多项式，得到一个比较美观的结果。
2.函数的句柄可以直接对符号变量做加法

syms x
f=@(x)(x+1);
f=f+x+1;  %输出f=2*x+1 是可以运行的

3.符号计算的几个函数 $subs(), vpa()$等，之后的博客中也会详细的梳理。

附上$vpa()$的帮助文档

vpa(x) uses variable-precision floating-point arithmetic (VPA) to evaluate each element of the symbolic input x to at least d significant digits, where d is the value of the digits function. The default value of digits is 32.
vpa(x,d) uses at least d significant digits, instead of the value of digits.

对多项式的系数进行操作

%lagrange3.m
function y=lagrange3(x0,y0,x)
n=length(x0);y=0;
for i=1:n
    l=1;
    for j=1:n
        if j~=i
            l=conv(l,poly (x0(j))/(x0(i)-x0(j)));
        end
    end
    y=y+l*y0(i);
end
y=poly2sym(y,'x')

RMK:
1.该函数充分利用多项式的特点，只对多项式的系数进行处理，结合MATLAB强大的多项式函数ployval()等，简化程序
2.附上$conv()$的简单说明

w = conv(u,v) returns the convolution of vectors u and v. If u and v are vectors of polynomial coefficients, convolving them is equivalent to multiplying the two polynomials.
w = conv(u,v,shape) returns a subsection of the convolution, as specified by shape. For example, conv(u,v,’same’) returns only the central part of the convolution, the same size as u, and conv(u,v,’valid’) returns only the part of the convolution computed without the zero-padded edges.