[数值算法]求解微分方程的龙格---库塔方法

龙格---库塔方法是求解微分方程比较常用的方法,在理解数学上是怎么一回事后,编制这个程是相当容易的,就是个迭代的过程.步长的选取也是很有讲究的,过小的步长反而会导致误差累积过大.

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

/*标准四阶龙格,库塔公式*/

/*Made by EmilMatthew              

05/8/8*/

#include "RangeKuta.h"

#include "MyAssert.h"

#include "Ulti.h"

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

 

void classicRangeKuta4(Type x0,Type y0,Type xEnd,Type h,

Type (*arguF1)(Type,Type),FILE* outputFile)

{

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

       Type k1,k2,k3,k4;/*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;

             

              k1=(*arguF1)(x0,y0,);

              k2=(*arguF1)(x0+h/2,y0+h*k1/2);

              k3=(*arguF1)(x0+h/2,y0+h*k2/2);

              k4=(*arguF1)(x0+h,y0+h*k3);

      

              y1=y0+h*(k1+2*k2+2*k3+k4)/6;

             

              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);

}

 

/*Test Program*/

#include "MyAssert.h"

#include "RangeKuta.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*************/

              classicRangeKuta4(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;

}