拉格朗日插值算法例子、二次n次插值（matlab）

菜鸟粥粥

已于 2023-12-24 21:24:51 修改

阅读量1.9k

点赞数 19

文章标签：算法 matlab

于 2023-12-22 22:17:15 首次发布

本文链接：https://blog.csdn.net/2301_76767110/article/details/135161880

版权

拉格朗日插值算法思想

首先我们需要建立插值基函数，根据插值点序列 $\left(x_k, y_k\right)(k=0,1, \cdots, n)$ ，建立的函数如下

$l_k(x)=\frac{\left(x-x_0\right)\left(x-x_1\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_n\right)}{\left(x_k-x_0\right)\left(x_k-x_1\right) \cdots\left(x_k-x_{k-1}\right)\left(x_k-x_{k+1}\right) \cdots\left(x_k-x_n\right)}$

根据插值基函数构造拉格朗日插值多项式
$L_n(x)=\sum_{k=0}^n l_k(x) y_k$
最后把x的具体值代入，即求得相应的函数值。
输入插值节点数 $(n+1)$ ，插值点序列 $\left(x_k, y_k\right)$ ，然后计算函数点x，输出参数。
拉格朗日插值算法步骤
第一步：输入插值点序列 $\left(x_k, y_k\right)(k=0,1, \cdots, n)$ ，令 $L_n(x)=0$ 。
第二步：对 $k=0,1,...,n$ ，计算 $l_k(x)=\prod_{\substack{i=0 \\ i \neq \hat{k}}}^n \frac{x-x_i}{x_k-x_i}$ 。
第三步：计算 $L_n(x)=L_n(x)+l_k(x) y_k$ 。
第四步：输出 $L_n(x)$ 的近似值或计算失败信息，结束程序。
拉格朗日插值算法N-S流程图如下：

已知函数表如下表所示，用二次插值函数求 $f(1.54)$ 的近似值（计算结果取五位小数）

x	1.2	1.3	1.4	1.5	1.6	1.7
$f(x)$	1.244	1.406	1.602	1.837	2.121	2.465

程序：

Lagrange二次插值：

clear
clc
x = [1.4 1.5 1.6];
y = [1.602 1.837 2.121];
xs=1.54;%待求值
f=0;%拉朗二次插值结果
ftemp=1;%中间变量，临时保存各插值基函数中的项的值
ft1 = 1;%中间变量，临时保存各插值基函数的值
for i=0:2%第i个插值基函数
    for j=0:2%第j个节点
        if i~=j
            ftemp = (xs-x(j+1))/(x(i+1)-x(j+1));
        end
        ft1 = ft1 * ftemp;
        ftemp = 1;
    end
    f = f + y(i+1)*ft1;
    ft1 = 1;
end

Lagrange n次插值：

%%Lagrange拉格朗日n次插值（函数封装）
clear
clc
%输入数据
x = [1.2 1.3 1.4 1.5 1.6 1.7];
y = [1.244 1.406 1.602 1.837 2.121 2.465];
xs=1.54;%待求值
n=2;%指定插值次数，可选。
f = Lag02(xs,x,y)%不指定插值次数
%Lag02(xs,x,y,n)%指定插值次数
plot(x,y,'--',xs,f,'*')

function f = Lag02(xs,x,y,n)
%Lagrange插值
%   xs为待求值
%   x,y为列表向量
%   n为插值次数，可选，若为空则自动按“点个数-1”次插值，若不为空，自动从xs值左右选取n+1个值
%确保输入了足够的参数
if nargin < 4
    n = length(x)-1;
    if nargin <3
        error(message('参数个数不符合要求'));
    end
    if length(x) ~= length(y)
        error(message('x与y的长度不相等，快去检查是不是数据输错了'));
    end
else
    x1 = [x(find(x<xs,ceil(n/2)+1,'last')),x(find(x>xs,n-ceil(n/2)))];
    y1 = [y(find(x<xs,ceil(n/2)+1,'last')),y(find(x>xs,n-ceil(n/2)))];
    x=[];y=[];x=x1;y=y1;
end
f=0;%声明变量，拉朗插值结果
ftemp=1;%中间变量，临时保存各插值基函数中的项的值
ft1 = 1;%中间变量，临时保存各插值基函数的值
for i=0:n%第i个插值基函数
    for j=0:n%第j个节点
        if i~=j
            ftemp = (xs-x(j+1))/(x(i+1)-x(j+1));
        end
        ft1 = ft1 * ftemp;
        ftemp = 1;
    end
    f = f + y(i+1)*ft1;
    ft1 = 1;
end
end

Lagrange n次插值输出多项式：

%%Lagrange拉格朗日n次插值（函数封装）%输出插值多项式
clear
clc
%输入数据
x = [1.2 1.3 1.4 1.5 1.6 1.7];
y = [1.244 1.406 1.602 1.837 2.121 2.465];
xs=1.54;%待求值
n=2;%指定插值次数，可选。
f = Lag03(xs,x,y)%不指定插值次数
%Lag02(xs,x,y,n)%指定插值次数
%plot(x,y,'--',xs,f,'*')

function f = Lag03(p,x,y,n)
%Lagrange插值
%   xs为待求值
%   x,y为列表向量
%   n为插值次数，可选，若为空则自动按“点个数-1”次插值，若不为空，自动从xs值左右选取n+1个值
syms p;
%确保输入了足够的参数
if nargin < 4
    n = length(x)-1;
    if nargin <3
        error(message('参数个数不符合要求'));
    end
    if length(x) ~= length(y)
        error(message('x与y的长度不相等，快去检查是不是数据输错了'));
    end
else
    x1 = [x(find(x<p,ceil(n/2)+1,'last')),x(find(x>p,n-ceil(n/2)))];
    y1 = [y(find(x<p,ceil(n/2)+1,'last')),y(find(x>p,n-ceil(n/2)))];
    x=[];y=[];x=x1;y=y1;
end
f=0;%声明变量，拉朗插值结果
ftemp=1;%中间变量，临时保存各插值基函数中的项的值
ft1 = 1;%中间变量，临时保存各插值基函数的值
for i=0:n%第i个插值基函数
    for j=0:n%第j个节点
        if i~=j
            ftemp = (p-x(j+1))/(x(i+1)-x(j+1));
        end
        ft1 = ft1 * ftemp;
        ftemp = 1;
    end
    f = f + y(i+1)*ft1;
    ft1 = 1;
    simplify(f);
end
end

结果：