21、多项式插值，最小二乘法和牛顿迭代法

1、多项式插值

利用函数f (x)在某区间中已知的若干点的函数值，作出适当的特定函数，在区间的其他点上用这特定函数的值作为函数f (x)的近似值，这种方法称为插值法。如果这特定函数是多项式，就称它为插值多项式。

//*x、*fx分别为已知点数组（x0，f（x0））..（xn，f（xn））；n为已知点个数
// *z、*fz分别为待求点数组（z0，f（z0））..（zn，f（zn））；m为待求点个数
int interpol(const double *x, const double *fx, int n, double *z, double *pz, int m)
{
double term, *table, *coeff;
int i, j ,k;

//table存临时差商，每组从一阶到n阶更新，对应上图的每行更新 ，
//如 [x0 x1]，[x1x2],...，[xn-1 xn]；更新到[x0 x1 x2],...,[xn-2, xn-1, xn]
//coeff存最终差商 ，机每次取table更新的第一项，如[x0 x1]到[x0 x1...xn]
if((table = (double *)malloc(sizeof(double) * n)) == NULL)
return -1;
if((coeff = (double *)malloc(sizeof(double) * n)) == NULL)
{
free(table);
return -1;
}

memcpy(table, fx, sizeof(double) * n);
coeff[0] = table[0];

//k对应行
for(k = 1; k < n; k++)
{
//第k行只有n-k个差商
for(i = 0;i < n - k; i++)
{
table[i] = (table[i + 1] - table[i]) / (x[i + k] -x[i])
}

//我们要的差商只是每行第一个即[x0 x1],...,[x0,x1,...,xn]
coeff[k] = table[0];
}
free(table);

//求m个未知点
for(k = 0;k < m; k++)
{
pz[k] = coeff[0];

//j对应公式中每个加数的长度。
for(j = 1; j < n; j++)
{
term = coeff[j];

//下面term为第j个加数对应的an(x-x0)(x-x1)...(x-xn-1) ,需要累乘
for(i = 0; i < j; i++)
term = term * (z[k] - x[i]);

//把这些加数都加起来
pz[k] = pz[k] + term;
}
}
free(coeff);
return 0;
}

2、最小二乘估计法

最小二乘法（又称最小平方法）通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。就是求一条线到所有点距离（x相同y的差值）都近。

//*x，*y对应已知的坐标数组；n为已知点个数；a、b为待求方程y = a + bx的系数
void lsqe(const double *x, const double *y, int n, double n, double *b, double *a)
{
//sumx2位xi的平方和
double sumx, sumy, sumx2, sumxy;
int i;
sumx = 0.0;
sumy = 0.0;
sumxy = 0.0;
sumx2 = 0.0;

for(i = 0; i < n; i++)
{
sumx = sumx + x[i];
sumy = sumy + y[i];
sumx2 = sumx2 + pow(x[i],2);
sumxy = sumxy +(x[i] * y[i]);
}

*b = (sumxy -(sumx * sumy) / (double)n)  / (sumx2 - (pow(sumx, 2.0) / (double)n));
*a = (sumy - (*b) * sumx) / (double)n;
return;
}

3、牛顿迭代法

1、牛顿迭代法通过一系列操作使得结果不断逼近方程的实根。



2、1阶导数表示单调性，2阶导数表示凹凸性。x0的选择要求：（1）x0在[a b]内只有一个跟；（2）x0=a或者x0=b时f（x0）和二阶导数 在[a b]上符号相同。这样保证每次迭代逐步逼近。

//g（x）为f（x）的导数，x为数组存放每次迭代的数，n为最大迭代次数或最终迭代次数，delta为允许误差
int root(double (*f)(double x), double (*g)(double *x), double *x, int *n, double delta)
{
int satisfied, i;
i = 0;
satisfied = 0;

while(!satisfied && i+1 < *n)
{
x[i+1] = x[i] - (f(x[i]) / g(x[i]));
if(fabs(x[i+1] - x[i]) < delta)
satisfied = 1;
i++;
}

//如果没有达到最大迭代次数就符合要求，使n等于此时的迭代次数
if(i == 0)
*n = 1;
else
*n = i + 1;

if(satisfied)
return 0;
else
return -1;
}