摘要：
简介 维基百科简介：内维尔方法在数学中，内维尔算法是一种计算插值多项式的方法，由数学家埃里克·哈洛德·内维尔提出。从一个给定的n+1节点，有一个幂≤n的唯一多项式，它通过一个给定点。满足以下迭代关系[egin{eqnarray}egin{aligned}&p{i，i}=y_icr&p_{i，j}=frac{p_{i，j-1}+p{i+1，j}}{x_j-x_i}，quad0leilejlenend{aligned}end{eqnarray}]。以n=4的节点为例，迭代过程为[egin{eqnarray}egin{aligned}&p{1,1}=y_1，cr&p{2,2}=y_2，p_{1,2}，cr&p{3,3}=y_3，p_{2,3}，p_{1,3}，cr&p{4,4}=y_4，p_{3,4}，p_{2,4}，p_｛1,4}crend{aligned}end{eqnarray}]代码伪代码，因为计算插值点是矢量，而不是多级循环嵌套，每个元素都被重写为向量形式，每个元素都存储插值点处多项式的函数值。外部循环的数量每次都会减少。此时，它将从最后一行丢弃。只有forirow=1:%，分别被变量x1和x2替换。
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}) 满足以下迭代关系
[egin{eqnarray} egin{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 0le ile j le n end{aligned} end{eqnarray}]
以n=4的节点举例，其迭代过程为
[egin{eqnarray} egin{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