[数值算法]积分算法之Simpson算法

[数值算法]积分算法之Simpson算法

       Simpson算法来源于Newton-Cotes公式,是一种插值型求积算法.在高阶时数值的稳定性较差,故常用的是比较低较的几种情况.

       算法的具体实现还是相当简单的:(在这里列举三阶和四阶的情况)

       /*Simpson Algorithm

                            Made by EmilMatthew 05/8/12

*/

Type2 simpson3(Type2 buttom,Type2 up,Type2 (*arguF)(Type2))

{

       Type2 x;/*The node answer,should be print out*/

       Type2 ans;/*The ans*/

       int iteratorNum=0;

       int iterLimit=4;

       int k;/*The argu of iter num*/

       int sList[]={1,3,3,1};

 

       assertF(arguF!=NULL,"in simpson3 arguF is NULL");

      

       /*init the data*/

              x=0;

              ans=0;

              k=0;

       /*Core Program*/

       while(iteratorNum<iterLimit)

       {

              x=buttom+k*(up-buttom)/3;

              ans+=sList[k]*(*arguF)(x);         

              k++;

              iteratorNum++;

       }

       ans*=(up-buttom)/8;

      

       return ans;     

}

 

Type2 simpson4(Type2 buttom,Type2 up,Type2 (*arguF)(Type2))

{

       Type2 x;/*The node answer,should be print out*/

       Type2 ans;/*The ans*/

       int iteratorNum=0;

       int iterLimit=5;

       int k;/*The argu of iter num*/

       int sList[]={7,32,12,32,7};

 

       assertF(arguF!=NULL,"in simpson3 arguF is NULL");

      

       /*init the data*/

              x=0;

              ans=0;

              k=0;

       /*Core Program*/

       while(iteratorNum<iterLimit)

       {

              x=buttom+k*(up-buttom)/4;

              ans+=sList[k]*(*arguF)(x);         

              k++;

              iteratorNum++;

       }

       ans*=(up-buttom)/90;

      

       return ans;     

}

 

/*Test Program*/

#include "MyAssert.h"

#include "Simpson.h"

#include <time.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

 

#define p1  0.0306

#define miu 0.0075

#define fi  0.0122

Type2 testF1(Type2 x);

 

char *inFileName="inputData.txt";

/*

       input data specification

       a,b, the limitation of the jifen qujian

*/

 

char *outFileName="outputData.txt";

#define DEBUG 1

 

void main(int argc,char* argv[])

{

       FILE *inputFile;/*input file*/

       FILE *outputFile;/*output file*/

 

       double startTime,endTime,tweenTime;/*time callopsed info*/

       float a,b;/*read in limitation of the jifen qu jian*/

       double tmpAns;

       float i;/*iter num*/

 

       /*input file open*/

       if(argc>1)strcpy(inFileName,argv[1]);

       assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");

       printf("input file open success/n");

      

       /*outpout file open*/

       if(argc>2)strcpy(outFileName,argv[2]);

       assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");

       printf("output file open success/n");

      

       /*Read info data*/

       fscanf(inputFile,"%f,%f,",&a,&b);

       printf("read in data info:%f,%f./n",a,b);

      

#if  DEBUG

       printf("/n*******start of test program******/n");

       printf("now is runnig,please wait.../n");

       startTime=(double)clock()/(double)CLOCKS_PER_SEC;

       fprintf(outputFile,"[");

       /******************Core program code*************/

                            tmpAns=simpson4(-2,1,&testF1);

                            printf("%f",tmpAns);

                            fprintf(outputFile,"%f ",tmpAns);

       /******************End of Core program**********/

       fprintf(outputFile,"];");

       endTime=(double)clock()/(double)CLOCKS_PER_SEC;

       tweenTime=endTime-startTime;/*Get the time collapsed*/

       /*Time collapsed output*/

       printf("the collapsed time in this algorithm implement is:%f/n",tweenTime);

       fprintf(outputFile,"the collapsed time in this algorithm implement is:%f/r/n",tweenTime); 

       printf("/n*******end of test program******/n");

#endif

 

       printf("program end successfully,/n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer./n");

       getchar();/*Screen Delay Control*/

       return;

}

 

#define p1  0.0306

#define miu 0.0075

#define fi  0.0122

 

Type2 testF1(Type2 x)

{

       return 3*x-9*x*x;

}

 

在用一些积分函数进行严格测试后,竟发现Float类型在精度计算时根本毫无招架之力.必用用Double.

<The Art of C Numberic Programing>一书的第一章中提到数值的计算精度一直是程序员和科学工作者之间的鸿沟,今天实践后方有感悟.