微分方程常用的两种数值解法：欧拉方法与龙格—库塔法
四川师范大学本科毕业论文
微分方程常用的两种数值解法：欧拉方法与龙
格—库塔法
学生姓名XXX
院系名称数学与软件科学学院
专业名称信息与计算科学
班级2006级 4 班
学号20060640XX
指导教师Xxx
四川师范大学教务处
二○一○年五月
微分方程常用的两种数值解法：欧拉方法与龙格—库塔法
学生姓名：xxx 指导教师：xx
【内容摘要】微分方程是最有生命力的数学分支，在自然科学的许多领域中，都
会遇到常微分方程的求解问题。当前计算机的发展为常微分方程的应用及理论研究提供了非常有力的工具，利用计算机解微分方程主要使用数值方法，欧拉方法和龙格——库塔方法是求解微分方程最典型常用的数值方法。本文详细研究了这两类数值计算方法的构造过程，分析了它们的优缺点，以及它们的收敛性，相容性，及稳定性。讨论了步长的变化对数值方法的影响和系数不同的同阶龙格—库塔方法的差别。通过编制C程序在计算机上实现这两类方法及对一些典型算例的结果分析比较，能更深切体会它们的功能，优缺点及适用场合，从而在实际应用中能对不同类型和不同要求的常微分方程会选取适当的求解方法。
关键词：显式单步法欧拉（Euler）方法龙格—库塔（Runge—Kutta）方法截断误差收敛性
Two commonly used numerical solution of differential
equations：Euler method and Runge - Kutta method Student Name: Xiong Shiying Tutor：Zhang Li
【Abstract】The differential equation is the most vitality branch in mathematics. In many domains of natural science, we can meet the ordinary differential equation solution question. Currently, the development of computer has provided the extremely powerful tool for the ordinary differential equation application and the fundamental research, the computer solving differential equation mainly uses value method. The Euler method and the Runge—Kutta method are the most typical commonly value method to solve the differential equation. This article dissects the structure process of these two kinds of values commonly value method to solve the analyses their good and bad points, to their astringency, the compatibility, and the stability has made the proof. At the same time, the article discuss the length of stride to the numerical method changing influence and the difference of the coefficient different same step Runge—kutta method. Through establishing C program on the computer can realize these two kind of methods, Anglicizing some models of calculate example result can sincerely realize their function, the advantage and disadvantage points and the suitable situation, thus the suitable solution method can be selected to solve the different type and the
different request ordinary differential equation in the practical application . Keywords:Explicit single-step process Euler method Runge—Kutta method truncation error convergence
目录
微分方程常用的两种数值解法：欧拉方法与龙格—库塔法
前言
常微分方程的形成与发展是和力学、天文学、物理学以及其他科学技术的发展密切相关的。数学其他分支的新发展，如复变函数、群、组合拓扑学等，都对常微分方程的发展产生了深刻的影响，当前计算机的发展更是为常微分方程的应用及理论研究提供了非常有力的工具。牛顿在研究天体力学和机械力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律。后来，法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王星的位置。这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的数学分支。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值，利用计算机解微分方程主要使用数值方法。 本文主要讨论一阶常微分方程初值问题
?????==0
0)()
,(y x y y x f dx
dy
（1.1） 在区间],[b a 上的数值解法，其中),(y x f 为关于x ，y 的已知函数，0y 为给定的初始值，将上述问题的精确解记为)(x y 。该问题常用的数值解法有：欧拉（Euler ）方法、龙格—库塔（Runge —Kutta ）方法及一些常用的线性多步法。本文重点介绍欧拉
（Euler ）方法和龙格—库塔（Runge —Kutta ）方法。并对这两种方法编制程序，体会它们的功能、优缺点及适用场合，对不同类型常微分方程会选取适当的求解方法。
１ 基本概念和准备知识
一阶常微分方程初值问题是：
(1.1.2)
)((1.1.1)
),(00?????==y x y y x f dx
dy
其中),(y x f 是平面上某个区域D 上的连续函数，式（1.1.1）的微分方程一般有无穷多个解，式（1.1.2）是确定解的初始条件，如果一元函数)(x y 对一切b x a ≤≤满
足
（1）D x y x ∈))(,(； （2）00)(y x y =；
（3）'y 存在而且))(,()('x y x f x y =； 则称)(x y 是初值问题（1.1）在区间],[b a 上的解。
误差：假定在计算1+n y 时，用到的前一步的值n y 是准确的(即))(n n x y y =，把用)(n x y 计算得到的近似值记为1+-
n y ，估计误差：1+n E = y (x n +1)1+-
-n y ，这种误差称为局部截断误差。如果不作这一假定，在每一步计算时除局部截断误差以外，还有由于前一步不准确而引起的误差，称为总体截断误差。
收敛性：对于解初值问题的数值方法，我们希望它产生的数值解收敛于初值问题的准确解，“收敛性”的一般定义为：
对于所有满足引理 1.1条件的初值问题（1.1），如果有一种显式单步法：
),,(1h y x h y y n n n n ?+=+
产生的近似解，对于任意固定的nh x x n +=0，均有0
lim →h )(n n x y y =，则称该显式单步
法是收敛的。
相容性：显式单步法（1.2.1）称为与原微分方程相容，如果
))(,()0),(,(x y x f x y x =? （1.2.3）
成立。并称式（1.2.3）为相容性条件。
稳定性：在实际计算中，一方面初始值0y 不一定精确，往往带有一定误差，同时由于计算机字长有限，在计算过程中有舍入误差，而且这种误差在式（1.2.1）的递推过程中会传递下去，对以后的结果产生影响。因此要考虑舍入误差的积累是否会得到控制，也即要考虑数值方法的稳定性。当∞→n 时，若舍入误差引起的后果是有限的，则可以认为该方法是数值稳定的。
2 欧拉方法
2.1欧拉方法简介
对常微分方程初值问题（1.1）用数值方法求解时，我们总是认为（1.1）的解存在且唯一。
欧拉方法是解初值问题的最简单、最原始的数值方法，它是显式单步法。下面介绍几种导出欧拉法的途径，每个途径皆可以推导出更为有效的数值方法。
（1）Taylor 展开
在n x 点将)()(1h x y x y n n +=+作Taylor 展开得：
)
(!
2)(')()
()(2
1n n n n n y h x hy x y h x y x y ε''++=+=+ ),(1+∈n n n x x ε 当h 充分小时略去误差项)("2
2
n y h ε，并注意到))(,()('n n n x y x f x y =，得))(,()()(1n n n n x y x hf x y x y +≈+，以n y 近似代替)(n x y ，以1+n y 近似代替)(1+n x y ，且用
“=”代替“≈”得差分方程初值问题：
???=+=+)(),(001x y y y x hf y y n n n n nh x x n +=0
，1,,2,1,0-=N n ，=N h
a
b - (2.1.1) 用式(2.1.1)求解初值问题（1.1）的方法称为欧拉方法。
（2）数值微分
由导数的定义知，对于充分小的h
))(,()(')
()(1n n n n n x y x f x y h
x y x y =≈-+
整理得))(,()()(1n n n n x y x hf x y x y +≈+，对此作相似的处理也可以得到欧拉方法(2.1.1)。
（3）数值积分
在区间],[1+n n x x 对))(,('x y x f y =积分得
dx x y x f x y x y n n
x x n n ?
++=+1
))(,()()(1 （2.1.2）
用数值积分的左矩形公式计算式（2.1.2）右端的积分，得
))(,()()(1n n n n x y x hf x y x y +≈+，于是同样可以得到欧拉方法(2.1.1)。
（4）多项式插值
利用解和其导数),(y x f 在n x 点的值n y ，),(n n y x f 作一次埃尔米特插值，得到关于h 的插值多项式：),()(n n n y x hf y h p +=，用)(h p 近似代替1+n y 就得到欧拉方法。
（5）待定系数法
在第n 步，已知n y 和),('n n n y x f y =，利用这两个值估计出下一步的1+n y ，将已知的值与估计值作线性组合：'1n n n y y y βα+=+，其中α，β 为待定系数。为确定这两个参数，要求这个估计值对c x y ≡)(和cx x y =)(（c 为常数）精确成立。如果
c x y ≡)(，则0)('=x y ，得到方程：0βα+=c c ，得1=α。如果cx x y =)(，则c x y =)('，
这样有：)(1ββα+=+=+n n n x c c cx cx ，说明h x x n n =-=+1β，这样估计值为：
),('1n n n n n n y x hf y hy y y +=+=+，即为欧拉方法。
欧拉方法的几何意义：
由点斜式得切线方程
()
),()()
(0000,0000y x f x x y dx
dy x x y y y x -+=-+=
等步长为h ，则h x x =-01，可由切线方程算出1y ：),(0001y x hf y y +=，逐步计算出)(x y y =在1+n x 点的值：),(1n n n n y x hf y y +=+， ,2,1,0=n ，用分段的折线逼近函数为“折线法”而非“切线法”，除第一个点是曲线上的切线，其它都不是。
2.2欧拉方法的截断误差，收敛性，相容性，稳定性
设)(n n x y y =，把)(1+n x y 在n x 处展开成泰勒级数，即
)(!
2)(')()(2
1n n n n y h x hy x y x y ε''++=+ ()1,+∈n n n x x ε
再由欧拉方法：
)(')())(,()(),(1n n n n n n n n n x hy x y x y x hf x y y x hf y y +=+=+=+
两式相减得欧拉方法的局部截断误差为：)("2
12
1n n y h E ε=
+，若)(x y 在],[b a 上充分
光滑，且令=M []
b a x ,max ∈)("x y ，则1+n E )(2
122
h M h ο=≤
，故欧拉方法是一阶方法或具有一阶精度。
欧拉方法的增量函数),,(h y x ?就是f ，由引理 1.3、引理 1.4知当f 满足Lipschitz 条件时欧拉方法是收敛的而且是相容的。
用欧拉方法求解典型方程（1.2.4）的计算公式为：n n n n y h hy y y )1(1λλ+=+=+，有 n n h δλδ)1(1+=+。要让n n δδ<+1，必须有11<+h λ，因此欧拉方法的绝对稳定域为 ：11<+h λ，当λ为实数时，绝对稳定区间为)0,2(-。在复平面上，11<+h λ是以1为半径、以1-为圆心的内部。
3 龙格—库塔法
3.1 龙格—库塔法的基本思想
为了导出龙格—库塔法的一般公式，我们取如下的线性组合形式：
i s
i i n n k b h y y ∑=++=11
（3.2.3）
其中 ),(1
1
j i j ij n i n i k a h y h c x f k ∑-=++= （3.2.4）
即 ???????+++=++==
),(),()
,(232131331
21221k ha k ha y h c x f k k ha y h c x f k y x f k n n n n n n
1b ,s b b ,,2 ；s c c c c ,,,,321 ；a 21,a 31,, 1-ss a 除c 1=0外均为待定系数。显然用公式（3.2.3）每计算一个新值1+n y 要计算函数f 的值s 次，又因每个j k 都能以一种明显的方式由1k ,2k ,, 1-j k 计算出来，故将公式（3.2.3）称为s 级显式龙格—库塔法。s 级显式龙格—库塔法又可以写成下面既简洁又直观的阵列形式：
2c 21a 3c 31a 32a
s c 1s a 2s a 1-ss a
1b 2b 1-s b s b
3.2 二、三、四级龙格—库塔法
常见的二级二阶龙格—库塔法有：
中点方法（取2
1
,0,122112====c a b b ）
???
?
???++==+=+)2,2()
,(12121k h y h x f k y x f k hk y y n
n n n n n 二阶Heun 方法（取3
2
,41,4322112====c a b b ）
???
?
???++==++=+)
32,32(),()3(4
121411
hk y h x f k y x f k k k h y y n n n n n n 取s = 3，完全仿照上述方法推导出三阶龙格—库塔公式。这时参数满足下列
条件
???
??
??
??
??==+=+=++613
121
1322323322
23322321a c b c b c b b c b c b b b 这个方程组解也不唯一，从而可以得到不同的三级三阶龙格—库塔法。两个较为常
见的三级三阶龙格—库塔法是：
Kutta 方法（取2c =21
，3c = 1，a 21 = 21，a 31 = 1-，a 32 = 2，1b = 6
1，2b =
64，3b = 6
1） 将它代入（3.2.3）得
??????
??
?
+-+=++==+++=+)
2,()
2
,2()
,()4(6
213
1213211
hk hk y h x f k k h y h x f k y x f k k k k h y y n n n
n n n n n (3.3.2) 三阶Heun 方法（取2c =31，3c = 32，a 21 = 21
，a 31 = 0，a 32 = 32，1b = 4
1，2b
= 0，3b =4
3
）将它代入（3.2.3）得
???
??
?
??
???
++=++==++=+)32,32()3,3(),()4
341(23121311
hk y h x f k k h y h
x f k y x f k k k h y y n n n n n n n n 通常人们所说的龙格—库塔法是指四阶而言的。取s=4，我们同样可以仿照二
阶的情形推导出此公式。这样就得到如下含12个未知数但仅有10个方程的方程组
???
?????
???
?????
?????=
++=
++=++=++=
++=++=+++=++=+81)(121)(6
1
)(413121
1
43342244323234323422243222343342243223344333322244233222442322432144342413
3231a c a c c b a c c b a c a c b a c b a c a c b a c b c b c b c b c b c b c b c b c b c b b b b b c
a a a c a a
此方程组的解同样不是唯一的，从而有许多不同的四级四阶龙格—库塔法，最常用的两个四阶公式是：
经典龙格—库塔法
????
?
?????
???
++=++=++==++++=+)
,()2
,21()2,21(),()22(6
3423
12143211
hk y h x f k k h y h x f k k h y h x f k y x f k k k k k h y y n n n n n
n n n n n
（3.3.3）
RKG(Runge-Kutta-Gill)方法
?????
?
?
??
????
++-+=-+-++=++==+++-++=+)
22222,()
222212,2()2,2()
,(])22()22([6
3242
1312
143211
hk hk y h x f k hk hk y h x f k k h y h x f k y x f k k k k k h y y n n n n n n n n n n
(3.3.4)
经典龙格—库塔方法的几何解释：
将（3.3.3）式中的j k 改记为j f ，从而得到经典龙格—库塔方法的计算公式为
)(6
43211f f f f h
y y n n ++++
=+ 只考虑区间],[10x x 上的解曲线)(x y y =。由(3.3.3)式可知，曲线)(x y y =在0x 处的斜率))(,()('000x y x f x y =1f =，在h x x +=01处的斜率4111))(,()('f x y x f x y ≈=，而
2f 和3f 则是在中点2
021h
x x +
=处的斜率))(,()('1211x y x f x y =的两个近似值，如图2：
另一方面有：?=-1
0))(,()()(01x x dx x y x f x y x y ，用Simpson 积分公式逼近此式右
端的积分，得
))](,())(,(4))(,([6
)()(112110001x y x f x y x f x y x f h
x y x y ++≈-
上式右端要计算),(y x f 的三个值。由前面的结果，可取100))(,(f x y x f =，
411))(,(f x y x f ≈，)(2
1
))(,(312121f f x y x f +≈
，见图3：
代入后得到
]2
)(4[6432101f f f f h
y y ++++=)22(643210f f f f h y ++++= （3.3.5）
另外，从(3.3.3)式容易看出，当),(y x f 与y 无关时，公式(3.3.5)实际上就是用标准的Simpson 积分公式计算积分dx y x f x x ?1
0),(得到的，因此，可以将经典龙格——库
塔方法视为是用Simpson 积分公式的推广形式计算积分?+1
),(n n
x x dx y x f 而得到的数值
方法。
经典龙格—库塔法的稳定性：
增量函数：)22(6
1
),,(4321k k k k h y x +++=?，其中
y
h hy y k hy y k y
k 23232214
1
2121
λλλλλλ++=+==
y h y h hy y k 34232424
1
21λλλλ+++=
代入后得
=),,(h y x ?y h y h y h hy y 3
42323224
16121λλλλλ+
+++ 于是有
y
???=3423232241
6121h h h h λλλλλ++++
从而得出龙格—库塔方法的绝对稳定区间为 3
423232241
6121h h h h λλλλλ++++≤1
0785.2≤≤-h λ
经典龙格—库塔方法的算法框图：
4 应用举例建模及其分析
4.1 欧拉方法解题及其数学模型
4.1.1 问题提出
例4.1.1 假定某公司的净资产因资产本身产生了利息而以4%的年利率增长，同时，该公司以每年100万的数额支付职工工资。净资产的微分方程为
10004.0-=w dt
dw
（t 以年为单位） 分别以初始值,3500)0(,2500)0(,1500)0(===z y x
问题：用Euler 公式预测公司24后的净资产趋势。
4.1.2 模型建立
分析：这是求微分方程的数值积分，为的是预测公司24年后的净资产趋势。 确定变量：设净资产是时间的微分函数，不妨设变量t 为时间（以年为单位）。 设n w 为第n 年的净资产，1+n w 为第n+1年的净资产，利息以每年4%的速度增长，且公司每年支付职工工资为100万，则第n+1年的净资产增长数额为)10004.0(-n w ，由于已得第n 年的净资产为n w ，则第n 年的净资产加上第n+1年的净资产增长数额就得到第n+1年的净资产1+n w 。
归纳公式：.1,10004.1)10004.0(1=-=-+=+h w w h w w n n n n （4-1） 确定其微分方程的标准形式，这就是例4.1.1的微分方程模型。
4.1.3 解决问题
0w 分别以3500,2500,1500
000===z y x 代入，x 表示利息赢利低于工资支出的数额，y 表示利息赢利与工资支出平衡的数额，z 表示利息赢利高于工资支出的数额，
计算结果见表1.
资产为负值；当利息赢利与工资的支出平衡时，公司的净资产每年保持不变；当利息赢利超过工资的支出，公司的净资产稳步增长。
再用欧拉方法求解，每计算一步
y，依次需要计算1次、2次和4次函数f的
i
值，为了比较在计算量相同的条件下近似解的精度，步长分别取0.025、0.05、0.1。我对欧拉的C程序做了一些优化，加入了计算误差，计算结果如下（][i
y表示近似解、][1i
y表示准确解、][i e表示误差）：
This is Euler method:
please input a b and y0：0 1 1 please input N：40
x[0]=0.000000 y[0]=1.000000 y1[0]=1.000000 e[0]=0.000000
x[1]=0.025000 y[1]=1.025000 y1[1]=1.024695 e[1]=-0.000305
x[2]=0.050000 y[2]=1.049405 y1[2]=1.048809 e[2]=-0.000597
x[3]=0.075000 y[3]=1.073258 y1[3]=1.072381 e[3]=-0.000878
x[4]=0.100000 y[4]=1.096596 y1[4]=1.095445 e[4]=-0.001151
x[5]=0.125000 y[5]=1.119451 y1[5]=1.118034 e[5]=-0.001417
x[6]=0.150000 y[6]=1.141854 y1[6]=1.140175 e[6]=-0.001679
x[7]=0.175000 y[7]=1.163832 y1[7]=1.161895 e[7]=-0.001937
x[8]=0.200000 y[8]=1.185410 y1[8]=1.183216 e[8]=-0.002194
x[9]=0.225000 y[9]=1.206609 y1[9]=1.204159 e[9]=-0.002450
x[10]=0.250000 y[10]=1.227451 y1[10]=1.224745 e[10]=-0.002706
x[11]=0.275000 y[11]=1.247954 y1[11]=1.244990 e[11]=-0.002964
x[12]=0.300000 y[12]=1.268134 y1[12]=1.264911 e[12]=-0.003223
x[13]=0.325000 y[13]=1.288009 y1[13]=1.284523 e[13]=-0.003486
x[14]=0.350000 y[14]=1.307593 y1[14]=1.303841 e[14]=-0.003753
x[15]=0.375000 y[15]=1.326900 y1[15]=1.322876 e[15]=-0.004024
x[16]=0.400000 y[16]=1.345941 y1[16]=1.341641 e[16]=-0.004300
x[17]=0.425000 y[17]=1.364730 y1[17]=1.360147 e[17]=-0.004583
x[18]=0.450000 y[18]=1.383278 y1[18]=1.378405 e[18]=-0.004873
x[19]=0.475000 y[19]=1.401594 y1[19]=1.396424 e[19]=-0.005170
x[20]=0.500000 y[20]=1.419689 y1[20]=1.414214 e[20]=-0.005475
x[21]=0.525000 y[21]=1.437572 y1[21]=1.431782 e[21]=-0.005790
x[22]=0.550000 y[22]=1.455251 y1[22]=1.449138 e[22]=-0.006113
x[23]=0.575000 y[23]=1.472735 y1[23]=1.466288 e[23]=-0.006447
x[24]=0.600000 y[24]=1.490032 y1[24]=1.483240 e[24]=-0.006792
x[25]=0.625000 y[25]=1.507149 y1[25]=1.500000 e[25]=-0.007149
x[26]=0.650000 y[26]=1.524093 y1[26]=1.516575 e[26]=-0.007518
x[27]=0.675000 y[27]=1.540872 y1[27]=1.532971 e[27]=-0.007900
x[28]=0.700000 y[28]=1.557490 y1[28]=1.549193 e[28]=-0.008297
x[29]=0.725000 y[29]=1.573955 y1[29]=1.565248 e[29]=-0.008708
x[30]=0.750000 y[30]=1.590273 y1[30]=1.581139 e[30]=-0.009134
x[31]=0.775000 y[31]=1.606449 y1[31]=1.596872 e[31]=-0.009577
x[32]=0.800000 y[32]=1.622489 y1[32]=1.612452 e[32]=-0.010037
x[33]=0.825000 y[33]=1.638397 y1[33]=1.627882 e[33]=-0.010515
x[34]=0.850000 y[34]=1.654180 y1[34]=1.643168 e[34]=-0.011013
x[35]=0.875000 y[35]=1.669842 y1[35]=1.658312 e[35]=-0.011530
x[36]=0.900000 y[36]=1.685388 y1[36]=1.673320 e[36]=-0.012068
x[37]=0.925000 y[37]=1.700823 y1[37]=1.688194 e[37]=-0.012629
x[38]=0.950000 y[38]=1.716151 y1[38]=1.702939 e[38]=-0.013212
x[39]=0.975000 y[39]=1.731376 y1[39]=1.717556 e[39]=-0.013820
x[40]=1.000000 y[40]=1.746504 y1[40]=1.732051 e[40]=-0.014453
this is improve Euler method
please input a b and y0:：0 1 1 please input N：20
x[0]=0.000000 y[0]=1.000000 y1[0]=1.000000 e[0]=0.000000
x[1]=0.050000 y[1]=1.048869 y1[1]=1.048809 e[1]=-0.000060
x[2]=0.100000 y[2]=1.095561 y1[2]=1.095445 e[2]=-0.000116
x[3]=0.150000 y[3]=1.140345 y1[3]=1.140175 e[3]=-0.000169
x[4]=0.200000 y[4]=1.183437 y1[4]=1.183216 e[4]=-0.000221
x[5]=0.250000 y[5]=1.225017 y1[5]=1.224745 e[5]=-0.000272
x[6]=0.300000 y[6]=1.265236 y1[6]=1.264911 e[6]=-0.000325
x[7]=0.350000 y[7]=1.304219 y1[7]=1.303841 e[7]=-0.000378
x[8]=0.400000 y[8]=1.342074 y1[8]=1.341641 e[8]=-0.000433
x[9]=0.450000 y[9]=1.378896 y1[9]=1.378405 e[9]=-0.000492
x[10]=0.500000 y[10]=1.414766 y1[10]=1.414214 e[10]=-0.000553
x[11]=0.550000 y[11]=1.449756 y1[11]=1.449138 e[11]=-0.000618
x[12]=0.600000 y[12]=1.483927 y1[12]=1.483240 e[12]=-0.000687
x[13]=0.650000 y[13]=1.517337 y1[13]=1.516575 e[13]=-0.000762
x[14]=0.700000 y[14]=1.550035 y1[14]=1.549193 e[14]=-0.000841
x[15]=0.750000 y[15]=1.582066 y1[15]=1.581139 e[15]=-0.000927
x[16]=0.800000 y[16]=1.613472 y1[16]=1.612452 e[16]=-0.001020
x[17]=0.850000 y[17]=1.644289 y1[17]=1.643168 e[17]=-0.001121
x[18]=0.900000 y[18]=1.674551 y1[18]=1.673320 e[18]=-0.001230
x[19]=0.950000 y[19]=1.704288 y1[19]=1.702939 e[19]=-0.001349
x[20]=1.000000 y[20]=1.733529 y1[20]=1.732051 e[20]=-0.001478
4.1.4 结论
这道题用普通的微分方程也能列式求解，关键在于如何预测若干年后的净资产趋势，用普通的微分方程就无法进行预测，且人工计算量相当大，这里我使用欧拉方法可以计算出精度以及误差，通过电脑运行程序就可以预测出若干年后的净资产情况，欧拉方法的使用使得解题更加方便且精确。
4.2 龙格—库塔解题及其数学模型
4.2.1 问题提出
例4.2.1 两种果树寄生虫，其数量分别是),
u
=其中一种寄生虫以吃
t
u=
v
(
v
),
(t
另一种寄生虫为生，两种寄生虫的增长函数如下列常微分方程组所示：
???
?
??
???==--=--=2.1)0(6.1)0(001
.0)151(06.045.0)20
1(09.0v u uv v v d d uv u u d d t v
t
u
问题：预测3年后这一对寄生虫的数量。
4.2.2 模型建立
分析：这是一个典型的常微分方程例题，要求预测出3年后这对寄生虫的数量。 确定变量：假定时间是寄生虫数量的积分函数，不妨设变量t 为时间（以年为单位）。
由题知有两种寄生虫u 和v ，u 寄生虫年增长函数为uv
u u d d t u 45.0)201(09.0--=，v 寄生虫年增长函数为uv
v v d d t v 001.0)151(06.0--=，初始值：u 寄生虫为1.6，v 寄生虫
为1.2，由于其中一种寄生虫以吃另一种寄生虫为生，我们可建立u 和v 的关联函
数f(u,v),g(u,v)。
归纳后得到公式：
???
????
--=--=uv
v v v u g uv u u v u f 001.0)151(06.0),(45.0)20
1(09.0),( （4-2）
（4-2）即为例4.2.1所述问题的微分方程模型。
4.2.3 解决问题
在本题中),,(),,(),,(),,(v u g v u t g v u f v u t f ==用Euler 预估—校正公式
????
?
?+???? ??=???? ??++),(),(11n n n n n n n n v u g v u f h v u v u ???????????? ??+???? ??+???? ??=???? ??++++++)11()11(),(),(11,,2n n n n n n n n n n n n v u g v u f v u f v u f h v u v u 取h=1，计算结果如表2.
please input a b and y0:：0 1 1 please input N ： 10
x[0]=0.000000 y[0]=1.000000 y1[0]=1.000000 e[0]=0.000000 x[1]=0.100000 y[1]=1.095446 y1[1]=1.095445 e[1]=-0.000000 x[2]=0.200000 y[2]=1.183217 y1[2]=1.183216 e[2]=-0.000001 x[3]=0.300000 y[3]=1.264912 y1[3]=1.264911 e[3]=-0.000001 x[4]=0.400000 y[4]=1.341642 y1[4]=1.341641 e[4]=-0.000001 x[5]=0.500000 y[5]=1.414215 y1[5]=1.414214 e[5]=-0.000002 x[6]=0.600000 y[6]=1.483242 y1[6]=1.483240 e[6]=-0.000002 x[7]=0.700000 y[7]=1.549196 y1[7]=1.549193 e[7]=-0.000003 x[8]=0.800000 y[8]=1.612455 y1[8]=1.612452 e[8]=-0.000004 x[9]=0.900000 y[9]=1.673324 y1[9]=1.673320 e[9]=-0.000004 x[10]=1.000000 y[10]=1.732056 y1[10]=1.732051 e[10]=-0.000005
4.2.3 结论
从上面的结果可以分析，用每一种方法计算节点x =0.1、0.2 0.9、1.0，上的值y 都需要计算4次),(y x f 的值，即它们的计算量基本相同，其结果是经典的龙格—库塔方法的精度最好，龙格—库塔方法的推导是基于Taylor 级数方法的，因而在使用高阶龙格—库塔方法计算时可以很精确的推算出寄生虫每一年的数量。
4.3 分别使用二阶、三阶龙格—库塔法解初值问题
对一些特殊的微分方程，使用欧拉方法和低阶的龙格—库塔方法也能达到很高的精度，例如：微分方程的解析解是一次函数，则用欧拉方法求得的数值解与准确解相符，微分方程的解析解是二次函数，则用二阶龙格—库塔方法求得的数值解与准确解相符。微分方程的解析解是三次多项式，用三阶龙格—库塔方法求得的数值解与准确解相符。
4.3.1 建立模型
建立初值问题模型：
?????=+=1
)0(1
2y x dx
dy
)10(≤≤x (4.1)
4.3.2 用不同的龙格—库塔法解（4.1）初值问题模型]1[]2[ （1）用二阶龙格库塔方法求解初值问题（4.1） 结果如下：
this is the second-order Runge-Kutta(Heun) method please input a b and y0： 0 1 1 please input N ：10
龙格库塔解微分方程
《数值分析》课程实验报告 一、实验目的 1．掌握用MA TLAB求微分方程初值问题数值解的方法； 2．通过实例学习微分方程模型解决简化的实际问题； 3．了解龙格库塔方法的基本思想。 二、实验内容 用龙格库塔方法求下列微分方程初值问题的数值解 y’=x+y y(0)=1 0 end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 四、实验小结 通过实验结果分析可知，龙格库塔方法求解常微分方程能获得比较好的数值解，龙格库塔方法的数值解的精度较高，方法比较简便易懂。 姓名：樊元君学号：02 日期： 一、实验目的 掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。掌握使用MATLAB程序求解常微分方程问题的方法。 ： 二、实验内容 1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。 实验中以下列数据验证程序的正确性。 求，步长h=。 * 2、实验注意事项 的精确解为，通过调整步长，观察结果的精度的变化 ^ ) 三、程序流程图： ●改进欧拉格式流程图： ~ | ●四阶龙格库塔流程图： ] 四、源程序： ●改进后欧拉格式程序源代码： function [] = GJOL(h,x0,y0,X,Y) format long h=input('h='); … x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); X=input('X=');Y=input('Y='); n=round((Y-X)/h); \ i=1;x1=0;yp=0;yc=0; for i=1:1:n x1=x0+h; yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);% · yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);% y1=(yp+yc)/2; x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y); : end end ●四阶龙格库塔程序源代码： function [] = LGKT(h,x0,y0,X,Y) 。 format long h=input('h='); x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); " X=input('X=');Y=input('Y='); n=round((Y-X)/h); i=1;x1=0;k1=0;k2=0;k3=0;k4=0; for i=1:1:n ~ x1=x0+h; k1=-x0*y0^2;%k1=y0-2*x0/y0;% k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);% … y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);% x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y); end · end matlab编的4阶龙格库塔法解微分方程的程序 2010-03-10 20:16 function varargout=saxplaxliu(varargin) clc,clear x0=0;xn=1.2;y0=1;h=0.1; [y,x]=lgkt4j(x0,xn,y0,h); n=length(x); fprintf(' i x(i) y(i)\n'); for i=1:n fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end function z=f(x,y) z=-2*x*y^2; function [y,x]=lgkt4j(x0,xn,y0,h) x=x0:h:xn; n=length(x); y1=x; y1(1)=y0; for i=1:n-1 K1=f(x(i),y1(i)); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4); end y=y1; 结果： i x(i) y(i) 1 0.0000 1.0000 2 0.1000 0.9901 3 0.2000 0.9615 4 0.3000 0.9174 5 0.4000 0.8621 6 0.5000 0.8000 7 0.6000 0.7353 8 0.7000 0.6711 9 0.8000 0.6098 10 0.9000 0.5525 11 1.0000 0.5000 12 1.1000 0.4525 13 1.2000 0.4098 1.经典四阶龙格库塔法解一阶微分方程组 1.1运用四阶龙格库塔法解一阶微分方程组算法分析 ),,(1k k k y x t f f =， )2,2,2(112g h y f h x h t f f k k k +++= )2 ,2,2(223g h y f h x h t f f k k k +++= ),,(334hg y hf x h t f f k k k +++= ),,(1k k k y x t g g = )2,2,2(112g h y f h x h t g g k k k +++= )2,2,2(223g h y f h x h t g g k k k +++= ),,(334hg y hf x h t g g k k k +++= ) 22(6 )22(6 43211 43211g g g g h y y f f f f h x x k k k k ++++=++++=++ 1k k t t h +=+ 经过循环计算由 推得 …… 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为() N O h ,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精 准，稳定，且易于编程。 000,,t x y ()()111222,,,,t x y t x y （1-1） （1-2） （1-3） （1-4） （1-5） （1-6） （1-7） （1-8） （1-9） （1-10） 1.2经典四阶龙格库塔法解一阶微分方程流程图 图1-1 经典四阶龙格库塔法解一阶微分方程流程图 1.3经典四阶龙格库塔法解一阶微分方程程序代码： #include 龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下 ()()()11234121324 3226,,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 +?=++++???=????=++? ???????=+? ?????=++? 1龙格-库塔法解一阶ODE 对于形如()()0, dy f x y a x b dx y a y ?=<≤???=? 的一阶ODE 初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例 ()2 0101dy x y x dx y y ?=-<≤???=? 取步长h=0.1 ，此处由数学知识可得该方程的精确解为y =。在这里利用MATLAB 编程，计算数值解并与精确解相比，代码如下： (1)写出微分方程，便于调用和修改 function val = odefun( x,y ) val = y-2*x/y; end (2)编写runge-kutta 方法的函数代码 function y = runge_kutta( h,x0,y0 ) k1 = odefun(x0,y0); k2 = odefun(x0+h/2,y0+h/2*k1); k3 = odefun(x0+h/2,y0+h/2*k2); k4 = odefun(x0+h,y0+h*k3); y = y0+h*(k1+2*k2+2*k3+k4)/6; end (3)编写主函数解微分方程,并观察数值解与精确解的差异clear all h = 0.1; x0 = 0; y0 = 1; x = 0.1:h:1; y(1) = runge_kutta(h,x0,y0); for k=1:length(x) x(k) = x0+k*h; y(k+1) = runge_kutta(h,x(k),y(k)); end z = sqrt(1+2*x); plot(x,y,’*’); This program is to solve the initial value problem of following system of differential equations: dx/dt=x+2*y,x(0)=0, dy/dt=2*x+y,y(0)=2, x and y are to be calculated ****************************************************************************/ #include<> #include<> #define steplength //步长?è可¨|根¨′据Y需¨¨要°a调ì??整; #define FuncNumber 2 //FuncNumber为a微?é分¤方¤程¨?的ì数oy目; void main() { float x[200],Yn[20][200],reachpoint;int i; x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初值|ì条?件t; reachpoint=; //所¨′求¨?点ì可¨|根¨′据Y需¨¨要°a调ì??整; void rightfunctions(float x ,float *Auxiliary,float *Func); void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]); Runge_Kutta(x ,reachpoint, Yn); printf("x "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",x[i]); printf("\nY1 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[0][i]); printf("\nY2 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[1][i]); getchar(); } void rightfunctions(float x ,float *Auxiliary,float *Func)//当ì?à右?¨°方¤程¨?改变à时o?à，ê需¨¨要°a改变à; { Func[0]=Auxiliary[0]+2*Auxiliary[1]; Func[1]=2*Auxiliary[0]+Auxiliary[1]; } void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]) { int i,j; float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber]; for(i=0;i<=(reachpoint-x[0])/steplength;i++) { for(j=0;j /*************************************************************************** This program is to solve the initial value problem of following system of differential equati ons: dx/dt=x+2*y,x(0)=0, dy/dt=2*x+y,y(0)=2, x(0.2) and y(0.2) are to be calculated ****************************************************************************/ #include 四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程：cos ，初始值y(0)=0，求 y t 解区间[0 10]。 MATLAB程序： %%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost %%%%%%%%%%% y(0)=0, 0≤t≤10，h=0.01 %%%%%%%%%%% y=sint h=0.01; hf=10; t=0:h:hf; y=zeros(1,length(t)); y(1)=0; F=@(t,y)(cos(t)); for i=1:(length(t)-1) k1=F(t(i),y(i)); k2=F(t(i)+h/2,y(i)+k1*h/2); k3=F(t(i)+h/2,y(i)+k2*h/2); k4=F(t(i)+h,y(i)+k3*h); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h; end subplot(211) plot(t,y,'-') xlabel('t'); ylabel('y'); title('Approximation'); span=[0,10]; p=y(1); [t1,y1]=ode45(F,span,p); subplot(212) plot(t1,y1) xlabel('t'); ylabel('y'); title('Exact'); 图1 ②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。 MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程 %%%%%%%%%%% x'(t)=(t*x-x^2)/t^2 %%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128 %%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt) h=1/128; %%%%% 步长 tf=3; t=1:h:tf; x=zeros(1,length(t)); x(1)=2; %%%%% 初始值 《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一 例应用之。 【实例】采用龙格-库塔法求微分方程： ?? ?==+-=0 , 0)(1 '00 x x y y y 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为： )22(6 3211k k k h y y n n +++=+ 其中： ?????????++=++=++ ==) ,() 21 ,21()21 ,21() ,(34 23121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图： 2.1、四阶龙格－库塔（R-K ）方法流程图： 2.2、实例求解流程图： 3、源程序代码 3.1、四阶龙格－库塔（R-K）方法源程序： function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin) %Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x)) %此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式 %函数 f(x,y): fun %自变量的初值和终值：x0, xt %y0表示函数在x0处的值，输入初值为列向量形式 %自变量在[x0,xt]上取的点数：PointNum %varargin为可输入项，可传适当参数给函数f(x,y) %x：所取的点的x值 %y：对应点上的函数值 if nargin<4 | PointNum<=0 PointNum=100; end if nargin<3 y0=0; end y(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长 x=x0+[0:(PointNum-1)]'*h; %得x向量值 for k=1:(PointNum)%迭代计算 f1=h*feval(fun,x(k),y(k,:),varargin{:}); f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:}); f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:}); f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:}); f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end 3.2、实例求解源程序： %运行四阶R-K法 蒙特卡洛法求定积分： 整体思路，蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对函数进行求值，进而除以所取次数并乘以其区间，即为所积值。 Step 1： 在执行程序前，打开matlab，执行以下操作打开File—>New—>Function,并点击Function，弹出Editor窗口，将其中内容修改为 function [ y ] = f( x ) y=cos(x)+2.0; end （如图所示）。 Step 2： 在Editor窗口中点击File—>Save As，保存至所建立的文件夹内，保存名称必须为f.m。 Step 3： 回到Matlab程序中，将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。Step 4： 在Command Window里输入以下源程序。 源程序： >> n=10^6; %用来表示精度，表示需要执行次数。 >> y=0; %初始化y=0，因为f(x)=cos(x)+2.0，下面出现y=y+f(x),故尔y用来统计计算出的所有f(x)值的和。 >> count=1; %从1开始计数，显然要执行n次。 >> while count<=n %当执行次数少于n次时，继续执行 x=unifrnd(0,4); %在matlab中，unifrnd是在某个区间内均匀选取实数（可为小数或整 数），表示x取0到4的任意数。 y=y+f(x); %求出到目前为止执行出的所有f(x)值的和，并用y表示 count=count+1; %count+1表示已执行次数再加一，将来通过与n的比对，判断是否执行下一次 end %如果执行次数已达n次，那么结束循环 z=y/n*4 %用y除以执行次数n,那么取得平均值,并用它乘以区间长度4,即可得到z 的近似值 z=7.2430 由于matlab中不能使用θ字符，故用z代替，即可得到θ=7.2395。 在执行过程中，发现每一次执行结果都会不一样，这是因为是随机选取，所以不同次数的计算结果会有偏差，但由于执行次数很多，从概率的角度来讲都是与真实值相近似的值。 《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一 例应用之。 【实例】采用龙格-库塔法求微分方程： ? ? ?==+-=0 , 0)(1 '00x x y y y 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分 方程的初值问题。其算法公式为： )22(6 3211k k k h y y n n +++ =+ 其中： ?????????++=+ + =++==) ,() 21 ,21()21 ,21(),(34 23121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图： 2.1、四阶龙格－库塔（R-K ）方法流程图： 2.2、实例求解流程图： 3、源程序代码 3.1、四阶龙格－库塔（R-K）方法源程序： function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin) %Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x)) %此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式 %函数 f(x,y): fun %自变量的初值和终值：x0, xt %y0表示函数在x0处的值，输入初值为列向量形式 %自变量在[x0,xt]上取的点数：PointNum %varargin为可输入项，可传适当参数给函数f(x,y) %x：所取的点的x值 %y：对应点上的函数值 if nargin<4 | PointNum<=0 PointNum=100; end if nargin<3 y0=0; end y(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长 x=x0+[0:(PointNum-1)]'*h; %得x向量值 for k=1:(PointNum) %迭代计算 f1=h*feval(fun,x(k),y(k,:),varargin{:}); f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:}); f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:}); f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:}); f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end 3.2、实例求解源程序： %运行四阶R-K法 一.作业： 用四阶龙格—库塔法求下列方程： ()()()1010100y x y x y '=-≤≤???=??步长0.1h =。并与解析解2 51x y e -=-比较。 二.程序 #include 左侧为数值解,右侧为解析解. 求解常微分方程组初值问题的 龙格库塔法分析及其C 代码 1、概 述 由高等数学的知识可知，一些特殊类型的常微分方程（组）能够求出给定初始值的解析解，而在科学与工程问题中遇到的常微分方程（组）往往是极其复杂的，要想求得其给定初始值的解析解就变得极其困难，甚至是得不到解析解。尽管如此，在研究实际问题时，往往只需要获得若干点上的近似值就行了，这就说明研究常微分方程（组）的数值解法是很有必要的。求解常微分方程（组）的数值解法有多种，比如欧拉法、龙格库塔法、线性多步法等等，其中龙格库塔法是这几种方法中比较常用的，因为其计算精度较高且具有自启动的特点。 对于用龙格库塔法求解单个常微分方程和求解常微分方程组的思路基本相似（注意一点一个微分方程组是常微分方程组即表明微分方程中的各阶导数都是对同一个变量求导，例如可以把各个量对时间求导得到一个常微分方程组，如果一个微分方程组中的有对不同变量的导数那么这个方程组就成了偏微分方程组），都是根据泰勒展开得到其迭代计算形式，其基本思想都是按照一定原则求取当前点附近一些点的斜率，通过这些斜率的线性组合作为当前点处的斜率，进行递推求解。 2、数学模型 2.1 常微分方程初值问题的数学模型 000 (,) ()()dy f x y x x dx y x y ?=>???=? （1） 式中(,) f x y 为,x y 的已知函数，0y 为给定的初始值。 常微分方程的数值解法的任务就是要求出函数值y 在微分自变量x 取如下序列： 012n x x x x <<< 时的值() (1,2,3,)i y x i n = 的近似值(1,2,3,)i y i n = ，一般情况下都采取等距结点的方式，即：0 (1,2,)i x x ih i n =+= ，其中h 为相邻两结点的距离，称为步长。 2.2 常微分方程组初值问题的数学模型 解微分方程的欧拉法，龙格-库塔法简单实例比较 欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER 法、后退EULER 法、改进的EULER 法。 缺点： 欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。 改进欧拉格式（向前欧拉公式）： 为提高精度，需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 算法： 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。对于常微分方程： (,)dy f x y dx = [,]x a b ∈ 0()y a y = 可以将区间[,]a b 分成n 段，那么方程在第i x 点有'()(,())i i i y x f x y x =，再用向前差商近似代替导数则为： ((1)()) (,())i i i i y x y x f x y x h +-= 在这里，h 是步长，即相邻两个结点间的距离。因此可以根据i x 点和i y 的数值计算出1i y +来： 1(,)i i i i y y h f x y +=+?0,1,2,i L = 这就是向前欧拉公式。 改进的欧拉公式： 将向前欧拉公式中的导数(,)i i f x y 改为微元两端导数的平均，即上式便是梯形的欧拉公式。 可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式： 数值分析中，龙格－库塔法（Runge-Kutta ）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为(,)n n f x y ,而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想： 在区间1[,]n n x x +内多取几个点，将他们的斜率加权平均，作为导数的近似。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 '(,)y f t y =00()y t y = 则，对于该问题的RK4由如下方程给出： 11234(22)6 n n h y y k k k k +=++++ 其中 1(,)n n k f t y = 21(,)22 n n h h k f t y k =++ 32(,)22 n n h h k f t y k =++ 43(,)n n k f t h y hk =++ 陕西科技大学 数值计算课程设计任务书 理学院信息与计算科学/应用数学专业信息08/数学08 班级学生： 题目：典型数值算法的C++语言程序设计 课程设计从2010 年 5 月17日起到2010 年 6 月18 日 1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）： 每人需作10个算法的程序、必做6题、自选4题。 对每个算法要求用C++语言进行编程。 必选题： 1、经典四阶龙格库塔法解一阶微分方程组 2、高斯列主元法解线性方程组 3、牛顿法解非线性方程组 4、龙贝格求积分算法 5、三次样条插值算法（压紧样条）用C++语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。 6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。 自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次. 2、对课程设计成果的要求〔包括图表、实物等硬件要求〕： 1）提交课程设计报告 按照算法要求，用C++语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档(包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程序代码为附录构成的课程设计报告。 2）课程设计报告版式要求 打印版面要求：A4纸，页边距：上2cm，下2cm，左2.5cm、右2cm；字体：正文宋体、小四号；行距：固定值20；页眉1.5cm ，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”，页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印（段前段后各0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。换行后以小四号宋体打印正文。节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。 当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用（1）、（2）……或1）、2）……顺序表示。字体为小四号宋体。 对条文内容采用分行并叙时，其编号用（a）、（b）……或a）、b）……顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范 封面 数值计算课程设计任务书 目录 数值计算设计课程设计报告正文 设计体会及今后的改进意见 参考文献（资料） 左边缘装订 3 指导教师：日期： 教研室主任：日期： 四阶龙格库塔法解一阶二元微分方程 //dxi/dt=c*(xi-xi^3/3+yi)+K*(X-xi)+c*zi //dyi/dt=(xi-b*yi+a)/c //i=1,2,3 //X=sum(xi)/N #include 龙格—库塔方法求解微分方程初值问题 （数学1201+41262022+陈晓云） 初值问题： y x x -+=2dx dy ，10≤≤x 1)0(y = 四阶龙格-库塔公式： ()y x K n n ,f 1= ????? ? ??+=+K h y x K n h n 122f ,2 ??? ??++=K y x f K h n h n 232,2 ()K h y h x f K n n 34,++= ()K K K K y y h n 4 3211n 226++++=+ 程序： 1）建立四阶龙格-库塔函数 function [ x,y ] = nark4( dyfun,xspan,y0,h ) % dyfun 为一阶微分方程的函数；y0为初始条件；xspan 表示x 的区间；h 为区间的步长； x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+k2*2+2*k3+k4)/6; end x=x;y=y; 2）执行程序（m文件） dyfun=inline('x^2+x-y'); [x,y1]=nark4(dyfun,[0,1],1,0.1); x=0:0.1:1; Format long y2=x.^2-x+1 R4=y2-y1 [x',y1',y2',R4'] y2=dsolve('Dy=x^2+x-y','y(0)=1','x') plot(x,y1,'b*-') hold on y3=inline('x^2-x+1') fplot(y3,[0,1],'ro-') legend('R-K4','解析解') 3）执行结果 ans = X RK4近似值解析值 0 1.000000000000000 1.000000000000000 0.100000000000000 0.910000208333333 0.910000000000000 0.200000000000000 0.840000396841146 0.840000000000000 0.300000000000000 0.790000567410084 0.790000000000000 0.400000000000000 0.760000721747255 0.760000000000000 0.500000000000000 0.750000861397315 0.750000000000000 0.600000000000000 0.760000987757926 0.760000000000000 0.700000000000000 0.790001102093746 0.790000000000000 0.800000000000000 0.840001205549083 0.840000000000000 0.900000000000000 0.910001299159352 0.910000000000000 1.000000000000000 1.000001383861433 1.000000000000000MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
matlab编的4阶龙格库塔法解微分方程的程序
经典四阶龙格库塔法解一阶微分方程组
MATLAB龙格-库塔方法解微分方程
龙格库塔算法解微分方程组-c语言
龙格库塔算法解微分方程组 c语言
四阶龙格库塔法解微分方程
龙格库塔法求微分方程2
蒙特卡洛求定积及龙哥库塔解微分方程
龙格-库塔法求微分方程2
四阶龙格-库塔法解微分方程(C++)
求解常微分方程组初值问题的龙格库塔法分析及其C代码
欧拉法与龙格库塔法比较分析
1、经典四阶龙格库塔法解一阶微分方程组
四阶龙格库塔法解一阶二元微分方程
龙格库塔法求微分方程matlab