[数值算法]常微分方程的尤拉方法:
尤拉方法是求解常微分方程的入门级的方法,精度并不算高,但它具有较大的理论价值。
一些较好的算法,如龙格.库塔方法等都是在这个方法的基础上实现的。
(相关的理论请参考相关的数值算法的书籍,我这里只给出关键的函数及主程序段,其余相关的细节就不再一一罗列了.)
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级别中,而龙格库塔所得结果则要更接近实际值.