牛顿插值法matlab程序解析.pdf

牛顿插值法在 MATLAB 中的实现 经过 n+1 个不同的插值点12n+1,x xx,，构造牛顿插值公式 1211231212n+112n=,(),()(),()()()Nf x xxxf x xxxxxxf x xxxxxxxx（x）注：牛顿插值法中，用到了插值公式%我们以二次牛顿插值公式为例解析牛顿插值法的 matlab 程序 functionc,d=newpoly(x,y)%这里 x 为 3 个节点的横坐标组成的向量，即123,xx x x，y 为纵坐标的组成向量，即 123,yf xf xf x%c 为所得的牛顿插值多项式的系数组成的向量 n=length(x);%测量向量 x 的长度，即向量 x 中元素ix的个数，赋值给 n，所以 n=3，注：这里的“n”仅为变量，和公式中的次数 n 不一样 d=zeros(n,n);d=zeros(3,3)%把变量 d 定义为一个 n 行，n 列的零矩阵,此矩阵用来储存各阶差商，格式完全等同于书中 21 页的表 d(:,1)=y;%此句是把向量 y 的转置，即123()()()f xyf xf 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)表示1kx，这种计算差商的方法是根据差商表的排列位置而得来，具体解释见下面。end end%下面以二次牛顿插值公式为例解析双循环构造差商表，让我们先来看看构造好的差商表 121232312333()(),(),Xf xdf xf x xf xf x xf 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 xf xddkdf xxxxxxjf xf xddkdf xxxxxxfxf xxddjkdf x xxxxxx，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 和以kx为单根的多项式相乘起来，即()kcxx，k 循环结束后即为11n()()()kcxxxxxx，注初值 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,1 3,(3,3),(3)(3,3)*(2,2)(3,3)*(2,2),(),1,1 2,(2)(2,2)(3)(1,1)(3,3)*()(1,1),()(),km kdf x x xcccx xddx xdf x x xx xf x xkm kcccdcx xdf x x xx xx xf x x 初值为 c=,11()()x xf x