码迷,mamicode.com
首页 > 编程语言 > 详细

数值计算算法-多项式插值算法的实现与分析

时间:2018-05-20 10:44:44      阅读:163      评论:0      收藏:0      [点我收藏+]

标签:之间   管理   初始化   eof   介绍   int   col   有一个   double   

数值计算是指在数值分析领域中的算法。数值分析是专门研究和数字以及近似值相关的数据问题,数值计算在数值分析的研究中发挥了特别重要的作用。

多项式插值是计算函数近似值的一种方法。其中函数值仅在几个点上已知。

该算法的基础是建立级数小于等于n的一个插值多项式pn(z),其中n+1是已知函数值的点的个数。

多项式插值法

许多问题都可以按照函数的方式来描述。然而,通常这个函数是未知的,我们只能通过少量的已知点来推断函数的大致模型。为了实现这个目的,在已知点之间做插值处理。如图所示,关于函数f(x),已知的点为x0...x8,在图中以黑色的圆点表示。通过插值法的帮助,我们能够获取函数在z0、z1、z2处的值,在图中以白色小方块表示。本节主要讨论多项式插值法。

技术分享图片

多项式插值法的根本点就是要建立一个特殊形式的多项式,称为插值多项式

为了深入理解插值多项式的意义,我们先来看看多项式的一些基本法则:

首先,多项式是具有如下形式的函数:

p(x) = a0 + a1x + a2x2 + ... + anxn

这里的a0,...,an是系数当an为非零整数时,这种形式的多项式称为n阶多项式。这是多项式的指数形式,在数学问题中尤为常见。但是,在某些特定的环境中其他形式的多项式则更为简便。比如,在有关多项式插值问题中,牛顿插值多项式就是一个很好的例子

p(x) = a0 + a1 (x - c1) + a2 (x - c1)(x - c2) + ... + an(x - c1)(x - c2)...(x - cn)

这里a0,...,an系数,而c0,...,cn中值。注意到,当c0,...,cn全为0时,牛顿插值多项式就退化为前面定义的n阶多项式。

构建插值多项式

下面我们来看看如何对函数f(x)构建一个插值多项式。

为了对函数f(x)进行插值,要构建一个阶数小于等于n的多项式pn(z),而这又需要用到函数f(x)的n+1个已知点:x0,...,xn这些已知点x0,...,xn就称为插值点。通过插值多项式pn(z)可以计算出函数f(x)在x=z处的近似值。插值法需要满足点z在[x0,xn]内。可以采用如下公式来构建插值多项式pn(z)。

pn(z) = f[x0] + f[x0,x1](z-x0) + f[x0,x1,x2](z-x0)(z-x1) + ... + f[x0,...,xn](z-x0)(z-x1)...(z-xn-1)

其中函数f(x)在点x0,...,xn处的值已知,而f[x0],...,f[x0,...,xn]则称为差商

差商可以通过点x0,...,xn以及函数f(x)在这些点处的值来计算得出。这就是牛顿插值多项式的计算公式。注意到该公式同牛顿插值多项式的相同点。差商的计算公式为:

f[xi,...,xj] = f(xi)                                               如果i=j

f[xi,...,xj] = ( f[xi+1,...xj] - f[xi,...xj-1] ) / (xj - xi) 如果i < j

 通过这个公式不难看出,当i < j时,必须预先计算出其他的差商值。例如,要计算f[x0,x1,x2,x3],就需要先计算出f[x1,x2,x3]和f[x0,x1,x2]的值。幸运的是,可以通过一个差商表来帮助我们以一种系统的方式来计算差商值。如下图。

差商表由多行组成。最顶行保存已知点x0,...,xn 的值。第二行保存f[x0],...,f[xn]的值。要计算出表中其他的差商值,从每个待求的差商值处画一条对角线,使其回到f[xi]和f[xj](如下图中差商f[x1,x2,x3]处的虚线)要得到分母中的xi和xj,直接通过xi和xj求得。分子中的两个差商就是前一阶段计算出来的结果

当完成了整个差商表的计算后,插值多项式的系数就是从第二行开始,每行最左边的那一项

技术分享图片

计算插值多项式

一旦确定了插值多项式的系数,对于函数f,如果我们想知道某个点处的函数值,只需要对多项式求值即可

比如,已知函数f在4个点处的函数值:x0=-3.0,f(x0)=-5.0;x1=-2.0,f(x1)=-1.1;x2=2.0,f(x2)=1.9;x3=3.0,f(x3)=4.8;现在要求出点z0=-2.5,z1=0.0,z2=1.0,z3=2.5处的函数值。由于已经知道函数f在4个点处的值,因此插值多项式为3阶。下图是3阶插值多项式p3(z)的差商表。

技术分享图片

一旦从差商表中得到了系数,就可以采用前面介绍过的牛顿公式来构建插值多项式p3(z)

p3(z)=-5.0 + 3.9(z+3.0) + (-0.63)(z+3.0)(z+2.0) + 0.1767(z+3.0)(z+2.0)(z-2.0)

下一步可以通过该多项式计算出每个点z处的函数值。比如,在点z=-2.5处,经过如下计算得到

p3(z)=-5.0 + 3.9(-2.5+3.0) + (-0.63)(-2.5+3.0)(-2.5+2.0) + 0.1767(-2.5+3.0)(-2.5+2.0)(-2.5-2.0) = -2.694

点z1、z2、z3处的函数值可以通过相似的方法计算得出。最终结果以表格和函数图像的方式表达。如下图。

技术分享图片

和任何其他近似算法一样,通常会有一些与插值多项式相关的误差出现,理解这一点很重要。定性的来讲,如果要使误差降至最小,构建的插值多项式必须要在函数f(x)上获取足够多的已知点才行。并且点与点之前的距离要适当,这样最终得到的多项式才能精确地表示出函数的特性。

多项式插值的接口定义

interpol


int interpol ( const double *x, const double *fx, int n, double *z, double *pz, int m );

返回值:如果插值操作成功,返回0;否则返回-1;

描述:采用多项式插值法来求得函数在某些特定点上的值。

由调用者在参数x处指定函数值已知的点集。每个已知点所对应的函数值都在fx中指定。对应待求的点由参数z来指定,而z所对应的函数值将在pz中返回x和f(x)中的元素个数由参数n来表示。z中的待求点的个数(以及pz中返回值个数)由参数m来表示。由调用者管理x、fx、z以及pz所关联的存储空间。

复杂度:O(mn2),这里m代表待求值的个数,而n代表已知点的个数。

多项式插值的实现与分析

多项式插值法主要基于对一系列期望点的插值多项式的确定。要得到这个多项式,首先必须通过计算差商来确定该多项式的系数。

首先,为差商以及待确定的系数分配存储空间。注意,由于差商表中每一行的每一项都仅依赖于其前一行的计算结果,因此,并不需要一次性保留所有的表项。所以,只为最占用空间的行分配空间即可该行将有n个条目

接下来,用fx中的值来初始化差商表的第一行。这是为计算差商表中的第三行做准备。(前两行不需要计算,因为这两行中的条目都已经保存在x和fx中)。

初始化的最后一步是在coeff[0]中保存fx[0]的值,因为这是插值多项式的第一个系数

计算差商的过程涉及一个嵌套循环,我们在循环中根据前面介绍过的公式来计算差商。在外层循环中,k用来统计正在计算的是哪一行(排除x和fx所代表的行)。在内层循环中i表示在当前行中正在计算的是哪一个条目一旦计算完一行的条目,table[0]中的值就成为插值多项式的下一个系数。因此,保存该值到coeff[k]中一旦得到插值多项式的所有系数,就可以计算出z中每个目标点的值,最后将这些值保存在pz中。

 我们命名这个函数为interpol,它的时间复杂度为O(mn2这里m代表z中的元素个数(也是pz中值的个数),n代表x中(也是fx中)的元素个数。复杂度因子n2是这样得到的,变更j控制循环中的每次迭代,当前迭代中需要乘以的因子比上一轮要多一个。也就是说,当j=1时,term需要做1次乘法,当j=2时,term需要做2次乘法,持续这个过程直到j=n-1时,term需要做n-1次乘法。实际上,这就成了对1~n-1的整数求和,得到的计算时间为T(n)=(n(n+1)/2)-n,再乘以某段固定的时间。(这是由计算等差数列的著名公式得到的)。在大O记法中可以简化为O(n2)。O(mn2)中的因子m来源于针对z中的每个点计算多项式插值的过程。在第一个嵌套循环中,计算出所有的差商,其复杂度为O(n2)。因此,最终的复杂度有一个额外的因子m,该因子对实际的复杂度没有多大的影响。

示例:多项式插值的实现

/*interpol.c*/
#include <stdlib.h>
#include <string.h>

#include "nummeths.h"

/*interpol  */
int interpol(const double *x, const double *fx, int n, double *z, double *pz, int m)
{
    double term, *table, *coeff;
    int i,j,k;
    
    /*为差商和待确定的系数分配空间*/
    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];
    for(k=1; k<n; k++)
    /*外层循环k用来统计正在计算的是哪一行*/
    {
        for(i=0; i<n-k; i++)
        /*内层循环i表示在当前行中正在计算的是哪一个条目(随之行数的增加,条目减少)*/
        {
            j=i+k;
            /*当前行的每一项中分子的差商就是其前一阶段计算的结果*/
            table[i] = (table[i+1] - table[i]) / (x[j] - x[i]);
        }
        /*当前行的首个条目计算结果即是多项式的下一个系数*/
        coeff[k]=table[0];
    }
    free(table);
    
    /*在指定点上对插值多项式进行求值(循环构造插值多项式)*/
    for(k=0; k<m; k++)          /*最外层:遍历z点数组*/
    {
        /*插值多项式的第一个因子*/
        pz[k] = coeff[0];       
        for(j=1; j<n; j++)      /*嵌套:构造多项式(新算式等于上一步的结果加上新因子)*/
        {
            term = coeff[j];    /*因子构成:以多项式系数为基础*/
            for(i=0; i<j; i++)  /*嵌套:新因子以上一步的结果乘以(z[k] - x[i]*/
                term=term*(z[k] - x[i]);
            
            pz[k]=pz[k] + term;
        }
    }
    
    free(coeff);
    return 0;
}

 

数值计算算法-多项式插值算法的实现与分析

标签:之间   管理   初始化   eof   介绍   int   col   有一个   double   

原文地址:https://www.cnblogs.com/idreamo/p/9039000.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!