1、隆格库塔法求解常微分方程 隆格库塔法求解常微分方程 摘要 科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本 章将着重考察的一阶方程的初值问题. 虽然 求解常微分方程有各种各样的解析方法，但 解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数 值解法求解 . 本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的 经典的四阶隆格库塔法 , 然后通过对实例运用Matlab 编程进行计算求解， 为了体现计算 结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为 比较. 通过比较三种方法的计算精度，发现四阶经典龙格
2、-库塔方法的误差最小，预估较 正法其次，欧拉方法误差则比较大 .最后通过选取不同的步长， 研究了不同的步长对隆格 库塔法求解常微分方程初值问题的计算精度的影响. 总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解 法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度 也最高. 关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；Matlab 隆格库塔法求解常微分方程 Runge Kutta Method For Solving Ordinary Differential Equations ABSTRACT Problem solving
3、ordinary differential equations are often needed in science andtechnology. 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 analytical methods for solving ordinary differential equations, the a
4、nalytical method can only be used to solve some special types of equations.differential equations 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 high ac
5、curacy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of me
6、thods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value probl
7、ems 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 convergence, the Runge Kutta method of orde
8、r four of the highest efficiency and its precision is the highest. Key words:Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab 隆格库塔法求解常微分方程 目录 1 问题的提出.1 1.1问题背景.1 1.2问题的具体内容.1 2 问题假设. 2 3 符号系统. 2 4 问题的分析.3 4.1欧拉格式.3 4.2预估较正
9、法.3 4.3四阶隆格库塔法的格式.4 5 模型的建立与求解.4 5.1隆格库塔法的基本原理.4 5.1.1 Taylor级数.4 5.1.2隆格库塔法的基本思想.4 5.1.3四阶的隆格库塔法.5 5.2 其他求解常微分方程边值问题算法的简介.6 5.3模型求解.8 5.3.1运用MATLA软B 件对模型求解结果及析. 8 6 模型的评价.16 7 课程设计的总结与体会.16 参考文献.17 附录.18 第页/ 共 221页 隆格库塔法求解常微分方程 1.1 问题背景： 一、问题的提出 科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过 于一阶常微分方程的初值问题，即：
10、 dy f (x, y) dx y(a)y0 axb （1） 其中 a，b 为常数 . 虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于 求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解. 因为一 阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法. 本文着重 讨论一类高精度的单步法隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问 题， 并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较. 1.2 问题的具体内容 实例一：在区间0,1上采用经典的四阶隆格库塔方法求解微分方程 dy yx1
11、 dx ，其 y(0)1 精确解为 yxe，步长为 0.5, 然后用欧拉法，预估校正法分别求解，且将计算结果 与精确解进行比较，对三个算法的收敛性的进行分析比较. 实例二：在区间 0,1 上用经典的四阶龙格库塔方法求解初值问题 dy xe x y dx， 其 y(0)1 精确解为 1 (x 2 2 2)e x ，然后用欧拉法， 预估校正法分别求解， 且将计算结果与精确解进 行比较，对三个 算法的收敛性 的进行 分 析比 较. 最后 在区间 0,1上分别取 步长 h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响. x 第页/ 共 222页 隆格库塔法求解常微分方程 二、
12、问题假设 2.1 假设数值方法本身的计算是准确的. 2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解. 2.3 假设步长的增加不会导致舍入误差的严重积累. 三、符号系统 3.1 符号说明 符号 h 含义 选取的步长 K * p 平均斜率 精度的阶数 前后两次计算结果的偏差 yn 第 n 个节点的实验值 y (xn) 第 n 个节点的精确值 实验值与精确值的绝对误差 第页/ 共 223页 隆格库塔法求解常微分方程 四、问题的分析 问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实 例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优 越
13、性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个 算法的收敛性的进行分析比较. 最后在区间 0,1 上分别取步长 h=0.1;0.05时进行计算， 分析在选取不同的步长时，求解结果的精度变化如何. 下面是欧拉法，预估校正法以及经 典的四阶隆格库塔法的计算公式. 4.1 欧拉格式 （1）显式欧拉格式 局部截断误差 yn 1yn h 2 hf (xn, yn) h 2 （2） y(x)yy ( )y (x )o(h 2 )（3） （2）隐式欧拉格式 n 1n 1n 22 局部截断误差 yn 1ynhf (xn h 2 1, yn 1) （4） y(x)yy (x )o(h
14、 2 )(5) 4.2 预估校正法 预估 n 1n 1 yn 1 2 ynhf n (xn, yn)(6) 校正yy h f (x , y ) f (x, y) (7) n 1nnnn 1n 1 2 统一格式 yy h f( x,y)f( x,h yh(f,x )y)(8) n 1nnnnnnn 2 ypynh(f nx ,ny ), 平均化格式 ycynh(f nx1,yp),(9) yn 1 1 ( ypyc). 2 第页/ 共 224页 n ,y 2 ) 2p 4.3 四阶龙格库塔方法的格式（经典格式） h 隆格库塔法求解常微分方程 yn 1yn(K1 6 2K22K3K4), K1f(
15、xn,yn), K2f(xn h ,y h 22 K1),(10) K3f hh (xn nK ), 22 K4f(xnh, ynhK3). 5.1 隆格库塔法的基本原理 五、模型的建立与模型求解 隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联 系. 5.1.1 Taylor级数 yf (x, y) 设初值问题 y(x0)y0 有解，按泰勒展开，有 y( x)y( x ) h y( x h y(x)；. . . .（11） n 1nnn 2 其中y( x)的各阶导数依据所给方程可以用函数f 来表达， 下面引进函数序列f j (x, y) 来描叙求导过程，即 f
16、( 0 ) f ( 1 ) yff ( 0,)y ff xx ( 1 ) （12） f ( 1 ) f ( 1 ) ( 2 ) yff xy f ( j2 ) f j (2 ) y ( j ) ff j(1 ) xy （13） 根据上式，结果导出下面Taylor 格式 yyhy h y . h y p () （14） n 1nnnn 2!p! 其中一阶 Taylor 格式为：yyhy n 1nn 2 第页/ 共 225页 ,y 1 2 隆格库塔法求解常微分方程 提高 Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶 Taylor 格式的局部截断误差为： h p 1 yyy p
17、1 （ 15） 因此它具有p 阶精度 . 5.1.2隆格库塔方法的基本思想 (n 1)n 1 ( p1)! 隆格库塔法实质就是间接地使用Taylor 级数法的一种方法 , 考察差商 根据微分中值定理，这有 y(xn 1) h y( xn) y(xn 1) h y( xn) y (xnh) （16） 利用所给方程 yf (x, y) y(xn 1)y(xn)hf (xnh, y(xnh)（17） 设 平均斜率 K*f (xnh, y(xnh) ，由此可见，只要对平均斜率一种算法，便相应 地可以导出一种计算格式 . 再考察改进的 Euler 格式，它可以改写成平均化的形式： yy h (KK )
18、n 1n 2 12 K1f(xn, yn) （18） K2f(xn 1, ynhK1) 这个处理过程启示我们，如果设法在(xn,xn1) 内多预测几个点的斜率值，然后将它们加 权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基 本思想. 5.1.3四阶的隆格库塔法 为了方便起见， 本文主要运用经典的隆格库塔算法- 四阶隆格库塔格式 . 其格式如 下： yy h (K2K 2KK ), n 1n 6 1234 K1f(xn, yn), K2f hh (xn nK ), （19） 22 K 3 f (xn h , y 2 n h K ), 2 K4f(xnh, ynhK3
19、). 下面为其具体的算法流程图： 第页/ 共 226页 隆格库塔法求解常微分方程 开始 输入 x0,y0,h,N n=1 x1=x0+h k1=f(x0,y0),k2=f(x0+h/2,y0+hk1/2) k3=f(x0+h/2,y0+hk2/2),k4=f(x1,y0+hk3) y1=y0+h(k1+2k2+2k3+k4)/6 输出 x1,y1 n=N? 是 结束 否 n=n+1 x0=x1,y0=y1 图 5.1龙格 -库塔法流程图 5.2 其他求解常微分边值问题算法的简介 5.2.1 欧拉数值算法（显式） 微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即： dy f dx y
20、(a) (x, y) y0 axb (20) 其中 a，b为常数. 因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法. 所有算法中的 f 第页/ 共 227页 n 就是代表上式中 f ( x,y)，而fy表示 f ( x, y) y ， fx表示 f (x, y) . x 隆格库塔法求解常微分方程 简单欧拉法是一种单步递推算法. 简单欧拉法的公式如下所示： yn 1ynhf(xn, yn)(21) 简单欧拉法的算法过程介绍如下: 给出自变量 x的定义域a,b，初始值y0及步长h. 对 k0,1,( ba) / h，计算yk 1ykhf(xk, yk) 5.2.2欧拉数值算法（隐式）
21、 隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示： yn 1ynhf (xn 1, yn 1)(22) 隐式欧拉法是一阶精度的方法，比它精度高的公式是： yy h f (x , y ) f (x, y) (23) n 1nnnn 1n 1 2 隐式欧拉的算法过程介绍如下. 给出自变量 x的定义域a,b，初始值y0及步长h. 对k0,1,( ba) / h，用牛顿法或其他方法求解方yk 1ykhf (xk 1, yk 1) 得出 yk 1. 5.2.3欧拉预估 - 校正法 改进欧拉法是一种二阶显式求解法，其计算公式如下所示： hh tyn hf (xn, yn) yn 1ynhf xn 2 ,
22、 yn 2 f (xn, yn) yy h f (x , y )f (x,t) (24) n 1nnnn 1 2 四阶龙格 - 库塔法有多种形式，除了改进的欧拉法外还有中点法. 中点法计算公式为： yn 1 ynhf xn h , y h 22 f (xn,yn) (25) 第页/ 共 228页 初 值 Euler 法 相对误 差 预估校 正法 相对误 差 经典四 阶库 相对误差精确值 0-1.00000 0.11.00910 0.00424 1.00500 0.00016 1.004840.000001.00484 0.21.02646 0.00759 1.01903 0.00029 1.0
23、18730.000001.01873 0.31.05134 0.01011 1.04122 0.00038 1.040820.000001.04082 0.41.08304 0.01189 1.07080 0.00045 1.070320.000001.07032 0.51.12095 0.01303 1.10708 0.00049 1.106530.000001.10653 0.61.16451 0.01366 1.14940 0.00052 1.148810.000001.14881 0.71.21319 0.01388 1.19721 0.00052 1.196590.000001.19
24、659 0.81.26655 0.01378 1.24998 0.00052 1.249330.000001.24933 0.91.32414 0.01344 1.30723 0.00050 1.306570.000001.30657 1.01.38558 0.01294 1.36854 0.00048 1.367880.000001.36788 隆格库塔法求解常微分方程 5.3 模型求解 5.3.1 运用 MATLAB 软件对模型求解结果及分析 用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题， 并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：
25、实例一. dy dx yx1 y(0)1 由题可知精确解为：yxe x ，当 x=0 时， y(x)=0.在这里取步长h 为 0.1 ， 通过 MATLAB 程序的计算，相应的结果如下： 表 5-1 步长为 0.1 时各方法的预测值与精确值的比较(精确到 6 位小数 ) 第页/ 共 229页 Y 1.4 1.3 1.2 步长为 0.1 时的精确值与 预 测 值 的 比 较 精确值 欧拉法 隆格库塔法求解常微分方程 1.1 1 0.9 0.8 轴 0.7 Y 0.6 0.5 0.4 0.3 0.2 0.1 0 改进欧拉格式 四阶龙格库塔 00.10.20.30.40.50.60.70.80.91
26、1.11.21.31.4 X 轴 图 5.2 步长为 0.1时各方法的预测值与精确值的比较 原 函 数 图 像 1.4 1.35 1.3 1.25 轴 1.2 1.15 1.1 1.05 1 00.10.20.30.40.50.60.70.80.91 X轴 图 5.3步长为 0.1 时原函数图像 在这里取步长 h 为 0.05，通过 MATLAB 程序的计算，相应的结果如下： 第页/ 共 221页0 表 5-2 h=0.05时三个方法与精确值的真值表 隆格库塔法求解常微分方程 步长Euler 法 相对误预估校正相对误经典四相对误 精确值 1.00000 0.051.00250 0.009111
27、.001250.01035 1.00123 0.01037 1.00484 0.101.00738 0.017111.004880.01954 1.00484 0.01958 1.01873 0.151.01451 0.024051.010760.02765 1.01071 0.02770 1.04082 0.201.02378 0.030011.018800.03473 1.01873 0.03479 1.07032 0.251.03509 0.035071.028890.04085 1.02880 0.04093 1.10653 0.301.04834 0.039301.040920.04
28、610 1.04082 0.04619 1.14881 0.351.06342 0.042771.054800.05053 1.05469 0.05063 1.19659 0.401.08025 0.045551.070440.05422 1.07032 0.05432 1.24933 0.451.09874 0.047721.087750.05724 1.08763 0.05734 1.30657 0.501.11880 0.049331.106660.05964 1.10653 0.05975 1.36788 0.551.14036 0.050451.127090.06150 1.1269
29、5 0.06161 1.00484 0.601.16334 0.051131.148950.06286 1.14881 0.06298 1.01873 0.651.18768 0.051431.172190.06379 1.17205 0.06391 1.04082 0.701.21329 0.051391.196740.06433 1.19659 0.06445 1.07032 0.751.24013 0.051061.222520.06453 1.22237 0.06465 1.10653 0.801.26812 0.050481.249490.06443 1.24933 0.06455
30、1.14881 0.851.29722 0.049691.277570.06408 1.27742 0.06419 1.19659 0.901.32735 0.048711.306730.06349 1.30657 0.06361 1.24933 0.951.35849 0.047591.336900.06272 1.33674 0.06283 1.30657 1.001.39056 0.046341.368040.06178 1.36788 0.06189 1.36788 差法差阶库差 0- 第页/ 共 221页1 隆格库塔法求解常微分方程 1.5 1.4 欧拉格式 1.3 改进欧拉格式 1
31、.2 四阶龙格库塔 1.1 精确值 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 00.10.20.30.40.50.60.70.80.91 轴 Y 00.10.20.30.40.50.60.70.80.91 X轴 图 5.5步长为 0.05 时原函数图像 实例 2 dy xe x y dx y(0)1 由题可知精确解为： 1 (x 2 2 2)e x 图 5.4 步长为 0.05 时各方法的预测值与精确值的比较 1.4 原 函 数 图 像 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1 第页/ 共 221页2 步长精确值 法差正法差阶库差
32、 0.10.90000.03180.90960.02140.90940.02160.9295 0.20.81920.06110.83590.04200.83560.04240.8726 0.30.75440.08760.77610.06140.77570.06190.8268 0.40.70270.11090.72770.07930.72720.07990.7903 0.50.66170.13100.68860.09560.68810.09630.7615 0.60.62940.14800.65720.11030.65670.11100.7387 0.70.60400.16210.63200.
33、12340.63150.12410.7209 0.80.58410.17370.61160.13490.61110.13550.7069 0.90.56860.18300.59510.14500.59460.14560.6959 1.00.55630.19050.58150.15380.58110.15440.6872 隆格库塔法求解常微分方程 当 x=0时，y(x)=0. 在这里取步长 h 为 0.1 ，通过 MATLAB 程序的计算，相应的结果 如下： 表 5-3 步长为 0.1 时各方法的预测值与精确值的比较(精确到 5 位小数 ) Euler相对误预估校相对误经典四相对误 第页/ 共
34、221页3 隆格库塔法求解常微分方程 原函数图像 1 0.95 0.9 0.85 0.8 轴 Y 0.75 0.7 0.65 0.6 0.55 00.10.20.30.40.50.60.70.80.91 X轴 图 5.6步长为 0.1 时原函数图像 0.95 0.9 0.85 各方法的预测值与精确值的比较 欧拉格式 改进的格式 四阶龙格库塔 精确值 0.8 轴 0.75 Y 0.7 0.65 0.6 0.55 0.10.20.30.40.50.60.70.80.91 X 轴 图 5.7步长为 0.1时各方法的预测值与精确值的比较 在这里取步长 h 为 0.05 ，通过 MATLAB 程序的计算
35、，相应的结果如下： 第页/ 共 221页4 隆格库塔法求解常微分方程 表 5-4 步长为 0.1 时各方法的预测值与精确值的比较(精确到 5 位小数 ) 步长Euler 法 相 对 误 差 预 估 校 正法 相 对 误 差 经 典 四 阶库 相 对 误 差 精确值 0.050.950000.01342 0.952450.01088 0.952420.01090 0.96292 0.100.904900.02650 0.909470.02158 0.909420.02163 0.92953 0.150.864280.03916 0.870670.03206 0.870610.03213 0.89
36、951 0.200.827740.05137 0.835670.04228 0.835590.04237 0.87256 0.250.794910.06307 0.804140.05219 0.804050.05230 0.84842 0.300.765450.07423 0.775760.06176 0.775660.06188 0.82682 0.350.739040.08482 0.750230.07096 0.750130.07110 0.80754 0.400.715410.09482 0.727300.07977 0.727190.07991 0.79035 0.450.69427
37、0.10422 0.706720.08817 0.706600.08832 0.77505 0.500.675390.11302 0.688250.09615 0.688130.09630 0.76146 0.550.658550.12123 0.671680.10370 0.671560.10386 0.74940 0.600.643520.12886 0.656830.11084 0.656710.11100 0.73871 0.650.630120.13593 0.643510.11756 0.643400.11773 0.72925 0.700.618190.14245 0.63157
38、0.12388 0.631450.12405 0.72087 0.750.607540.14847 0.620850.12981 0.620730.12998 0.71347 0.800.598050.15400 0.611210.13537 0.611100.13553 0.70691 0.850.589570.15908 0.602540.14058 0.602430.14073 0.70109 0.900.581980.16374 0.594700.14545 0.594600.14560 0.69593 0.950.575170.16802 0.587610.15001 0.58751
39、0.15016 0.69132 1.000.569040.17194 0.581170.15429 0.581070.15443 0.68719 第页/ 共 221页5 0.6 0.55 00.10.20.30.40.50.60.70.80.91 X轴 图 5.8步长为 0.05 时各方法的预测值与精确值的比较 第页/ 共 221页6 隆格库塔法求解常微分方程 六、模型的评价 本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库 塔法在求解初值问题的优越性.从上面的实例比较可知， 在计算精度上，四阶经典龙格 -库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所
40、以四阶经典龙 格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方 法也是最大， 改进欧拉方法其次， 欧拉方法计算量最小 .这样的结果， 说明了运用以上三 种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中， 可以根据实 际情况，选择这 3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙 格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一 般情况；而欧拉方法更主要的是适用于对y 的估计上，而精度则有所欠缺，以上各方法 的选择，都取决于具体的情况. 七、课程设计的总结与体会 本文着重采用隆格库塔法运用MATLAB
41、编程来求解微分方程，相比于欧拉法以及 预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量 不大、算法组织容易 .其次，每一次的课程设计总是让我学到了更多的知识，不论是 C+、 SPASS还是 MATLAB 软件，这些都让我学到了如何解决实际问题的好工具，通过这些 工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会： （1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从 一点一点做起 .以往训练得少只是把握的不牢靠.所以做起来感到有点吃力 .所以，无论什 么学科，一定要打好基础 . （2）程序设计要靠多练，多见识，那样形成一种编程
42、思维，我想对我是有很大好处的. 尤其像我这种平时学得不扎实的人. （3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没 有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力， 但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉 . 第页/ 共 221页7 隆格库塔法求解常微分方程 参 考 文 献 1李庆杨，王能超，易大义编数值分析（第四版）M ：华中科技大学出版社， 2006. 2姜启源，谢金星，叶俊编数学模型（第三版）M 北京：高等教育出版社，
43、 2005 3刘琼荪，数学实验 M ，高等教育出版社， 2004 4王建伟， MATLAB7.X程序设计 M ，中国水利水电出版社， 2007 5王高雄，周之铭等编 . 常微分方程（第三版）M: 高等教育出版社， 2006 6何坚勇编著 . 运筹学基础（第二版） M. 北京：清华大学出版社， 2008. 第页/ 共 221页8 隆格库塔法求解常微分方程 附录 附录一：显示欧拉法matlab 程序 % 欧拉法 clear all clc x=; y=; y1=; h=0.1; x=0:h:1; n=length(x); for i=1:n y(i)=f1(x(i); end figure(1)
44、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=fopen(data.txt,wt); fprintf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2) DT=abs(y-y1); 第页/ 共 221页9 隆格库塔法求解常微分方程 plot(x,DT) %1.建立导数函数文件 function z=f(x,y) z=y-2*x/y; %2.建立原函数文件 function z1=f1(x) z
45、1=(2*x+1)(1/2); 迭代 n 次的后退的欧拉格式matlab 程序 %1.建立导数函数文件 function z=f(x,y) z=y-2*x/y; %2.建立原函数文件 function z1=f1(x) z1=(2*x+1)(1/2); %迭代 n 次的后退的欧拉格式 clear all clc N=input(请输入迭代次数 :); x=; y=; y1=; h=0.1; x=0:0.1:1; n=length(x); for i=1:n y(i)=f1(x(i); end figure(1) plot(x,y,g-); hold on 第页/ 共 222页0 隆格库塔法求解
46、常微分方程 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=fopen(data.txt,wt); fprintf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2) DT=abs(y-y1); plot(x,DT) 附录二：预估校正法matlab 程序 %1.建立导数函数文件 function z=f(x,y) z=y-2*x/
47、y; %2.建立原函数文件 function z1=f1(x) z1=(2*x+1)(1/2); %预估校正法 % 欧拉法 clear all clc x=; y=; 第页/ 共 222页1 隆格库塔法求解常微分方程 y1=; y1(1)=1; y2=; y2(1)=1; h=0.1; x=0:0.1:1; n=length(x); for i=1:n y(i)=f1(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
48、),y2(j-1)+f(x(j),y1(j); end Y=x;y2; fid=fopen(data.txt,wt); fprintf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y2,r-); figure(2) DT=y-y2; plot(x,DT) 附录三：四阶龙格库塔法matlab 程序 % 四阶龙格库塔法 clear all clc h=0.2; 第页/ 共 222页2 隆格库塔法求解常微分方程 x=0:h:1; n=length(x); y=; y1=; y1(1)=1; for i=1:n y(i)=f1(x(i); end figure(
49、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=fopen(data1.txt,wt); fprintf(fid,%6.2f %12.4fn,Y); fclose(fid); plot(x,y1,r-); figure(2) DT=abs(y-y1); plot(x,DT)
