牛顿插值算法举例，matlab求解

最新推荐文章于 2025-03-13 16:46:02 发布

菜鸟粥粥

最新推荐文章于 2025-03-13 16:46:02 发布

阅读量1.4k

点赞数 12

文章标签：算法 matlab 开发语言

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

版权

本文详细介绍了牛顿插值算法，包括计算步骤、使用牛顿插值法求解给定函数的近似值，以及提供了一个MATLAB函数实现，展示了如何构造差商表和n次插值多项式，最终给出了1.9441的插值结果作为实例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目：

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

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

牛顿插值算法思想：

已知插值点序列 $\left(x_k, y_k\right)(k=0,1, \cdots, n)$ ，首先计算 $k$ 阶差商

$f\left[x_i, x_{i+1}, \cdots, x_{i+k}\right]=\frac{f\left[x_{i+1}, x_{i+2}, \cdots, x_{i+k}\right]-f\left[x_i, x_{i+1}, \cdots, x_{i+k-1}\right]}{x_{i+k}-x_i}$

利用差商建立牛顿插值多项式

$\begin{aligned} N_n(x)= & f\left(x_0\right)+f\left[x_0, x_1\right]\left(x-x_0\right)+f\left[x_0, x_1, x_2\right]\left(x-x_0\right)\left(x-x_1\right)+\cdots \\ & +f\left[x_0, x_1, \cdots, x_n\right]\left(x-x_0\right)\left(x-x_1\right) \cdots\left(x-x_{n-1}\right)\end{aligned}$

最后把 $x$ 的具体值代入，求得相应的函数值。

步骤：

第一步：输入插值点序列 $\left(x_k, y_k\right)(k=0,1, \cdots, n)$ ，令 $N_n(x)=0$ 。

第二步：对 $k=0,1,...,n$ ，计算 $f(x)$ 的各阶差商 $f\left[x_0, x_1, \cdots x_k\right]$ 。

第三步：计算函数值

$\begin{gathered}N_n(x)=f\left(x_0\right)+f\left[x_0, x_1\right]\left(x-x_0\right)+f\left[x_0, x_1, x_2\right]\left(x-x_0\right)\left(x-x_1\right)+\cdots \\ +f\left[x_0, x_1, \cdots, x_n\right]\left(x-x_0\right)\left(x-x_1\right) \cdots\left(x-x_{n-1}\right)\end{gathered}$

第四步：输出 $N_n(x)$ 的近似值或计算失败信息，结束程序。

牛顿插值算法N-S流程图：

程序：

%%Newton牛顿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 = NewT03(xs,x,y)%不指定插值次数
%NewT03(xs,x,y,n)%指定插值次数
plot(x,y,'--',xs,f,'*')

function f = NewT03(xs,x,y,n)
%Newton插值
%   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
%%构造差商表
ZC=[x',y'];%将x和y按列放到一起
CS=zeros(n+1);%构造差商表的容器，包括y值列
CS(:,1)=ZC(:,2);
for i=2:n+1%行
    for j=2:n+1%列
        if i>=j
            CS(i,j)=(CS(i-1,j-1)-CS(i,j-1))/(ZC(i-j+1,1)-ZC(i,1));
        end
    end
end
%%写n次Newton插值多项式
%提取差商表矩阵对角线元素
Dcs=diag(CS);
%开始写多项式表达式
f = Dcs(1);%f保存插值结果
ftemp=1;%中间变量，临时保存各项之中(x-x_(N))乘积的值的值
for i=2:n+1
    for j=1:i-1
        ftemp = ftemp*(xs-ZC(j,1));
    end
    f = f + ftemp*Dcs(i);
    ftemp = 1;
end
end

多项式输出

%%Newton牛顿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 = NewT04(xs,x,y)%不指定插值次数
%NewT03(xs,x,y,n)%指定插值次数
%plot(x,y,'--',xs,f,'*')

function f = NewT04(p,x,y,n)
%输出插值多项式
%Newton插值
%   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
%%构造差商表
ZC=[x',y'];%将x和y按列放到一起
CS=zeros(n+1);%构造差商表的容器，包括y值列
CS(:,1)=ZC(:,2);
for i=2:n+1%行
    for j=2:n+1%列
        if i>=j
            CS(i,j)=(CS(i-1,j-1)-CS(i,j-1))/(ZC(i-j+1,1)-ZC(i,1));
        end
    end
end
%%写n次Newton插值多项式
%提取差商表矩阵对角线元素
Dcs=diag(CS);
%开始写多项式表达式
f = Dcs(1);%f保存插值结果
ftemp=1;%中间变量，临时保存各项之中(x-x_(N))乘积的值的值
for i=2:n+1
    for j=1:i-1
        ftemp = ftemp*(p-ZC(j,1));
    end
    f = f + ftemp*Dcs(i);
    ftemp = 1;
    simplify(f);
end
end

结果分析：

Newton n次插值

1.9441