牛顿插值法matlab程序解析

牛顿插值法在MATLAB 中的实现

经过n+1个不同的插值点12n+1,,x x x …,，构造牛顿插值公式

1211231212n+112n =[,]()[,,]()()[,,]()()()N f x x x x f x x x x x x x f x x x x x x x x x -+--++---（x ）……… 注：牛顿插值法中，用到了插值公式

%我们以二次牛顿插值公式为例解析牛顿插值法的matlab 程序

function[c,d]=newpoly(x,y)

%这里x 为3个节点的横坐标组成的向量，即()123,,x x x x =，y 为纵坐标的组成向量，即()()()()123,,y f x f x f x =

%c 为所得的牛顿插值多项式的系数组成的向量

n=length(x);

%测量向量x 的长度，即向量x 中元素i x 的个数，赋值给n ，所以n=3，注：这里的“n ”仅为变量，和公式中的次数n 不一样

d=zeros(n,n); d=zeros(3,3)

%把变量d 定义为一个n 行，n 列的零矩阵,此矩阵用来储存各阶差商，格式完全等同于书中21页的表 d(:,1)=y ’;

%此句是把向量y 的转置，即123()()()f x y f x f x ⎛⎫ ⎪= ⎪ ⎪⎝

⎭，赋值给零矩阵d 的第一列

%下面运用两个for 循环来构造书中21页的差商表

%第一个循环（父循环），循环变量为k

for k=2:n

%用来表示零矩阵d 中的第几行

%第二个循环（父循环），循环变量为k

for j=k:n

%用来表示零矩阵d 中的第几列

d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));

%差商公式，其中d(k,j)表示零矩阵d 中的第k 行，第j 列的元素，d(k,j-1)，d(k-1,j-1)等也类似，它们代表的元素随着双循环而变化，x(k-1)表示1k x -，这种计算差商的方法是根据差商表的排列位置而得来，具体解释见下面。

end

end

%下面以二次牛顿插值公式为例解析双循环构造差商表，让我们先来看看构造好的差商表

121232312333

()

()

[,]

()[,][,,]X f x d f x f x x f x f x x f x x x ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦

%然后我们依旧用大括号的形式表示构造各阶差商的步骤

[][]212121213223323223211233131()()(2,1)(1,1)2,(2,2)[,],2()()(3,1)(2,1)3,(3,2)[,]x ,,(3,2)(2,2)33(3,3)[,,]f x f x d d k d f x x x x x x j f x f x d d k d f x x x x x x f x f x x d d j k d f x x x x x x x ⎧--⎧====⎪⎪--⎪⎪=⎨⎪--⎪⎪====⎨⎪--⎩⎪⎪⎧--=====⎪⎨--⎪⎩⎩

， c=d(n,n); 即 c=d(3,3)

%此句是把零矩阵d 中的对角线元素给了向量c ， c 的向量长度即为3

%下面的循环为子循环，循环变量为k ，用于构造牛顿插值公式

for k=(n-1):(-1):1; k=2:-1:1 %这句话是说循环变量k 从2循环到1

c(k+1,k+1)=conv(c(k+1,k+1),poly(x(k))); %conv 为相乘函数，此句k 每循环一次就是把c 和以k x 为单根的多项式相乘起来，即()k c x x ⨯-，k 循环结束后即为11n ()()()k c x x x x x x +⨯-⨯⨯-⨯⨯-……，注初值c=d(3,3). m=length(c); %此句就是统计向量c 的长度，其实我们把它认定为n 循环到了哪次，所以m=k+1 c(m)=c(m)+d(k,k) %此句是循环中的重点，是一个著名的循环,其中c （m ）=c （k+1,k+1）； end

%下面我们以二次牛顿插值公式为例用大括号的形式表示循环构造牛顿插值公式的过程 ()()12322123212112321122,13,(3,3)[,,](3)(3,3)*(2,2)(3,3)*(2,2)[,,]()[,]1,12,(2)(2,2)(3)(1,1)(3,3)*()(1,1)[,,]()()[,]k m k d f x x x c c c x x d d x x d f x x x x x f x x k m k c c c d c x x d f x x x x x x x f x x ==+===-+=-+=⨯-+==+===+=-+=⨯-⨯-+⨯初值为c==,11()()x x f x ⎧⎨-+⎩