数值分析·学习 | 牛顿插值法matlab实现

不变因子

已于 2022-12-26 22:24:44 修改

阅读量1.7w

点赞数 17

文章标签：学习

于 2022-07-03 09:24:37 首次发布

本文链接：https://blog.csdn.net/qq_59235603/article/details/125574364

版权

前言

本篇内容为个人所学知识分享

一、牛顿插值法是什么？

由于利用插值基函数得到的拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为重要，但当插值节点增减时，计算要全部重新进行，甚为不便，为了计算方便，重新设计出逐次生成插值多项式的牛顿插值

1.均差下的牛顿插值

当n=1时，此时插值多项式记为 $P_1(x)$ ，它满足条件 $P_1(x_0)=f(x_0),P_1(x_1)=f(x_1)$

用点弦式表示为

$P_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0),$

它可看成是零次插值 $P_0(x)=f(x_0)$ 的修正，即 $P_1(x)=P_0(x)+a_1(x-x_0),$

其中 $a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}$ 是函数 $f(x)$ 的差商，在考察三个节点的二次插值, $P_2(x)$ ,

它满足条件

$P_2(x_0)=f(x_0),P_2(x_1)=f(x_1),P_2(x_2)=f(x_2),$

可表示为

$P_2(x_0)=P_1(x)+a_2(x-x_0)(x-x_1)$

显然它满足条件 $P_2(x_0)=f(x_0)$ 及 $P_2(x_1)=f(x_1)$ .令 $P_2(x_2)=f(x_2)$ ，则得

$a_2=\frac{P_2(x_2)-P_1(x_2)}{(x_2-x_0)(x_2-x_1)}=\frac{\frac{f(x_2)-f(x_0)}{x_2-x_0}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_1}.$

系数 $a_2$ 是函数f的“差商的差商”，一般情景已知f在插值点 $x_i(i=0,1,\cdot \cdot \cdot ,n)$ 上的值为 $f(x_i)(i=0,1,\cdot \cdot \cdot ,n)$ ,要求n次插值多项式 $P_n(x)$ 满足条件

$P_n(x)=f(x_i),i=0,1,\cdot \cdot \cdot ,n,$

则 $P_n(x)$ 可表示为

$P_n(x)=a_0+a_1(x-x_0)+\cdot \cdot \cdot +a_n(x-x_0)\cdot \cdot \cdot (x-x_{n-1}),$

其中 $a_0,a_1,\cdot \cdot \cdot ,a_n$ 为待定系数，与拉格朗日插值不同，这里的 $P_n(x)$ 是由基函数 $\left \{ 1,x-x_0, \cdot \cdot \cdot ,(x-x_0)\cdot \cdot \cdot (x-x_{n-1})\right \}$ 逐次递推得到的，

2.为了给出 $a_i$ 的表达式，引入均差的概念

$f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}$ 称为函数 $f(x)$ 关于点 $x_0,x_k$ 的一阶均差。

$f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}$ 称为 $f(x)$ 的二阶均差。

$f[x_0,x_1,\cdot \cdot \cdot ,x_k]=\frac{f[x_0,\cdot \cdot \cdot ,x_{k-2},x_k]-f[x_0,x_1,\cdot \cdot \cdot ,x_{k-1}]}{x_k-x_{k-1}}$ 为 $f(x)$ 的k阶均差

则牛顿插值在均差形式下的公式为

$P_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdot \cdot \cdot +f[x_0,x_1,\cdot \cdot \cdot ,x_n](x-x_0)\cdot \cdot \cdot (x-x_{n-1}),$

3.差分形式的牛顿插值公式（牛顿前插公式）

针对于等距节点，即 $x_k=x_0+kh(k=0,1,\cdot \cdot \cdot ,n)$ 的情形，

这里称h为步长,令 $f_k=f(x_k)(k=0,1,\cdot \cdot \cdot ,n),$

称 $\Delta f_k=f_{k+1}-f_k$ 为 $x_k$ 处以h为步长的一阶差分，

$\Delta ^{2}f_k=\Delta f_{k+1}- \Delta f_k$ 为 $x_k$ 处的二阶差分，

$\Delta ^{n}f_k=\Delta ^{n-1}f_{k+1}-\Delta ^{n-1}f_{k}$ 为 $x_k$ 处的n阶差分，

且均差和差分具有以下的关系

$f[x_k,x_{k+1}]=\frac{f_{k+1}-f_k}{x_{k+1}-x_k}=\frac{\Delta f_k}{h},$

$f[x_k,x_{k+1},x_{k+2}]=\frac{f[x_{k+1},x_{k+2}]-f[x_k,x_{k+1}]}{x_{k+2}-x_k}=\frac{1}{2h^2}\Delta ^2f_k,$

$f[x_k,\cdot \cdot \cdot ,x_{k+m}]=\frac{1}{m!}\frac{1}{h^m}\Delta^mf_k,m=1,2,\cdot \cdot \cdot ,n.$

使用差分代替均差，并令x=x_0+th,则得牛顿前插公式

$P_n(x_0+th)=f_0+t\Delta f_0+\frac{t(t-1)}{2!}\Delta^2 f_0+\cdot \cdot \cdot +\frac{t(t-1)\cdot\cdot\cdot(t-n+1)}{n!}\Delta^n f_0,$

三、matlab实现代码

1.生成牛顿均差表

function table=Newtonian_mean_difference_table(x,y)
%功能：生成牛顿均差表
%输入：x为插值节点，y为插值节点对应的值
%输出:均差表table,
table=[x',y'];
[M,N]=size(table);
for n=3:M+1%循环作用：从左到右依次生成均差表
    for m=n-1:M%循环作用：在每一列里自上而下生成均差表
    table(m,n)=(table(m,n-1)-table(m-1,n-1))/(table(m,1)-table(m-n+2,1));
    end
end
end

2.牛顿插值多项式

function y0=Newtonian_interpolated(x,y,x0)
%功能：牛顿插值多项式求解
%输入：x为插值节点，y为插值节点对应的值，x0为计算点集
%输出：x0处值的集合y0
y0=onex(1,length(x0));
for k=1:length(x0)
    table=Newtonian_mean_difference_table(x,y);%调用均差表函数
    y0(k)=y(1);
    for m=2:length(x)%循环作用：均差系数（table对角线的上斜对角线元素)×(x-x0)(x-x1)···(x-xn-1)
        product=1;
        for n=2:m%循环作用：构成(x-x0)(x-x1)···(x-xn-1)
            product=product.*(x0(k)-x(n-1));
        end
        y0(k)=y0(k)+table(m,m+1).*product;
    end
end
end

3.生成牛顿差分表

function table=Difference_table(x,h,t,y)
%功能:构造差分表
%输入：第一个插值点x,步长h,次数t,插值点对应的y值
%输出：差分表table
x=x:h:x+h*t;
table=[x',y'];
[M,N]=size(table);
for n=3:M+1%循环作用：从左到右依次生成差分表
    for m=n-1:M%循环作用：在每一列里自上而下生成差分表
    table(m,n)=(table(m,n-1)-table(m-1,n-1));
    end
end
end

4.牛顿前插公式

function y0=Newtonian_pre-insertion(x,h,t,y)
%功能：牛顿前插公式
%输入：最小的点x0,步长h,次数t
%输出：x0对应
table=Difference_table(x,h,t,y)%调用差分表函数
y0=y(1);
for m=2:length(y)%循环作用：差分系数（table对角线的上斜对角线元素)×t(t-1)*...(t-n+1)/n!
    product=1;
    for n=1:m-1%循环作用：构成t(t-1)*...(t-n+1)/n!
        product=product.*(t-n+1)/n;
    end
    y0=y0+table(m,m+1).*product;
end
end
end