一个简单的数值运算程序

                 对拉格朗日插值公式做的简单程序

功能：

用拉格朗日插值公式，对给定的n对离散数据进行差值计算。

<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />

 

算法简介：

对给定的n个插值结点x1，x2，……，xn，及其对应的函数值y1=f（x1），y2=f（x2），……，yn=f（xn）；使用拉格朗日插值公式，计算在x点处的对应的函数值f（x）；

                       

程序：

 

 double lagrange( double x0[],double y0[],double n,double x,double *y)

{
 

 int i, j;
 double p;
 *y=0;
 if ( n>1)
 {
  for ( i=0; i<n; i++)
  {
   p=1;
   for ( j=0; j<n; j++)
   {
    if ( i!= j )
     p= p*(x-x0[j]) / ( x0[i]-x0[j]);
   }
   *y = *y + p*y0[i];
   return (0);
  }
 }
 else
  return (-1);
 
}

 

使用说明：

（1） 参数说明

输入参数： x0[ ]  ——  n个元素的一维实数组，存放给定的插值结点x1，x2，……，xn；

           y0[ ] —— n个元素的一维实数组，存放与插值结点相对应的函数值y1，y2，……，yn；

            n —— 整型量，给定插值结点的个数；

            x —— 实型量，插值点。

输出参数： *y —— 实型指针，接受调用程序传送的一个实型量的地址，在程序结束时，在该实型量返回计算结果。

注意： 该实型量中原有内容将被破坏。

（2）调用说明

调用的格式为：

             rtn = lagrange ( x0, y0, n, x, y)；

其中rtn应为一个整型量。

本子程序是一个整型函数，因此在返回主程序一个整型代码于变量rtn中。代码的意义如下：

0 —— 程序正常结束，在y中有计算结果。

-1 —— 程序异常返回，在y中没有结果。异常的原因是n不大于1，使运算无法继续进行。

 

 

例题：

已知函数如下表所示，求X = 0.472处的函数值：

X	0.46	0.47	0.48	0.49
Y	0.484655	0.4903745	0.502750	0.511668

 

计算这个问题的程序如下：

#include "stdio.h"

#include <math.h>

 
void main()
{
 double x0[4]={0.46,0.47,0.48,0.49};
 double y0[4]={0.484655,0.493745,0.502750,0.511668};
 double x,y;
 double n,rtn;
 
 n=4;
 x=0.472;
 
 rtn=lagrange(x0,y0,n,x,&y);
 
 if (rtn==0)
 {
  printf("Y(0.472)=:%f\n",y);
 }
 else
 {
  printf("N must be larger than 1.\n");
 }
}

 

计算结果：Y（0.472）= -0.023263