计算数学
牛顿插值
拉格朗日插值多项式每增加一个插值基点，原先计算的插值多项式\({p_n}(x)\) 对 \({p_{n+1}}(x)\) 没有用， 这样必然增加计算工作量，尤其在没有计算机的情况下这个问题更加突出。 我们期望增加插值基点时原先的计算结果对后面的计算仍然有用. 在介绍牛顿插值之前先要了解一下差商的基本概念。 \[\begin{array}{l} f[{x_k}] = f({x_k})\\ f[{x_{k - 1}},{x_k}] = \frac{{f[{x_k}] - f[{x_{k - 1}}]}}{{{x_{k - 1}} - {x_k}}}\\ f[{x_{k - 2}},{x_{k - 1}},{x_k}] = \frac{{f[{x_{k - 1}},{x_k}] - f[{x_{k - 2}},{x_{k - 1}}]}}{{{x_{k - 2}} - {x_k}}}\\ {\rm{ }} \vdots \\ f[{x_{k - j}},{x_{k - j + 1}}, \cdots ,{x_k}] = \frac{{f[{x_{k - j + 1}}, \cdots ,{x_k}] - f[{x_{k - j}}, \cdots ,{x_{k - 1}}]}}{{{x_k} - {x_{k - j}}}} \end{array}\] 有了上面的准备，便可以给出牛顿插值多项式为： \[\begin{array}{l} {N_n}(x) = f[{x_0}] + f[{x_0},{x_1}](x - {x_0}) + f[{x_0},{x_1},{x_2}]\\ {\rm{ }}(x - {x_0})(x - {x_1}) + \cdots + f[{x_0},{x_1}, \cdots ,{x_n}]\\ {\rm{ }}(x - {x_0})(x - {x_1}) \cdots (x - {x_{n - 1}}) \end{array}\] 如果记\[{w_k} = (x - {x_0})(x - {x_1}) \cdots (x - {x_{k - 1}}),k = 1,2, \cdots ,n\] 则牛顿插值可以表达为： \[\begin{array}{l} {N_n}(x) = f[{x_0}] + f[{x_0},{x_1}]{w_1} + f[{x_0},{x_1},{x_2}]{w_2}\\ {\rm{ }} + \cdots + f[{x_0},{x_1}, \cdots ,{x_n}]{w_n}{\rm{ }} \end{array}\]下面给出牛顿插值的C#源代码。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cnewton
/*--------------------------------------class comment
Version : V1.0
Coded by : syz
Date : 2011-07-16 21:33:04 *星期六*
----------------------------------------------------
Desciption : 牛顿插值法
*
parameters :
* fout--------输出文件
*
Methods :
* newton-------牛顿插值法（对多数据一次性完成插值）
* newton0------对单个点插值法
* drive--------驱动程序
--------------------------------------------------*/
{
public StreamWriter fout = new StreamWriter("fout.txt");
public VEC newton(VEC x, VEC y, VEC xx)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-16 21:34:30
---------------------------------------------------
Desciption : 牛顿插值法
*
Post Script :
* 详细算法参看《Fortran科学计算与工程》
paramtegers :
* x--------样本自变量
* y--------样本函数值
* xx-------要插值的向量
* 返回值-------插值结果（维数同xx）
-------------------------------------------------*/
{
int N = x.dim;
int M = xx.dim;
VEC yy=new VEC(M);
for (int i = 0; i < M; i++)
yy[i] = newton0(x, y, xx[i]);
return yy;
}
public double newton0(VEC x, VEC y, double t)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-16 21:36:41
---------------------------------------------------
Desciption : 对单点牛顿插值方法
*
Post Script :
* 编程时需要注意向量与矩阵的指标
paramtegers :
* x,y-------样本自变量与应变量
* t---------要插值的点
-------------------------------------------------*/
{
//获取样本维数
int N = x.dim;
//声明矩阵
MAT Q = new MAT(N, N);
VEC b = new VEC(N);
//把样本函数值复制给Q第一列
for (int i = 0; i < N; i++)
Q[i, 0] = y[i];
for (int i = 1; i < N; i++)
for (int j = 1; j <= i; j++)
Q[i, j] = (Q[i, j - 1] - Q[i - 1, j - 1]) / (x[i] - x[i - j]);
b[N - 1] = Q[N - 1, N - 1];
for (int k = N - 1; k >= 1; k--)
b[k - 1] = Q[k - 1, k - 1] + b[k] * (t - x[k - 1]);
//返回运算结果
return b[0];
}
}