Neville 插值方法

简介

wikipedia: Neville's method

在数学上，Neville 算法是一种计算插值多项式方法，由数学家Eric Harold Neville提出。由给定的n+1个节点，存在一个唯一的幂次≤n的多项式存在，并且通过给定点。

算法

给定n+1个节点及其对应函数值 \((x_i, y_i)\)，假设 \(P_{i,j}\) 表示 \(j-i\) 阶多项式，并且满足通过节点 \((x_k, y_k) \quad k =i, i+1, \cdots, j\)。\(P_{i,j}\) 满足以下迭代关系

\[\begin{eqnarray} \begin{aligned} & p_{i,i}(x) = y_i \cr & P_{i,j}(x) = \frac{(x_j - x)p_{i,j-1}(x) + (x - x_i)p_{i+1,j}(x)}{x_j - x_i}, \quad 0\le i\le j \le n \end{aligned} \end{eqnarray}\]

以n=4的节点举例，其迭代过程为

\[\begin{eqnarray} \begin{aligned} & p_{1,1}(x) = y_1, \cr & p_{2,2}(x) = y_2, p_{1,2}(x), \cr & p_{3,3}(x) = y_3, p_{2,3}(x), p_{1,3}(x),\cr & p_{4,4}(x) = y_4, p_{3,4}(x), p_{2,4}(x), p_{1,4}(x)\cr \end{aligned} \end{eqnarray}\]

代码

伪代码

• 由于计算插值点为一向量，为避免过多层循环嵌套，将每个 \(P_{i,j}\) 都改写为向量形式，各元素分别储存多项式在插值点 \(x_0\) 处函数值。
• 只有每次当一列 \(P_{i,j}\) 计算完后，才能利用迭代公式计算下一列 \(P_{i,j}\) 多项式，因此外层循环为计算每列 \(P_{i,j}\) 多项式。
• 每列 \(P_{i,j}\) 个数是逐渐减少的，最开始有n个多项式，最终循环只有一个。

可将矩阵P[nRow,nCol]用于存储多项式 \(P_{i,j}(x)\)。其中每行为 \(P_{i,j}(x_k)\) 在 nCol 个插值点\(x_k\)处函数值。每次外层循环 \(P_{i,j}(x)\) 个数减少，此时从最后一行开始舍弃，每次只循环

for irow = 1: (nRow - icol) % 

\(x_i\)与\(x_j\)分别用变量x1与x2代替。迭代公式可表示为

for icol = 1:nRow - 1
    for irow = 1: (nRow - icol) % 
        x1 = nodes(irow); x2 = nodes(irow + icol);
        P(irow, :) = ( (x2 - x0).*P(irow, :) + (x0 - x1 ).*P(irow+1, :) )./( x2 - x1 );
    end% for
end% for

最终完整代码为

function evalPol = f1300000_Neville(x0, nodes, fnodes)
% Implement Neville's algorithm to evaluate interpolation polynomial at x0
% Input:
%   x0  - the point where we want to evaluate the polynomial
%   nodes   - vector containing the interpolation nodes
%   fnodes  - vector containing the values of the function
% Output:
%   evalPol - vector containing the value at x0 of the different 
%               the interpolating polynomials

if iscolumn(x0)
    x0 = x0'; % transfer to row vector
end

if isrow(fnodes)
    fnodes = fnodes';
end

nCol = length(x0);
nRow = length(nodes);

% P = zeros(nRow, nCol);
P = repmat(fnodes, 1, nCol);

for icol = 1:nRow - 1
    for irow = 1: (nRow - icol) % 
        x1 = nodes(irow); x2 = nodes(irow + icol);
        P(irow, :) = ( (x2 - x0).*P(irow, :) + (x0 - x1 ).*P(irow+1, :) )./( x2 - x1 );
    end% for
end% for

evalPol = P(1,:);
end