[数值算法]常微分方程的尤拉方法

[数值算法]常微分方程的尤拉方法:

       尤拉方法是求解常微分方程的入门级的方法，精度并不算高，但它具有较大的理论价值。

一些较好的算法，如龙格.库塔方法等都是在这个方法的基础上实现的。

       (相关的理论请参考相关的数值算法的书籍，我这里只给出关键的函数及主程序段，其余相关的细节就不再一一罗列了.)

       void yulaMethod(Type x0,Type y0,Type xEnd,Type h,Type (*arguF)(Type,Type),FILE* outputFile)

{

       Type x1,y1;/*The node answer,should be print out*/

       Type yp,yc;/*The tween value*/

       int iteratorNum=0;

       int iterLimit=0;

      

       assertF(xEnd-x0>=h,"in classicRangeKuta4  xEnd-x0<h");

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

       assertF(outputFile!=NULL,"in classicRangeKuta4 outputFile is NULL");

      

       iterLimit=(xEnd+0.1-x0)/h;

      

       fprintf(outputFile,"xn:% -12c yn:% -12c /r/n",' ',' ');

       /*Core Program*/

       while(iteratorNum<iterLimit)

       {

              x1=x0+h;

             

              yp=y0+h*(*arguF)(x0,y0);

              yc=y0+h*(*arguF)(x1,yp);

              y1=(yp+yc)/2;

             

              fprintf(outputFile,"% -16f % -16f /r/n",x1,y1);

              x0=x1;

              y0=y1;

 

              iteratorNum++;

       }

 

       /*Output Answer*/

       fprintf(outputFile,"total iterator time is:%d/r/n",iteratorNum);

}

/*tess program*/

#include "MyAssert.h"

#include "Lula.h"

#include <time.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

Type testF1(Type x,Type y);

 

char *inFileName="inputData.txt";

/*

       input data specification

       x0,xEnd means the x value is located about in [x1,x2]

       y0 is the init state of the equation,mapped with x0.

       h is the step or  the adder of x.

*/

 

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 x0,xEnd,y0,h;

       /*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,%f,%f,",&x0,&y0,&xEnd,&h);

       printf("read in data info:%f,%f,%f,%f./n",x0,y0,xEnd,h);

      

#if  DEBUG

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

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

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

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

              yulaMethod(x0,y0,xEnd,h,&testF1,outputFile);

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

       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;

}

 

Type testF1(Type x,Type y)

{

       return y-2*x/y;

}

 

测试结果:

输入：

0            ,1           ,1           ,      0.1

X起始  x终点     y值起点    x的步长.

输出:

xn:             yn:           

0.100000        1.095909       

0.200000        1.184097       

0.300000        1.266201       

0.400000        1.343360       

0.500000        1.416402       

0.600000        1.485955       

0.700000        1.552514       

0.800000        1.616474       

0.900000        1.678166        

1.000000        1.737867      

 

下面是相同输入数据用龙格库塔方法

(http://blog.csdn.net/EmilMatthew/archive/ 2005/08/08 /448496.aspx)所得的结果.

 xn:             yn:           

0.100000        1.095446       

0.200000        1.183217       

0.300000        1.264912       

0.400000        1.341642       

0.500000        1.414215       

0.600000        1.483242       

0.700000        1.549196       

0.800000        1.612455       

0.900000        1.673324       

1.000000        1.732056    

 

可看到,二者的计算精度差大概在10^-3----10^-2级别中,而龙格库塔所得结果则要更接近实际值.