龙格-库塔法,求解常微分方程.doc

1、隆格库塔法求解常微分方程 Matlab 隆格库塔法求解常微分方程 摘要 科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本 章将着 重考察的一阶方程的初值问题?虽然求解常微分方程有各种各样的解析方法，但 解析方法只能 用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数 值解法求解? 本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题， 采用了精度较高的 经典的 四阶隆格库塔法，然后通过对实例运用 Matlab 编程进行计算求解，为了体现计算 结果的精确 性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为 比较.通过比较三 种方法的计算精度，发

2、现四阶经典龙格 -库塔方法的误差最小，预估较 正法其次，欧拉方法 误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格 库塔法求解常微分方程初 值问题的计算精度的影响? 总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解 法， 隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度 也最高. 关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程; Runge Kutta Method For Solving Ordinary Differential Equations ABSTRACT Problem sol ving ordi nary

3、 differe ntial equati ons are ofte n n eeded in scie nee an dtech no logy, the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various an alytical methods for sol ving ordinary differential equations, the

4、 analytical method can only be used to solve some special types of equati on s.differe ntial equati ons can be summed up the actual problems which This paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with h

5、igh accuracy.for instanee through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correct ion method for in sta nee by calculati on through the calculati on precisi on. The comparis on of thre

6、e kinds of methods, found that the error of four order Runge Kutta method of mi nimum, predictio n correcti on method sec on dly, Euler method error is relatively large. Fin ally, by selecti ng differe nt step, study the affect the calculatio n accuracy of differe nt step of Runge Kutta method to so

7、lve initial value problems of ordinary differential equations. In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergenee, the Run

8、ge Kutta method of order four of the highest efficiency and its precision is the highest. Key words：Four order Runge Kutta method; Euler method; predicti on correction method; first order ordinary differential equations; Matlab 隆格库塔法求解常微分方程 1 问题的提出.1 1.1问题背景.1 1.2问题的具体内容.1 2 问题假设.2 3 符号系统.2 4 问题的分析.

9、3 5.1欧拉格式.3 5.2预估较正法.3 5.3四阶隆格库塔法的格式 .4 5 模型的建立与求解.4 5.3隆格库塔法的基本原理.4 5.1.1 Taylor 级数.4 5.1.2隆格库塔法的基本思想.4 5.1.3四阶的隆格库塔法.5 5.2 其他求解常微分方程边值问题算法的简介 .6 5.2.3 模型求解.8 运用 MATLA 软件对模型求解结果及析.8 6 模型的评价.16 7 课程设计的总结与体会 .16 参考文献.17 附录.18 一、问题的提出 1.1 问题背景： 科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过 于一阶常微分方程的初值问题，即： 学=f

10、 (x, y) a 辽 X 岂 b dx .y(a) = y?( 1) 其中 a，b为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于 求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一 阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重 讨论一类高精度的单步法一一隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问 题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的 Euler 格式求解出的结果进行精度比较. 1.2 问题的具体内容 =_y+ x + 1 实例一：在区间0,1上采用经典的四阶隆格库塔方法

11、求解微分方程 dx ，其 汁(0)=1 精确解为 x*e，步长为 o.5,然后用欧拉法，预估校正法分别求解，且将计算结果 与精确解进 行比较，对三个算法的收敛性的进行分析比较 虬 x_y 实例二：在区间0,1上用经典的四阶龙格库塔方法求解初值问题dx e y，其 y(0) = 1 (x2+2)e 精确解为 2 ，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进 行比较，对 三个算法的收 敛性的进行 分析比较.最后在区间0,1上分别取 步长 h=0.1;0.05 隆格库塔法求解常微分方程 第 页/共 22 侦 时进行计算，并且探究选取不同的步长对计算结果精度的影响. 问题假设 2.1 假

12、设数值方法本身的计算是准确的. 2.2 假设选取的步长趋于 0 时计算的结果会收敛到微分方程的准确解 2.3 假设步长的增加不会导致舍入误差的严重积累. 三、符号系统 3.1 符号说明 符号含义 h 选取的步长 * K 平均斜率 p精度的阶数 前后两次计算结果的偏差 yn第 n 个节点的实验值 y(Xn)第 n 个节点的精确值 0实验值与精确值的绝对误差 隆格库塔法求解常微分方程 第页/共223页 ynynhf(Xn, yn) 局部截断误差 h2h2 （2）隐式欧拉格式 ynynhf(xn1,yn1) 局部截断误差 y(Xn 1)- yn 1 : 4.2 预估校正法 预估 yn 1=ynhf

13、(Xn,yn) 校正 h 讣尹（XnS川） 统一格式 八* f(x,y) f( x hnyMU) (2) (3) (4) (8) 平均化格式 (9) 四、问题的分析 问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实 例，本 文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优 越性，同时 用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个 算法的收敛性的 进行分析比较.最后在区间0,1上分别取步长 h=0.1;0.05 时进行计算, 分析在选取不同的 步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及 经典的四阶隆格库塔法的计算

14、公式. 4.1 欧拉格式 （1）显式欧拉格式 yp = yn+ h(f X y, yc = % h(f x,y, y?2(yp 十 y). L.2 p 前2刃?，才隆格库塔法求解常微分方程 (14) 第页/共225页 (10) y(x1戶y( x ) h y(n y(nX)；.(11) 4.3 四阶龙格库塔方法的格式(经典格式) h yn=yn十6(心 +2K2+2K3+K4), Ki= f (Xn, yn), h h K2 = f(Xn+;,yn+；KJ, 2 2 hh K3 = f(Xn+ 二，yn+ 二心)， iK4=f(Xn+h,yn+hK3). 五、模型的建立与模型求解 5.1 隆格

15、库塔法的基本原理 隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联 系. Taylor 级数 ；y = f (x, y) 设初值冋题y(x0) = y0有解，按泰勒展开，有 其中y(x)的各阶导数依据所给方程可以用函数 f 来表达，下面引进函数序列 f j(x,y) 来描叙求导过程，即 其中一阶 Taylor 格式为: yn 1 = ynhyn 汗y (0 ) ( 1 )三 f X: X (1) :f (1) ff (1) ( y 讥迸 ex (j -2 )j 2 ) -三 f j(一1 ) (13) 根据上式，结果导出下面 Taylor 格式h2 yn产 ynh

16、yn刁 yn hp Pp() 2处几扌隆格库塔法求解常微分方程 (19) 第页/共227页 提高 Taylor 格式的阶数 p 即可提高计算结果的精度，显然, p 阶 Taylor 格式的局部截断误差为: y(n .1)_yn 1= hP 1 P 1 (P 1)!y (15) 因此它具有 p 阶精度. 隆格库塔方法的基本思想 隆格库塔法实质就是间接地使用 根据微分中值定理，这有 TaylOr级数法的一种方法,考察差商 y(xn1)hy(xn) y(Xn 1) - y(Xn) h =y (xndh) (16) (17) 利用所给方程y = f (x, y) y(xn 1) = y(Xn) hf

17、(Xnrh,y(XnTh) 设平均斜率二 f(Xn汕,y(XnTh)，由此可见，只要对平均斜率一种算法，便相应 地可以导 出一种计算格式. 再考察改进的 Euler 格式，它可以改写成平均化的形式： * K1= f(Xn,yn)( 18) 这个处理过程启示我们，如果设法在(Xn,Xm)内多预测几个点的斜率值，然后将它们加 权平均 作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基 本思想. 5.1.3 四阶的隆格库塔法 为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如 下： Yn+ = y-(K2K2KK4), 6 心=f (Xn, Yn), Ih

18、h 心=心七”七心)， 22 心(心,*呼2), 22 iK4= f(Xn+h,yn+h y1=； h=0.1; x=0:h:1; n=len gth(x); for i=1: n y(i)=f1(x(i); end figure(1) plot(x,y,g-); hold on y1(1)=1; for j=2: n y1(j)=y1(j-1)+h*f(x(j-1),y1(j-1); end Y=x;y1; fid=fope n(data.txt,wt); fprin tf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2)

19、 DT=abs(y-y1); 隆格库塔法求解常微分方程 第页/共 221 页 plot(x,DT) %1.建立导数函数文件 fun ctio n z=f(x,y) z=y-2*x/y; %2.建立原函数文件 function z1=f1(x) z1= (2*x+1)A(1/2); 迭代 n 次的后退的欧拉格式 matlab 程序 %1.建立导数函数文件 fun ctio n z=f(x,y) z=y-2*x/y; %2.建立原函数文件 function z1=f1(x) z1= (2*x+1)A(1/2); %迭代 n 次的后退的欧拉格式 clear all clc N=input(请输入迭代

20、次数:); x=； y=; y1=; h=0.1; x=0:0.1:1; n=len gth(x); for i=1: n y(i)=f1(x(i)； end figure(1) plot(x,y,g-); hold on 隆格库塔法求解常微分方程 第 页/共 222 页 y1(1)=1; for j=2: n y1(j)=y1(j-1)+h*f(x(j-1),y1(j-1); T=y1(j); for k=1:N T=y1(j-1)+h*f(x(j),T); end y1(j)=T; end Y=x;y1; fid=fope n(data.txt,wt); fprin tf(fid,%6.2

21、f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2) DT=abs(y-y1); plot(x,DT) 附录二：预估校正法 matlab 程序 %1.建立导数函数文件 fun ctio n z=f(x,y) z=y-2*x/y; %2 建立原函数文件 function z1=f1(x) z1= (2*x+1)A(1/2); %预估校正法 %拉法 clear all clc x=; y=; y1=; y1(1)=1; y2=; y2(1)=1; h=0.1; x=0:0.1:1; n=len gth(x); for i=1: n y(i)=f1(

22、x(i); end figure(1) plot(x,y,g-); hold on for j=2: n y1(j)=y2(j-1)+h*f(x(j-1),y2(j-1); y2(j)=y2(j-1)+(h/2).*f(x(j-1),y2(j-1)+f(x(j),y1(j); end Y=x;y2; fid=fope n(data.txt,wt); fprin tf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y2,r-); figure(2) DT=y-y2; plot(x,DT) 附录三：四阶龙格库塔法 matlab 程序 %四阶龙格库塔法 clea

23、r all clc h=0.2; 隆格库塔法求解常微分方程 第 页/共 2 虫页 x=0:h:1; n=len gth(x); y=; y1=; y1(1)=1; for i=1: n y(i)=f1(x(i); end figure(1) plot(x,y,g-); hold on for j=2: n K1= f(x(j-1),y1(j-1); K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1); K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2); K4=f(x(j-1)+h,y1(j-1)+h*K3); y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); end Y=x;y1; fid=fope n( data1.txt,wt); fprin tf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2) DT=abs(y-y1); plot(x,DT)