龙格---库塔方法是求解微分方程比较常用的方法,在理解数学上是怎么一回事后,编制这个程是相当容易的,就是个迭代的过程.步长的选取也是很有讲究的,过小的步长反而会导致误差累积过大.
相关的理论请参考相关的数值算法的书籍，我这里只给出关键的函数及主程序段，其余相关的细节就不再一一罗列了.
/*标准四阶龙格,库塔公式*/
/*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;
}