数值计算课程设计--经典四阶龙格-库达法解一阶微分方程组.doc

一、经典四阶龙格-库达法解一阶微分方程组 1．摘要：已知点近似精度求解微分方程组 2. 关键词：龙格-库达法(R-K法) 3. 正文 3.1 算法原理： 从点（t0，y0）开始迭代， yk+1=yk+h(f1+2f2+2f3+f4)/6， 其中： f1=f(tk,yk) f2=f(tk+h/2,yk+hf1/2) f3=f(tk+h/2,yk+hf2/2) f4=f(tk+h,yk+hf3) 3．2 程序流程图： 3.3 程序代码： #include iostream #includeiomanip using namespace std; float G(float pointT,float pointX,float pointY); float F(float pointT,float pointX,float pointY); void function(float * t,float * x,float * y,float h); struct aa { float t; float x; float y; struct aa*next; }; void statement(void) { cout******************endl; cout*计算微分方程组：*endl; cout*dx/dt=x+2y *endl; cout*dy/dt=3x+2y *endl; cout*的龙格-库塔解 *endl; cout******************endl; } float G(float pointT,float pointX,float pointY) { return(pointX+2*pointY); } float F(float pointT,float pointX,float pointY) { return(3*pointX+2*pointY); } void function(float * t,float * x,float * y,float h) { float g[4],f[4]; g[0]=G(*t,*x,*y); f[0]=F(*t,*x,*y); g[1]=G(*t+h/2,*x+h*g[0]/2,*y+h*f[0]/2); f[1]=F(*t+h/2,*x+h*g[0]/2,*y+h*f[0]/2); g[2]=G(*t+h/2,*x+h*g[1]/2,*y+h*f[1]/2); f[2]=F(*t+h/2,*x+h*g[1]/2,*y+h*f[1]/2); g[3]=G(*t+h,*x+h*g[2],*y+h*f[2]); f[3]=F(*t+h,*x+h*g[2],*y+h*f[2]); *t=*t+h; *x=*x+h*(g[0]+2*g[1]+2*g[2]+g[3])/6; *y=*y+h*(f[0]+2*f[1]+2*f[2]+f[3])/6; } void main() { int step=0,count=0; float a=0,b=0,h=0; struct aa * head=NULL,* point=NULL,*nextpoint=NULL; float * pointT=NULL,* pointX=NULL,* pointY=NULL; char flag=y; while (flag==y||flag==Y) { point=(struct aa *)malloc(sizeof(struct aa)); point-next=NULL; head=point; nextpoint=head; statement(); cout请输入微分区间[a,b]端点:; cinab; point-t=a; cout请输入函数在端点a=%d处的初值(x,y):; cinpoint-xpoint-y; cout请输入微分区间的区间段个数step:; cinstep; h=(b-a)/step; pointT=(float *)malloc(sizeof(float)); pointX=(float *)malloc(sizeof(float)); pointY=(float *)malloc(sizeof(float)); * pointT=point-t; * pointX=point-x; * pointY=point-y; for(count=0;countstep;count++) { point=(struct aa *)malloc(sizeof(struct