龙格库塔专业论文.docVIP

包头师范学院 2012 -2013学年第一学期期末考试课试卷 考试科目 计算物理 成绩 院系物理学院 专业 10 级本一 姓名张艳有 任课教师签名:院系负责人签名: 任课教师签名: 院系负责人签名: 四阶龙格一库塔算法及C语言程序设计 物理科学与技术学院10级物理系一班1009320036张艳有 摘要木文叙述了四阶龙格一库塔算法的C语言实现过程、数据存储及其结果 并以具体实例说明了这一过程。 关键词 龙格一库塔算法 数据存储例题程序设计计算结果 弓I言 口在科学技术中常常需要求解常微分方程的定解问题，这就需要一种合适的 数值解法求出常微分方程的解。在诸多数值算法屮龙格一库塔算法具有较高的精 确度，是一种优先选取的算法。 具体内容 一、四阶龙格一库塔算法简述 从常微分方程数值解法的几何意义看，欧拉方法取一点匚处的斜率 作为平均斜率，因此欧拉方法近似公式为 Qn+l = 0 + k\t 向后的欧拉方法则采用点0+1处的斜率“2 = f(fn+l^Qn+1)作为平均斜 率，即 Qn+l =Qn+ k? 所以这两种方法也称作矩形法。改进的欧拉方法则取点处和点/曲处斜率 心和灯的平均值作为平均斜率，即 2+1 = Qn +1（*1 +爲）卜 因此改进的欧拉方法乂称为梯形方法。可以预见，若去多点处斜率的加权平均值 作为平均斜率，误差会更小，这就是龙格-库塔方法。 最常用的是四阶龙格库卡近似计算公式，W： 0+i = Qn + + 2k 2 + 2k 3 + 褊)d 6 式中 k\ = f 亿，Qn ) n+—Azk2丿 n+— Az k2 丿 褊=/亿+1?+乂3） 其屮：Ar表示计算过程屮选取的步长； K1表示兀”点处的斜率； K2表示利用K1求得的（乙+空〕点处的斜率; I 2丿 K3表示利用K2求得的1点处的斜率； I 2） K4表示利用K3求得的（￡+）点处的斜率； 二、使用标准四阶龙格——库塔法求解初值问题实例 1、用标准四阶龙格——库塔法求解 yf=2xy yf=2xy(0 y(0)=l 解：取步长*02, J 0.2 2兀+ J 0.2 2兀+ I 2丿 ,J 0.2) “3 = 2 xn + \ /丿 0.2、 2 k2 ( 0.2, 儿 + ? k2 I z 丿 0.2,) 严2 9 (儿+0.2心) Q 2 儿+】=yn + W(“1 + 2*2 + 2*3 + kj 2、计算所用程序如下所示 ttinclude stdio. h〃 #inelude conio. h〃 float func(float x, float y) { return(2*x*y); } float rungc_kutta (float xO, float xn, floeit yO, int n) { ~ float x, y, yl, y2, h, xh; float dl, d2, d3, d4; int i; x=xO; y=yO; h=(xn-x0)/n; for (i=l;i=n;i++) { xh二x+h/2; dl=func (x, y); d2=func (xh,y+h*dl/2.0); d3=func(xh, y+h*d2/2. 0); d4=func(xh, y+h*d3); y二y+h*(dl+2*d2+2*d3+d4)/6. 0; x=x0+i*h; } return (y); } void main() { float xO, xn, yO, e; int n; printf（〃\ninput n:\n〃）； scanf（〃％cT, n）; printf（input xO,xn:\n〃）； scanf（〃％f%f〃, xO, xn）; printf（z，input yO:\n〃）； scanf （〃％f, yO）; e=runge_kutta （xO, xn, yO, n）; printf （z，y （%f） =%6. 6f〃，yO, e）; } 3、部分计算结果如图所示 DOSBox 0.72, Cpu Cycles: max. Frameskip 0f Program: TC 5 input xO^xii: input yO： 0 卩 npiit 5 input 0 xCLxn： input 1 input xO?xn: input input z k2.000000）=5.257334 Input n: yO： 三、 结束语 利用c语言可方便地实现四阶龙格一库塔算法， 这对于科学技术中所遇到的常微分方程的求解具有重要意义。 参考文献： 物理学与计算物理，物理学与计算机翻译组译，北京：科学出版社，1986 计算机数值计算方法，计九三译，北京清华大学出版社，1987