对于一列 x 0 , x 1 , … , x n x_0,x_1,\ldots,x_n x0​,x1​,…,xn​，若与 y 0 , y 1 , … , y n y_0,y_1,\ldots,y_n y0​,y1​,…,yn​通过映射 f f f一一对应，则定义比值

f [ x 0 , x 1 ] = f ( x 1 ) − f ( x 0 ) x 1 − x 0 f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0} f[x0​,x1​]=x1​−x0​f(x1​)−f(x0​)​

为 f ( x ) f(x) f(x)关于节点 x 0 , x 1 x_0,x_1 x0​,x1​的一阶差商。

类似地，

f [ x 0 , x 1 , x 2 ] = f [ x 0 , x 2 ] − f [ x 0 , x 1 ] x 2 − x 1 f[x_0,x_1,x_2]=\frac{f[x_0,x_2]-f[x_0,x_1]}{x_2-x_1} f[x0​,x1​,x2​]=x2​−x1​f[x0​,x2​]−f[x0​,x1​]​

为二阶差商，

f [ x 0 , x 1 , … , x k ] = f [ x 0 , x 1 , … , x k − 2 , x k ] − f [ x 0 , x 1 , … , x k − 1 ] x k − x k − 1 f[x_0,x_1,\ldots,x_k]=\frac{f[x_0,x_1,\ldots,x_{k-2},x_k]-f[x_0,x_1,\ldots,x_{k-1}]}{x_k-x_{k-1}} f[x0​,x1​,…,xk​]=xk​−xk−1​f[x0​,x1​,…,xk−2​,xk​]−f[x0​,x1​,…,xk−1​]​

为 f ( x ) f(x) f(x)关于节点 x 0 , x 1 , … , x k x_0,x_1,\ldots,x_k x0​,x1​,…,xk​的 k k k阶差商。

k k k阶差商可表示为函数值 f ( x 0 ) , f ( x 1 ) , … , f ( x n ) f(x_0),f(x_1),\ldots,f(x_n) f(x0​),f(x1​),…,f(xn​)的线性组合，即

f [ x 0 , x 1 , … , x k ] = ∑ j = 0 k f ( x j ) ( x j − x 0 ) … ( x j − x j − 1 ) ( x j − x j + 1 ) … ( x j − x k ) f[x_0,x_1,\ldots,x_k]=\displaystyle\sum_{j=0}^k\frac{f(x_j)}{(x_j-x_0)\ldots(x_j-x_{j-1})(x_j-x_{j+1})\ldots(x_j-x_k)} f[x0​,x1​,…,xk​]=j=0∑k​(xj​−x0​)…(xj​−xj−1​)(xj​−xj+1​)…(xj​−xk​)f(xj​)​

则

f [ x 0 , x 1 , … , x k ] = f [ x 1 , x 2 , … , x k ] − f [ x 0 , x 1 , … , x k − 1 ] x k − x 0 f[x_0,x_1,\ldots,x_k]=\frac{f[x_1,x_2,\ldots,x_k]-f[x_0,x_1,\ldots,x_{k-1}]}{x_k-x_0} f[x0​,x1​,…,xk​]=xk​−x0​f[x1​,x2​,…,xk​]−f[x0​,x1​,…,xk−1​]​

则可通过递归实现两个数组的差商

根据差商定义，可得

f ( x ) = f ( x 0 ) + f [ x , x 0 ] ( x − x 0 ) f [ x , x 0 ] = f [ x 0 , x 1 ] + f [ x , x 0 , x 1 ] ( x − x 1 ) ⋮ f [ x , x 0 , … , x n − 1 ] = f [ x 0 , x 1 , … , x n ] + f [ x , x 0 , x 1 , … , x n ] ( x − x n ) \begin{aligned} f(x) &= f(x_0) + f[x,x_0](x-x_0)\\ f[x,x_0] &= f[x_0,x_1] + f[x,x_0,x_1](x-x_1)\\ &\vdots\\ f[x,x_0,\ldots,x_{n-1}] &= f[x_0,x_1,\ldots,x_n] + f[x,x_0,x_1,\ldots,x_n](x-x_n) \end{aligned} f(x)f[x,x0​]f[x,x0​,…,xn−1​]​=f(x0​)+f[x,x0​](x−x0​)=f[x0​,x1​]+f[x,x0​,x1​](x−x1​)⋮=f[x0​,x1​,…,xn​]+f[x,x0​,x1​,…,xn​](x−xn​)​

即

f ( 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 ) + … + f [ x 0 , x 1 , … , x n ] ( x − x 0 ) ( x − x 1 ) … ( x − x n − 1 ) + f [ x 0 , x 1 , … , x n ] ω n + 1 ( x ) = N n ( x ) + R n ( x ) \begin{aligned} f(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)\\ &+\ldots+f[x_0,x_1,\ldots,x_n](x-x_0)(x-x_1)\ldots(x-x_{n-1})\\ &+f[x_0,x_1,\ldots,x_n]\omega_{n+1}(x)\\ =&N_n(x)+R_n(x) \end{aligned} f(x)==​f(x0​)+f[x0​,x1​](x−x0​)+f[x0​,x1​,x2​](x−x0​)(x−x1​)+…+f[x0​,x1​,…,xn​](x−x0​)(x−x1​)…(x−xn−1​)+f[x0​,x1​,…,xn​]ωn+1​(x)Nn​(x)+Rn​(x)​

其中，

N 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 ) + … + f [ x 0 , x 1 , … , x n ] ( x − x 0 ) ( x − x 1 ) … ( x − x n − 1 ) R n ( x ) = f [ x 0 , x 1 , … , x n ] ω n + 1 ( x ) \begin{aligned} N_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)+\\ &\ldots+f[x_0,x_1,\ldots,x_n](x-x_0)(x-x_1)\ldots(x-x_{n-1})\\ R_n(x)=&f[x_0,x_1,\ldots,x_n]\omega_{n+1}(x) \end{aligned} Nn​(x)=Rn​(x)=​f(x0​)+f[x0​,x1​](x−x0​)+f[x0​,x1​,x2​](x−x0​)(x−x1​)+…+f[x0​,x1​,…,xn​](x−x0​)(x−x1​)…(x−xn−1​)f[x0​,x1​,…,xn​]ωn+1​(x)​

N n ( x ) N_n(x) Nn​(x)即为牛顿插值公式，且

N n + 1 ( x ) = N n ( x ) + f [ x 0 , x 1 , … , x n + 1 ] ω n + 1 ( x ) N_{n+1}(x) = N_{n}(x) + f[x_0,x_1,\ldots,x_{n+1}]\omega_{n+1}(x) Nn+1​(x)=Nn​(x)+f[x0​,x1​,…,xn+1​]ωn+1​(x)

可快速写为Julia

当然，无论差商还是连乘，运算过程中都会产生大量重复，上述代码并未有效利用已经计算出的值，所以性能并不占优。因此，对Newton插值进行更改

调用并验证