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

隆格库塔法求解常微分方程
摘要
科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.
本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.
总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.
关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；Matlab
Runge Kutta Method For Solving Ordinary Differential Equations
ABSTRACT
Problem solving 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 analytical 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 accuracy.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 methods, 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 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 convergence, the Runge Kutta method of order 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 预估较正法 (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 运用MATLAB软件对模型求解结果及析 (8)
6 模型的评价 (16)
7 课程设计的总结与体会 (16)
参考文献 (17)
附录 (18)
一、问题的提出
1.1 问题背景：
科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：
0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）
其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.
1.2 问题的具体内容
实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1
dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.
实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1
x dy xe y dx y -⎧=-⎪⎨⎪=⎩， 其
精确解为21(2)2x
x e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.
二、问题假设
2.1 假设数值方法本身的计算是准确的.
2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.
2.3 假设步长的增加不会导致舍入误差的严重积累.
三、符号系统
3.1 符号说明
符号含义
h选取的步长
*
K平均斜率
p精度的阶数
∆前后两次计算结果的偏差
y第n个节点的实验值
n
()
y x第n个节点的精确值
n
δ实验值与精确值的绝对误差
四、问题的分析
问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.
4.1欧拉格式
（1）显式欧拉格式
1(,)n n n n y y hf x y +=+ （2） 局部截断误差
22
211()()()()22
n n n h h y x y y y x o h ξ++''''-=≈= （3） （2）隐式欧拉格式
111(,)n n n n y y hf x y +++=+ （4）
局部截断误差
2
211()()()2
n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法
预估 1(,)n n n n y y hf x y +=+ (6)
校正 111[(,)(,)]2
n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2
n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2
p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)
4.3 四阶龙格库塔方法的格式（经典格式）
112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩
(10)
五、模型的建立与模型求解
5.1 隆格库塔法的基本原理
隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.
5.1.1 Taylor 级数
设初值问题 '00(,)
()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有
2
'
''1()()()()....2n n n n h y x y x hy x y x +=+++； （11） 其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)
j f x y 来描叙求导过程，即
(0)(1)
'(0)''
(1)
(1)(1)'''(2)
,f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）
(2)(2)
()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13） 根据上式，结果导出下面 Taylor 格式
2
'
''()1...2!!p
p n n n
n n h h y y hy y y p +=++++ （14）
其中一阶Taylor 格式为： '1n n n y y hy +=+
提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：
1
1(1)1(1)!
p p n n h y y y p ζ+++-+=+ （15） 因此它具有p 阶精度.
5.1.2 隆格库塔方法的基本思想
隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商
1()()n n y x y x h
+- 根据微分中值定理，这有
'1()()()n n n y x y x y x h h θ+-=+ （16） 利用所给方程 '(,)y f x y =
1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17） 设 平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.
再考察改进的Euler 格式，它可以改写成平均化的形式：
11212
11()2(,)
(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩
（18） 这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.
5.1.3 四阶的隆格库塔法
为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：
112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩
（19）
下面为其具体的算法流程图：
5.2 其他求解常微分边值问题算法的简介
5.2.1欧拉数值算法（显式）
微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：
(,)()dy
f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩
(20) 其中a ，b 为常数.
开始
输入x 0,y 0,h,N x 1=x 0+h
k 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)
k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)
y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6
n=1
输出x 1,y 1
n=N? 结束
n=n+1
x 0=x 1,y 0=y 1
否
是
图5.1 龙格-库塔法流程图
因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x
∂∂. 简单欧拉法是一种单步递推算法.简单欧拉法的公式如下所示：
1(,)n n n n y y hf x y +=+ (21)
简单欧拉法的算法过程介绍如下:
给出自变量x 的定义域[,]a b ，初始值0y 及步长h .
对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+
5.2.2欧拉数值算法（隐式）
隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：
111(,)n n n n y y hf x y +++=+ (22)
隐式欧拉法是一阶精度的方法，比它精度高的公式是：
111[(,)(,)]2n n n n n n h
y y f x y f x y +++=++ (23)
隐式欧拉的算法过程介绍如下.
给出自变量x 的定义域[,]a b ，初始值0y 及步长h .
对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+
得出1k y +.
5.2.3 欧拉预估-校正法
改进欧拉法是一种二阶显式求解法，其计算公式如下所示：
1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2
n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ (24)
四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：
1[,(,)]22
n n n n n n h h y y hf x y f x y +=+++ (25)
5.3 模型求解
5.3.1运用MATLAB 软件对模型求解结果及分析
用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：
实例一. 1
(0)1dy
y x dx y ⎧=-++⎪⎨⎪=⎩
由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1， 通过MATLAB 程序的计算，相应的结果如下：
表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)
初
值 Euler 法 相对误差 预估校正法 相对误差 经典四
阶库 相对误差
精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.0
1.38558 0.01294 1.36854 0.00048 1.36788 0.00000
1.36788
步长为0.1时的精确值与预测值的比较
精确值
欧拉法
改进欧拉格式
四阶龙格库塔
轴
Y
00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4
X轴
图5.2 步长为0.1时各方法的预测值与精确值的比较
原函数图像
轴
Y
X轴
图5.3 步长为0.1时原函数图像
在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：
表5-2 h=0.05时三个方法与精确值的真值表
步长Euler法相对误
差
预估校正
法
相对误
差
经典四
阶库
相对误
差
精确值
0 -- -- -- -- -- --
1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.24933
0.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.30657
1.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788
欧拉格式 改进欧拉格式四阶龙格库塔精确值
图5.4 步长为0.05时各方法的预测值与精确值的比较
1
1.051.11.151.21.25
1.31.351.4X 轴
Y 轴
原函数图像
图5.5 步长为0.05时原函数图像
实例2 (0)1x
dy xe y
dx y -⎧=-⎪⎨⎪=⎩
由题可知精确解为：
2
1(2)2x x e -+
当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：
步长Euler
法
相对误
差
预估校
正法
相对误
差
经典四
阶库
相对误
差
精确值
0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.7069
0.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.6959
1.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872
原函数图像
轴
Y
00.10.20.30.40.50.60.70.80.91
X轴
图5.6 步长为0.1时原函数图像
各方法的预测值与精确值的比较
欧拉格式
改进的格式
四阶龙格库塔
精确值
轴
Y
0.10.20.30.40.50.60.70.80.91
X轴
图5.7 步长为0.1时各方法的预测值与精确值的比较
在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：
表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)
步长Euler法相对误
差
预估校
正法
相对误
差
经典四
阶库
相对误
差
精确值
0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.69593
0.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.69132
1.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.68719
0.10.20.30.4
0.50.60.70.80.91
0.550.60.650.70.750.8
0.850.90.951X 轴
Y 轴
各方法下的预测值与精确值的比较
欧拉格式
改进欧拉格式四阶龙格库塔精确值
图5.8 步长为0.05时各方法的预测值与精确值的比较
六、模型的评价
本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.
七、课程设计的总结与体会
本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：
（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.
（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.
（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.
参考文献
[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.
[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005
[3] 刘琼荪，数学实验[M]，高等教育出版社，2004
[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007
[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006
[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.
附录附录一：显示欧拉法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)
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.4f\n',Y);
fclose(fid);
plot(x,y1,'r-');
figure(2)
DT=abs(y-y1);
plot(x,DT)
%1.建立导数函数文件
function z=f(x,y)
z=y-2*x/y;
%2.建立原函数文件
function z1=f1(x)
z1=(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-');
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.4f\n',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/y;
%2.建立原函数文件
function z1=f1(x)
z1=(2*x+1)^(1/2);
%预估校正法
%欧拉法
clear all
clc
x=[];
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),y2(j-1))+f(x(j),y1(j))]; end
Y=[x;y2];
fid=fopen('data.txt','wt');
fprintf(fid,'%6.2f %12.4f\n',Y);
fclose(fid);
plot(x,y2,'r-');
figure(2)
DT=y-y2;
plot(x,DT)
附录三：四阶龙格库塔法matlab程序
%四阶龙格库塔法
clear all
clc
n=length(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=fopen('data1.txt','wt');
fprintf(fid,'%6.2f %12.4f\n',Y);
fclose(fid);
plot(x,y1,'r-');
figure(2)
DT=abs(y-y1);
plot(x,DT)。