您好,欢迎访问三七文档
当前位置:首页 > 商业/管理/HR > 经营企划 > 四阶Runge-Kutta方法
2012-2013(1)专业课程实践论文四阶Runge-Kutta方法李元东,0818180107,R数学08-1班一、算法理论常微分方程初值问题的数值解法是求方程(1)的解在点列1(0,1,)nnnxxhn上的近似值ny,这里nh是1nx到nx的步长,一般略去下标记为h。00(,)()dyfxydxyxy(1)经典的RK方法是一个四阶的方法,它的计算公式是:112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK(2)RK方法的优点是:单步法、精度高,计算过程便于改变步长,缺点是计算量较大,每前进一步需要计算四次函数值f。在用龙格库塔方法时,要注意n的选择要合适,n太大,会使计算量加大,n太小,h较大,可能会使误差增大。因此选择合适的n很重要。我们要在考虑精度的基础上,选择合适的n。在此,用C语言实现了龙格库塔方法。二、算法框图开始输入hbay,,,0x[i+1]=x[i]+hk1=function(x[i],y[i])k2=function(x[i]+h/2,y[i]+h*k1/2)k3=function(x[i]+h/2,y[i]+h*k2/2)k4=function(x[i]+h,y[i]+h*k3)y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/60i输出x[i+1]y[i+1]?ni结束1ii否是四阶Runge-Kutta流程图三、算法程序#includestdlib.h#includestdio.hintRungeKutta(doubley0,doublea,doubleb,intn,double*x,double*y,double(*function)(double,double)){doubleh=(b-a)/n,k1,k2,k3,k4;inti;x[0]=a;y[0]=y0;for(i=0;in;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h/2,y[i]+h*k2/2);k4=function(x[i]+h,y[i]+h*k3);y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}return1;}doublefunction(doublex,doubley){returny-2*x/y;}//例子求y'=y-2*x/y(0x1);y0=1;intmain(){inti;doublex[6],y[6];printf(用四阶龙格-库塔方法\n);RungeKutta(1,1,2,5,x,y,function);for(i=0;i6;i++)printf(x[%d]=%f,y[%d]=%f\n,i,x[i],i,y[i]);return1;}四、算法实现例1.求解方程110,0yxyxy步长2.0h;解:将程序中yxyreturn/*2改为yxreturn1.输入nbay,,,0的值:1;0;1;52.显示输出结果:000000.00x时,000000.10y;200000.01x时,242800.11y;400000.02x时,583636.12y;600000.03x时,044213.23y;800000.04x时,651042.24y;000000.15x时,436502.35y;例2.求解方程,1,10),1/(30yxxyy步长2.0h;解:将程序中yxyreturn/*2改为xyreturn1/*31.输入nbay,,,0的值:1;0;1;52.显示输出结果:000000.00x时,000000.10y;200000.01x时,727548.11y;400000.02x时,742951.22y;600000.03x时,094181.43y;800000.04x时,829211.54y;000000.15x时,996012.75y;
本文标题:四阶Runge-Kutta方法
链接地址:https://www.777doc.com/doc-3864695 .html