青岛科技大学
学士学位论文
（2010年～2014年）
题目：常微分方程数值解法在动力天文中的应用学院：数理学院
专业班级：信计 101 学生姓名：潘力杨学号：1006010129 指导教师：王斌
起讫日期： 2014年3月————2014年6月
常微分方程数值解法在动力天文中的应用
前言
常微分方程，作为数学分支的新发展在人类认识自然，改造自然方面发挥了巨大的力量。牛顿在研究天体力学和机械动力学时，一微分方程为工具，从而在理论上得出了行星运动的规律。后来，法国天文学家勒维烈和英国天文学家亚当斯都各自分别利用了常微分方程，演算出了当时为止的海王星位置及其运动规律。如今，常微分方程更是在许多学科领域内发挥着重要作用。小行星、人造天体的运动及其规律，导弹轨道的计算，飞机空气动力学的研究，化学反应过程问题等，这些问题都可以转化为求常微分方程的解，或研究解的性质的问题。
虽然常微分方程存在众多的研究和应用领域，但事实上，只有极少数诸如线性系数常微分方程和个别典型微分方程的解能用初等函数，特殊函数或它们的级数与积分表达。而在变系数常微分方程的解析求解上已经比较困难了，更别说一般的非线性常微分方程了。大多数情况下，常微分方程只能用近似方法求解，近似方法可大致分为两大类：用级数解法，逐次逼近法等求展开式的近似解析法；和在方程一些离散点上给出近似解得数值解法。
在动力天文的具体应用中，事实上由于所处的角度关系，我们很难宏观直接地去观测天体的具体运动状态，并直接以经典物理的力学方法去计算它的准确运动轨道及规律。由于涉及到太多的不确定因素以及运动状态的复杂性，我们只能尽量地以我们的实际观测值为参数，建立模型，运用数学的方法去研究，拟合它的运动轨迹，以此为基础才能进一步作出定性的分析。在天体运动方程的求解问题上，基本上来说也只有二体问题能给出严格解。因此为了得出满足具体观测精度的定量解，我们不得不求助于数值解法。
利用数值解法求解常微分方程初值问题，一般只需要在满足规定精度的条件下求得解在若干相应点上的近似解，或是求出解的近似表达式就行，而无需使用复杂的方法试图去拟合解的表达式，或花费大量的计算去得出形式解。由于它的理论简单，可操作性强且能适应各种复杂的问题条件，因此在数学建模领域有着独一无二的重要地位。
在过去技术尚未得到发展的时代，在运用数值解法求解天体运动方程的过程中，处于计算精度，方程稳定性和收敛性的需要，我们往往得做出大量繁杂重复的运算从而必须花费大量的时间，可靠性也得不到保证，这也反过来制约了数值解法的应用。但现如今，随着计算机的普及以及编程运算技术的成熟运用，关于数值解法的应用领域也出现了新的契机，设定步长，循环次数的程式化运算大大缩短了我们的运算时间同时也提高了运算的精度，例如只要通过Matlab就可以完成简单方程求解问题的建模、迭代、作图、误差分析等一系列的工作，大大减轻了我们的负担同时也提高了运算的效率。使得我们可以更从容地将数值解法运用到天体运动求解问题的更深邻域中去。
第一章 常微分方程数值解法介绍
常微分方程可分为线性、非线性、高阶方程与方程组等类，其中线性方程包含于非线性类中，高阶方程可化为一阶方程组，若方程组中的所有未知量看作一个向量，则方程组可写成向量形式的单个方程。因此研究一阶常微分方程的初值问题
(,)()dy
f x y a x b
dx
y a y ?=≤≤???=?， (9-1) 的数值解法具有典型性，其中方程的解为:y R R →。
只有保证问题(9-1)的解存在惟一的前提下，其数值解法的研究才有意义。由常微分方程的基本理论，我们有：
定理9.1 如果(9-1)中的),(y x f 满足条件
（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在D 上关于y 满足Lipschitz 条件，即存在常数L ，使得
y y L y x f y x f -≤-),(),(
则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
所谓数值解法，是通过常微分方程离散化而给出解在某些节点上的近似值。对于问题(9-1)，在区间[],a b 引入若干节点
012N a x x x x b =<<<
<=
1,(0,1,
,1)n n n h x x n N +=-=-称为由n x 到1+n x 的步长，当n h h =（常数）时称为定步长，否则称为变
步长。多数情况下，采用等步长，即,(1,2,,)n b a
h x a nh n N N
-==+=。记问题(9-1)的精确解为()y x , 记()n y x 的近似值为n y , 记(,y )n n f x 为n f ，),,2,1(N n y n =称为问题(9-1)的数值解。
求初值问题数值解的方法一般采用步进法，即在计算出,i y i n ≤后计算1n y +，方法分为单步法和多步法。单步法是指在计算1n y +时只利用n y ，而多步法在计算1n y +时不仅要利用n y ，还要利用前面已经计算出来的若干个,1,2,
,1n j y j l -=-。我们称要用到n y 和,1,2,
,1n j y j l -=-的多步方法为l 步方法。
不论单步法和多步法，它们都有显式方法和隐式方法之分。
显式单步法的计算公式可以写为：+1=(,,)n n n n y y h x y h φ+，此公式右端不含1n y +。
隐式单步法的计算公式可以写为：+11=(,,,)n n n n n y y h x y y h φ++，此公式右端含有1n y +，从而是1n y +的方程式，需要求解方程或者采用迭代法。
显式多步法的计算公式可以写为+1-11=(,,,
,,)n n n n n n l y y h x y y y h φ-++。
隐式多步法的计算公式可以写为+11-11=(,,,,
,,),1n n n n n n n l y y h x y y y y h k l φ+-++≥- 。
显式公式与隐式公式各有特点。显式公式的优点是使用方便，计算简单，效率高，其缺点是计算精度低，稳定性差。隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。
1.1 欧拉方法（显）
欧拉方法的建立可以通过以下三种方法： （一）用差商近似导数，用向前差商h
x y x y n n )
()(1-+代替)(n x y '，则得
)1,,1,0( ))(,()
()(1-=≈-+N n x y x f h
x y x y n n n n n
)(n x y 用其近似值n y 代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有
1(,) (0,1,,1)n n n n y y hf x y n N +=+=-L 这样，问题(9-1)的近似解可以表示为
10
0(,) (0,1,,1)
()n n n n y y hf x y n N y y x +=+=-??=?L (9-2)
按式此式由初值0y 经过N 步迭代，可逐次算出N y y y ,,21，此方程称为差分方程。式(9-2)称为求解
一阶常微分方程初值问题(9-1)的欧拉公式，也称显式欧拉公式。
需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 （二）用Taylor 多项式近似，把
1()n y x +在n x 点处Taylor 展开，取一次多项式近似，则得：
2
12
1()()()()
2!
()(,())() [,]
2!
n n n n n n n n h y x y x hy x y h
y x hf x y x y x x ξξξ++'''=++''=++∈ 设步长h 的值较小，略去余项，并以n y 代替()n y x ，便得：
1(,)n n n n y y hf x y +=+ （三）将问题(9-1)中的微分方程在区间],[1+n n x x 上两边积分，可得：
)1,,1,0( ))(,()()(1
1-==
-?
++N n dx x y x f x y x y n n
x x n n
用1+n y ,n y 分别代替)(1+n x y ,)(n x y ，若对右端积分采用取左端点的矩形公式，即：
),())(,(1
n n x x y x hf dx x y x f n n
≈?
+
同样可得到显式欧拉公式。
1.2 欧拉方法（隐）
在微分方程离散化时，用向后差商代替导数，即11()()
()n n n y x y x y x h
++-'≈
，则得到如下差分方程
11100(,) (0,1,
,1)
()
n n n n y y hf x y n N y y x +++=+=-??
=? (9-3)
此公式称为求解问题(9-1)所得的数值解称为隐式欧拉法。
隐式欧拉法与欧拉公式 (9-2)式上相似，但实际计算时却复杂得多。欧拉公式 (9-2)计算1+n y 的公式中不含1+n y ，但是隐式欧拉法计算1+n y 的公式中含有1+n y 。在求解1n y +时，1,n n y x +为已知，1n y +是方程
111(,)n n n n y y hf x y +++=+的根。一般说来，这是一个非线性方程，当不能精确求解得到1n y +的表达式时，
我们需要运用简单迭代法来求解。迭代格式为：
[0]1[1][]
111(,)
(,) (0,1,2,)
n n n n k k n n n n y y hf x y y y hf x y k +++++?=+??=+=?? 由于(,)f x y 满足Lipschitz 条件，所以
[1][][]
11111111(,)(,)k k k n n n n n n n n y y h f x y f x y hL y y +++++++++-=-≤-
由此可知，只要01hL <<，迭代法就收敛到解1n y +。 1.3
梯形公式
利用数值积分方法将微分方程离散化时，若用梯形公式计算式中右端积分，即
1
11(,())[(,())(,())]2
n n
x n n n n x h
f x y x dx f x y x f x y x +++≈+? 并用1,n n y y +代替1(),()n n y x y x +，则得计算公式
111[(,)(,)]2
n n n n n n h
y y f x y f x y +++=++ (9-4)
这就是求解初值问题(9-1)的梯形公式。
梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为
[0]1
[1][]
111(,)(,)(,) (0,1,)2
n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++?=+????=++=???? 由于函数(,)f x y 关于y 满足Lipschitz 条件，所以
[1][]
[]11111111(,)(,)22
k k k n n n n n n n n h hL y y f x y f x y y y +++++++++-=-≤- 其中L 为Lipschitz 常数。因此，当012
hL
<<时，迭代法就收敛到解1n y +。
1.4 改进的欧拉法
相比于显式欧拉法和隐式欧拉法，梯形方法提高了精度，但是梯形方法，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数(,)f x y 的值，而迭代又要反复进行若干次，计算量很大。当函数
(,)f x y 比较复杂时，这个问题会变得更加突出。为了控制计算量，我们先用欧拉公式求得一个初步的近似
值1n y +，称之为预测值，预测值1n y +的精度可能很差，再用梯形公式将它校正一次得1n y +，称为校正值。这样的预测校正系统通常称为改进的欧拉公式。即 预测：1(,)n n n n y y hf x y +=+ 校正：[]111(,)(,)2
n n n n n n h
y y f x y f x y +++=+
+ 为了便于编制程序上机，将上式改写成
y 1(,)(,)1()2p n n n q n n p n p q y y hf x y y y hf x h y y y y +?
=+???
=++??
?=+??
(9-5)
算法实现如表
算法9.5
① 输入 ,,(,)a b f x y ，h ，初值0y ② 0,0,,b a
N n x a y y h
-=
=== ③ 计算
(,)(,)
p q p y y hf x y y y hf x h y =+=++
1
(),2
p q y y y x h x +?+? 输出(,)x y
④ 若1n N <-，置1n n +?，转③；否则停机。
1.5 龙格－库塔（Runge-Kutta ）法
Runge-Kutta 法由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。设初值问题(9-1)的解
1()[,]y y x C a b =∈，由微分中值定理，我们知道，必存在1[,]n n x x ξ+∈，使
1*
()()()
(,())n n n n y x y x hy y hf y y hK ξξξ+'=+=+=+ 其中*
K 叫做()y x 在1[,]n n x x +上的平均斜率。对于平均斜率*
K ，只要提供一种计算方法，该公式就给出一种数值解公式。例如，用1(,)n n K f x y =计算*
K ，就得到一阶精度的欧拉公式；用211(,)n n K f x y ++=替代*
K ，就得到隐式欧拉公式；如果用12,K K 的算术平均值计算*
K ，则可得到二阶精度的梯形公式。由此可以设想，如果在1[,]n n x x +上能多预测几个点的斜率值，用它们的加权平均值来计算*
K ，就有望得到具有较高精度的数值公式，这就是Runge-Kutta 法的基本思想。
Runge-Kutta 公式的一般形式是
1
11(,)(1,,)s
i n i n i j j j s
n n i i
i K f x
c h y h a K i s y y h b K
=+=?
=++??=??=+??
∑∑ (9-6)
其中h 为步长，s 称为方法的级数，i K 是()y y x =在 (01)n i i x c h c +≤≤点的斜率预测值，,,i i ij c b a 均
为常数。这些常数的选取原则是使公式(9-13)具有尽可能高的精度。公式(9-13)叫做s 级Runge-Kutta 公式，简称RK 方法。
Runge-Kutta 公式根据系数ij a 选择的不同，可以分为显式方法和隐式方法。当=0ij a i j ≤，时，方法为显式方法，否则为隐式方法。显式Runge-Kutta 方法可用下面的形式表示
1
1
11(,)(1,,)i i n i n ij j j s
n n i i
i K f x c h y h a K i s y y h b K
-=+=?
=++??=??=+??
∑∑ (9-7)
2级（2s =）显式RK 方法的公式为
1222111
1122(,)(,)()
n n n n n n K f x y K f x c h y ha K y y h b K b K +=??
=++??=++? （9-8）
这里12221,,,b b c a 为待定常数。 常见而二阶公式，俗称中点公式为：
1211
2(,)
(0.5,0.5)n n n n n n K f x y K f x h y hK y y hK
+=??
=++??=+? （9-9）
类似地，对3s =和4s =的显式RK 公式，通过更复杂的计算，可以导出三阶和四阶RK 公式，其中常用的三阶和四阶RK 公式分别为
1213121123(,)
(,)22
(,2)(4)6n n n n n n n n
K f x y h h K f x y K K f x h y hK hK h y y K K K +=??
?=++??=+-+??=+++??
(9-10)
和
121324311234(,)
(,)
22(,)22(,)(22)6n n n n n n
n n n n K f x y h h K f x y K h h K f x y K K f x h y hK h y y K K K K +=??
?=++?
??=++??
=++??
?=++++??
(9-11) 式(9-22)称为经典的四阶RK 公式，通常说四阶RK 方法就是指用式(9-22)求解，具体算法实现如表 9.3。
表 9.3
算法9.2
② 输入 ,,(,)a b f x y ，h ，初值0y ② 0,0,b a
N n x a h
-=
== ③ 计算
10020013002400301234000
(,)
(,)
22(,)
22
(,)(22)
6
,K f x y h h
K f x y K h h
K f x y K K f x h y hK h
y y K K K K x h x y y ==++=++=++=+++++??
输出00(,)x y
④ 若1n N <-，则置1,n n +?转 ③；否则停机。
1.6 线性多步法
之前的各种数值解法，都是单步法，即在每一步计算时，只要前面一个值n y 已知的条件下就可以计算出1n y +。单步法的特点是，可以自成系统进行直接计算，因为初始条件只有一个已知0y ，由0y 可以计算1y ，不必借助于其它方法，这种单步法是自开始的。
如果考虑在计算值1n y +时，能够比较充分地利用前面的已知信息，如n y 和,1,2,,1n j y j l -=-，那么就可望使所得到的1n y +更加精确。这就是多步法的基本思想。我们称要用到n y 和,1,2,,1n j y j l -=-的多
步法为l 步方法。
多步法中最常用的是线性多步法，其一般公式是 10
1
r
r
n i n i
i n i i i y y
h f αβ+--==-=
+∑∑ (9-12)
其中,k k αβ均为常数，(),(,) (1,,)n k n k k k k y y x f f x y k n n r --===+-。当10β-=时，上式称为
显式，否则称为隐式。
定义：设()y x 是初值问题(9-1)的精确解，线性多步法(9-12)在1n x +上的局部截断误差定义为： 110
1
()r
r
n n k
n k k n
k k k T y x y h y α
β++--==-'=-+∑∑ 若1
1()p n T O h ++=，则称方法(9-12)是p 阶方法。
线性多步法的构造原理：利用Taylor 展开导出线性多步公式(9-12)的基本方法是将线性多步公式(9-12)在n x 处进行Taylor 展开，然后与1()n y x +在n x 处的Taylor 展开式相比较，要求它们前面的项重合，由此确
定参数,k k αβ。设初值问题(9-1)的解()y x 充分光滑，记()
()() (0,1,
)j j n n y y x j ==，把它在n x 处展开成
Taylor 级数，则有：
1()(1)
1
()()()!(1)!j p p
j p n n n n n j x x x x y x y y y j p ++=--=+++
+
∑
及
1()(1)
1
()()()(1)!!j p p
j p n n n n j x x x x y x y y j p -+=--'=++
-
∑
用0()n k x x n k h -=+-替代式以上两式中的x ，则得： 1()(1)
1
()()()!(1)!j p p
j p n k
n k n n n j kh kh y y x y y y j p ++--=--==+++
+
∑
把这两个式子代入式(9-12)，得：
1()(1)
1011()(1)111()01011()()!(1)!()()(1)!!()()!()(1)!j p p
r
j p n k n n n k j j p p r
j p k n n k j j p r r r j
j j k n k k n k j k k p p k kh kh y y y y j p kh kh h y y j p h y k j k y j h k p αβααβα+++==-+=-=-====-++??
--=+++
?+?
?
??
--+++ ?
-??
??=+-+-????
+-+∑∑∑∑∑∑∑∑1(1)11(1)()r r
p p k n k k p k y β+==-??++-+????
∑
∑
为了使式(9-12)具有p 阶精度，只须使上式的前1p +项与1()n y x +在n x 的Taylor 展开式
1()(1)
11()!
(1)!j p p
j p n n n n j h h y x y y y j p +++==++++
∑
的前1p +项系数对应相等即可。对比关于h 的同次项系数，得到确定,k k αβ的方程组：
11
11 (1,,)()()1r
k k r
r
j j k k k k j p k j k ααβ=-==-?=??=??-+-=??∑∑∑ （9-13）
可见，只要式(9-12)的系数,k k αβ满足上式，方法(9-12)就具有p 阶精度。由此可得局部截断误差为：
11(1)21111()(1)()()(1)!p r r
p p
p p n k k n k k h T k p k y O h p αβ+++++==-??=---+-+??+??
∑∑ 式(9-13)共有1p +个方程，23r +个待定常数：01101,,,,,,,r r αααββββ-，只要231r p +≥+，
即12
p
r ≥-，就可以由式(9-13)，定出具有p 阶精度的公式(9-12)的系数,k k αβ.
常用的线性多部法公式有：
一、阿达姆斯（Adams ）公式1123(5559379)24n n n n n n h
y y f f f f +---=+-+- 二、米尔恩（Milne ）公式 13124(22)3n n n n n h
y y f f f +---=+-+ 三、海明（Hamming ）公式121113
(9)(2)88
n n n n n n y y y h f f f +-+-=-++-
其中 阿达姆斯（Adams ）公式的构造原理为： 取4,3p r ==，有：
3
3
3
10
11 (1,2,3,4)()()1k k j j k k k k j k j k ααβ=-==-?=??=??-+-=??∑∑∑ (9-14) 此时方程有9个未知数，5个方程，有四个自由变量。特别取12310αααβ-====可解得：
1012355593791,,,,24242424
αββββ===-==- 相应的线性多步法公式为： 1123(5559379)24
n n n n n n h
y y f f f f +---=+
-+- 因为10β-=，上式称为Adams 显式公式，它是四阶公式，局部截断误差为：
533
5
4(5)61
115(5)
61()5()()5!251()720
n k k n k k n h T k k y O h h y O h αβ+==-??=----+????
=
+∑∑
如果令12330αααβ====，由式(9-14)可得解
010********,,,,24242424
αββββ-====-= 相应的线性多步法公式为: 1112(9195)24
n n n n n n h
y y f f f f ++--=+
+-+ (9-15) 因为10β-≠,式(9-15)称为Adams 隐式公式，它是四阶公式，局部截断误差为: 5(5)
6119()720
n n T h y O h +=-
+ 由于线性多步法公式只有在知道了前面的1,,,n n n r y y y --之后，才能使用。所以开头的12,,,r
y y y 的值必须先用同阶的单步法（例如Runge-Kutta 法）求出，然后才能用线性多步法公式。
第二章 常微分方程数值解法的程序实现
对于数值求解微分方程（组），MATLAB 有专门的求解器可以调用。MATLAB 中求微分方程的数值解命令的调用格式为：
[X,Y] = solver(odefun,tspan ,0y )
其中solver 为MATLAB 中求解微分方程的求解器：例如ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 等；odefun 是要求解的显式常微分方程：
(,)dy
f x y dx
，；0tspan=[,]end x x 为积分区间；0y 为初始条件。输出的X 为0tspan=[,]end x x 区间上点的列向量，Y 为对应于X 上各个点的数值解所组成的向量。要
获得问题在其他指定时间点012,,,end x x x x ，上的解，则令012tspan=[,,,,]end x x x x （要求是单调的）
。 因为没有一种算法可以有效地解决所有的求解常微分方程问题，为此，MATLAB 提供了多种求解器
Solver ，对于不同的常微分方程问题，采用不同的Solver ，一些常用的Solver 总结如下：
要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 MA TLAB 的常用程序，其中：
ode23 采用Runge-Kutta 2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度； ode45 则采用Runge-Kutta 4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度。
3.1欧拉法的Matlab 程序实现
隐式方法的程序去掉，整篇论文只给显式方法的程序
③改进欧拉法
Function [x,y]=naeuler2(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1)=y(0);
For n=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
y(n+1)=y(n)+h*k1;
k2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(k1+k2)/2;
End
x=x’;y=y’;
x1=0:0.2:1;y1=(1+2*x1).^0.5;
Plot(x,y,x1,y1)
2.2龙格库塔法的Matlab程序实现
三阶龙格—库塔公式：
function y = DELGKT3_kuta(f, h,a,b,y0,varvec)
format long;
N = (b-a)/h;
y = zeros(N+1,1);
y(1) = y0;
x = a:h:b;
var = findsym(f);
for i=2:N+1
K1 = Funval(f,varvec,[x(i-1) y(i-1)]);
K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);
K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);
y(i) = y(i-1)+h*(K1+4*K2+K3)/6;
end
format short;
DELGKT3_kuta
函数运行时需要调用下列函数：
function fv=Funval(f, varvec, varval)
var= findsym(f);
If length(var)<4
if var(1)==varvec(1)
fv=subs(f,varvec(1),varval(1));
else
fv=subs(f,varvec(2),varval(2));
end
else
fv=subs(f,varvec,varval);
End
四阶龙格—库塔法的计算公式为：
function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec)
format long;
N = (b-a)/h;
y = zeros(N+1,1);
y(1) = y0;
x = a:h:b;
var = findsym(f);
for i=2:N+1
K1 = Funval(f,varvec,[x(i-1) y(i-1)]);
K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);
K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]);
K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);
y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;
end
format short; 同理也要先调用Funval函数（见上）
2.3adams算法的matlab程序实现
选择Adams算法中的一个显式方法，给出matlab程序实现即可
第三章动力天文中的经典问题的数值求解
3.2行星运动模型
二体问题（two-body problem）是指研究两个可以视为质点的天体在其相互之间的万有引力作用下的动力学问题。二体问题作为天体力学中的一个最基本的近似模型，是各类天体真实运动的第一次近似结果，是研究天体精确运动的理论基础，也是天体力学中的一个基本问题，因此它具有很重要的意义。常见的应用有卫星绕着行星公转、行星绕着恒星公转、双星系统、双行星、一个经典电子绕着原子核运动等。
设太阳中心所在位置为复平面之原点O，在时刻t，行星位于
θi
)(（4.1）
Z=
t
re
所表示的点P。这里均是t的函数，分别表示的模和辐角。于是行星速度为：
??? ??+=+=dt d ir dt
dr
e dt d ire e dt dr dt dZ i i i θθθθθ 其加速度为：
???
????????? ??++???? ????? ??-=dt d dt dr dt r d r i dt d r dt r d e dt Z d i θθθ22222222 （4.2） 而太阳对行星的引力依万有引力定律，大小为
2
r mMG
方向由行星位置P 指向太阳的中心O ，故为 θi e r
mMG 2
-
，其中)(10989.130kg M ?=为太阳的质量,m 是行星的质量，)/(10672.62
211kg m N G ??=- 为万有引力常数。依Newton 第二定律，得到
222
dt
Z
d m
e r mMG i =-θ （4.3） 将（4.2）代入（4.3），然后比较实部和虚部，有
???????-=??? ??-=+2222
2
202r MG dt d r dt
r d dt d dt dr dt
d r θθθ 如设0=t 时，行星正处于远日点，远日点位于正实轴上，距原点O 为0r ,行星的线速度为0v
,那么就
有初值条件：
?????????
??========000000000r v dt
d dt dr r r t t t t θθ
方程(4.4) ~ (4.9)就是行星绕太阳运行的轨迹的数学模型。
将式（4.4）乘以r ，即得
02=??
? ??dt d r dt d θ 常数）( 12
C dt
d r
=?θ
，其中001v r C =。 这样，有向线段OP 在时间t ?内扫过的面积等于2
2112t C dt dt d r t
t t
?=?
?+θ
，而这正是Kepler 第二定律。将式(4.10)改写后代入式(4.5)
2
32
122
r MG
r dt r d C -
=- 由此得到行星运动的较为简单形式的数学模型：??????
?
?
??
?????====-=-===00000021
2
32
1
22t t t dt dr r r r C dt d r MG r dt
r d C θθ
为了计算二体运动，我们选择其中一个天体作为坐标系的中心，则另一个天体的运动轨迹为一圆锥曲线。我们用二维12(q ,q )表示第二个天体的位置，根据上述分析，由牛顿定律，经过一定的标准化化简，满足下面的微分方程
211
112223/2
12
2
22222223/212
,(0)1,(0)0(+),(0)0,(0)(+)d q q q e q dt q q d q q q q dt q q ?'=-=-=????'=-==??
用第三章给出程序的算法求解此问题，其中e 取为0.6，步长和求解区间你自己取，最后画出每种方法的数值解12(q ,q )。给出求解此问题的matlab 程序。 3.2
人造卫星运动模型
把这部分内容翻译完写进论文，利用第三章给出程序的算法求解此问题，得到数值解，步长和求解区间你自己取，最后画出每种方法的数值解12(q ,q )，给出求解此问题的matlab 程序。
3.3 带有扰动的轨道问题
这部分内容翻译完写进论文，利用第三章给出程序的算法求解此问题，得到数值解，步长你自己取，最后(q,q)，给出求解此问题的matlab程序。
画出每种方法的数值解
12
第四章总结
多写些参考文献
第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)；
(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<
高阶线性微分方程常用解法简介 关键词：高阶线性微分方程 求解方法 在微分方程的理论中，线性微分方程是非常值得重视的一部分内容，这不仅 因为线性微分方程的一般理论已被研究的十分清楚，而且线性微分方程是研究非线性微分方程的基础，它在物理、力学和工程技术、自然科学中也有着广泛应用。下面对高阶线性微分方程解法做一些简单介绍. 讨论如下n 阶线性微分方程：1111()()()()n n n n n n d x d x dx a t a t a t x f t dt dt dt ---++++= （1），其中()i a t （i=1,2,3,,n ）及f(t)都是区间a t b ≤≤上的连续函数，如果 ()0f t ≡，则方程（1）变为 1111()()()0n n n n n n d x d x dx a t a t a t x dt dt dt ---++++= （2），称为n 阶齐次线性微分方程，而称一般方程（1）为n 阶非齐次线性微分方程，简称非齐次线性微分方程，并且把方程（2）叫做对应于方程（1）的齐次线性微分方程. 1.欧拉待定指数函数法 此方法又叫特征根法，用于求常系数齐次线性微分方程的基本解组。形如 111121[]0,(3),n n n n n n n d x d x dx L x a a a x dt dt dt ---≡++++=其中a a a 为常数，称为n 阶常系数齐次线性微分方程。 111111111111[]()()()n t n t t t t n n n n n n n t t n n n n n n n d e d e de L e a a a e dt dt dt a a a e F e F a a a n λλλλλλλλλλλλλλλλ---------≡++++=++++≡≡++++其中=0(4)是的次多项式. ()F λ为特征方程，它的根为特征根. 1.1特征根是单根的情形 设12,,,n λλλ是特征方程111()0n n n n F a a a λλλλ--≡++++=的n 个彼此不相等的根，则应相应地方程（3）有如下n 个解：12,,,.n t t t e e e λλλ（5）我们指出这n 个解在区间a t b ≤≤上线性无关，从而组成方程的基本解组. 如果(1,2,,)i i n λ=均为实数，则（5）是方程（3）的n 个线性无关的实值 解，而方程（3）的通解可表示为1212,n t t t n x c e c e c e λλλ=+++其中12,,,n c c c 为任意常数. 如果特征方程有复根，则因方程的系数是实常数，复根将称对共轭的出现.设1i λαβ=+是一特征根，则2i λαβ=-也是特征根，因而于这对共轭复根
第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(y g x f dx dy = 当0)(≠y g 时，得到 dx x f y g dy )() (=，两边积分即可得到结果； 当0)(0=ηg 时，则0)(η=x y 也是方程的解。 例1.1、 xy dx dy = 解：当0≠y 时，有 xdx y dy =，两边积分得到)(2ln 2为常数C C x y += 所以)(112 12 C x e C C e C y ±==为非零常数且 0=y 显然是原方程的解； 综上所述，原方程的解为)(12 12 为常数C e C y x = ②、形如0)()()()(=+dy y Q x P dx y N x M 当0)()(≠y N x P 时，可有 dy y N y Q dx x P x M ) () ()()(=，两边积分可得结果； 当0)(0=y N 时，0y y =为原方程的解，当0(0=） x P 时，0x x =为原方程的解。 例1.2、0)1()1(2 2 =-+-dy x y dx y x 解：当0)1)(1(2 2 ≠--y x 时，有 dx x x dy y y 1 122-=-两边积分得到 )0(ln 1ln 1ln 22≠=-+-C C y x ，所以有)0()1)(1(22≠=--C C y x ； 当0)1)(1(2 2 =--y x 时，也是原方程的解； 综上所述，原方程的解为)()1)(1(2 2 为常数C C y x =--。 ⑵可化为变量可分离方程的方程： ①、形如 )(x y g dx dy = 解法：令x y u =，则udx xdu dy +=，代入得到)(u g u dx du x =+为变量可分离方程，得到
常微分方程解题方法总结 来源：文都教育 复习过半， 课本上的知识点相信大部分考生已经学习过一遍 . 接下来， 如何将零散的知 识点有机地结合起来， 而不容易遗忘是大多数考生面临的问题 . 为了加强记忆， 使知识自成 体系，建议将知识点进行分类系统总结 . 著名数学家华罗庚的读书方法值得借鉴， 他强调读 书要“由薄到厚、由厚到薄”，对同学们的复习尤为重要 . 以常微分方程为例， 本部分内容涉及可分离变量、 一阶齐次、 一阶非齐次、 全微分方程、 高阶线性微分方程等内容， 在看完这部分内容会发现要掌握的解题方法太多， 遇到具体的题 目不知该如何下手， 这种情况往往是因为没有很好地总结和归纳解题方法 . 下面以表格的形 式将常微分方程中的解题方法加以总结，一目了然，便于记忆和查询 . 常微分方程 通解公式或解法 ( 名称、形式 ) 当 g( y) 0 时，得到 dy f (x)dx ， g( y) 可分离变量的方程 dy f ( x) g( y) 两边积分即可得到结果； dx 当 g( 0 ) 0 时，则 y( x) 0 也是方程的 解 . 解法：令 u y xdu udx ，代入 ，则 dy 齐次微分方程 dy g( y ) x dx x u g (u) 化为可分离变量方程 得到 x du dx 一 阶 线 性 微 分 方 程 P ( x)dx P ( x) dx dy Q(x) y ( e Q( x)dx C )e P( x) y dx
伯努利方程 解法：令 u y1 n，有 du (1 n) y n dy ， dy P( x) y Q( x) y n（n≠0,1）代入得到du (1 n) P(x)u (1 n)Q(x) dx dx 求解特征方程：2 pq 三种情况： 二阶常系数齐次线性微分方程 y p x y q x y0 二阶常系数非齐次线性微分方程 y p x y q x y f ( x) (1)两个不等实根：1, 2 通解： y c1 e 1x c2 e 2x (2) 两个相等实根：1 2 通解： y c1 c2 x e x (3) 一对共轭复根：i , 通解： y e x c1 cos x c2 sin x 通解为 y p x y q x y 0 的通解与 y p x y q x y f ( x) 的特解之和. 常见的 f (x) 有两种情况： x （ 1）f ( x)e P m ( x) 若不是特征方程的根，令特解 y Q m ( x)e x；若是特征方程的单根，令特 解 y xQ m ( x)e x；若是特征方程的重根， 令特解 y*x2Q m (x)e x； （2）f (x) e x[ P m ( x) cos x p n ( x)sin x]
第三章一阶微分方程解的存在定理 [教学目标] 1.理解解的存在唯一性定理的条件、结论及证明思路，掌握逐次逼近法，熟练 近似解的误差估计式。 2.了解解的延拓定理及延拓条件。 3.理解解对初值的连续性、可微性定理的条件和结论。 [教学重难点] 解的存在唯一性定理的证明，解对初值的连续性、可微性定理的 证明。 [教学方法] 讲授，实践。 [教学时间] 12学时 [教学内容] 解的存在唯一性定理的条件、结论及证明思路，解的延拓概念及延 拓条件，解对初值的连续性、可微性定理及其证明。 [考核目标] 1.理解解的存在唯一性定理的条件、结论，能用逐次逼近法解简单的问题。 2.熟练近似解的误差估计式，解对初值的连续性及可微性公式。 3.利用解的存在唯一性定理、解的延拓定理及延拓条件能证明有关方程的某些性质。 §1 解的存在性唯一性定理和逐步逼近法 微分方程来源于生产实践际，研究微分方程的目的就在于掌握它所反映的客 观规律，能动解释所出现的各种现象并预测未来的可能情况。在第二章介绍了一 阶微分方程初等解法的几种类型，但是，大量的一阶方程一般是不能用初等解法 求出其通解。而实际问题中所需要的往往是要求满足某种初始条件的解。因此初 值问题的研究就显得十分重要，从前面我们也了解到初值问题的解不一定是唯一的。他必须满足一定的条件才能保证初值问题解的存在性与唯一性，而讨论初值 问题解的存在性与唯一性在常微分方程占有很重要的地位，是近代常微分方程定 性理论，稳定性理论以及其他理论的基础。 例如方程 过点(0,0)的解就是不唯一，易知0 y=是方程过(0,0)的解，此外，容易验证，2 =或更一般地，函数 y x 都是方程过点(0,0)而且定义在区间01 <<的任一数。 c ≤≤上的解，其中c是满足01 x
淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析 学院、专业数学科学学院数学与应用数学 研究方向计算数学 学生姓名李娜 学号 20091101070 指导教师姓名陈昊 指导教师职称讲师 年月日
常微分方程数值解法的误差分析 李娜 （淮北师范大学数学科学学院，淮北，235000） 摘要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler 法、后退Euler法、梯形方法、Runge—Kutta方法、投影法和多步法，等等．本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。 关键词：常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差
Error Analysis of Numerical Method for Solving the Ordinary Differential Equation Li Na (School of Mathematical Science, Huaibei Normal University, Huaibei, 235000) Abstract In nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution,there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on.Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential. Keywords:Ordinary differential equations, numerical solution methods, s ingle ste p methods, l inear multi-step methods, local truncation error
i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
一阶常微分方程的解法 摘要：常微分方程是微积分学的重要组成部分，广泛用于具体问题的研究中，在整个数学中占有重要的地位。本文对一阶常微分方程的解法作了简要的总结，并举例加以分析了变量可分离方程，线性微分方程，积分因子，恰当微分方程，主要归纳了一阶微分方程的初等解法，并以典型例题加以说明。 关键词：变量分离;积分因子;非齐次微分方程;常数变易法 Solution of first-order differential equation Abstract: Differential equations, important parts of calculus, are widely used in the research of practical problems, which also play important role in mathematics. The solution of a differential equation is summarized briefly, and illustrates the analysis of variable separable equation, linear differential equation, integral factor, exact differential equation, mainly summarizes the elementary solution of first order differential equations, and the typical examples to illustrate. Keywords: variable separation; integral factor; non-homogeneous differential equation; constant variation method 1. 引言 一阶常微分方程初等解法,就是把常微分方程的求解问题转化为积分问题, 能用这种方法求解的微分方程称为可积方程. 本文通过对一阶微分方程的初等解法的归纳与总结，以及对变量分离，积分因子，微分方程等各类初等解法的简要分析，同时结合例题把常微分方程的求解问题化为积分问题，进行求解. 2. 一般变量分离 2.1 变量可分离方程 形如 ()()dy f x g y dx = （1.1） 或 1122()()()()M x N y dx M x N y dy = （1.2） 的方程，称为变量可分离方程。分别称（1.1）、（1.2）为显式变量可分离方程和 微分形式变量可分离方程[1] . (1) 显式变量可分离方程的解法 在方程（1.1）中， 若()0g y ≠，（1.1）变形为 ()() dy f x dx g y =
常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate
各类微分方程的解法 1.可分离变量的微分方程解法 一般形式:g(y)dy=f(x)dx 直接解得∫g(y)dy=∫f(x)dx 设g(y)及f(x)的原函数依次为G(y)及F(x),则G(y)=F(x)+C为微分方程的隐式通解 2.齐次方程解法 一般形式:dy/dx=φ(y/x) 令u=y/x则y=xu,dy/dx=u+xdu/dx,所以u+xdu/dx=φ(u),即du/［φ(u)-u］=dx/x 两端积分,得∫du/［φ(u)-u］=∫dx/x 最后用y/x代替u,便得所给齐次方程的通解 3.一阶线性微分方程解法 一般形式:dy/dx+P(x)y=Q(x) 先令Q(x)=0则dy/dx+P(x)y=0解得y=Ce－ ∫P(x)dx,再令y=u e－∫P(x)dx代入原方程解得u=∫Q(x) e∫P(x)dx dx+C,所以y=e－∫P(x)dx［∫Q(x)e∫P(x)dx dx+C］ 即y=Ce－∫P(x)dx +e－ ∫P(x)dx∫Q(x)e∫P(x)dx dx为一阶线性微分方程的通解 4.可降阶的高阶微分方程解法 ①y(n)=f(x)型的微分方程 y(n)=f(x) y(n-1)= ∫f(x)dx+C1 y(n-2)= ∫［∫f(x)dx+C1］dx+C2 依次类推,接连积分n次,便得方程y(n)=f(x)的含有n个任意常数的通解②y”=f(x,y’) 型的微分方程 令y’=p则y”=p’,所以p’=f(x,p),再求解得p=φ(x,C1) 即dy/dx=φ(x,C1),所以y=∫φ(x,C1)dx+C2 ③y”=f(y,y’) 型的微分方程
令y ’=p 则y ”=pdp/dy,所以pdp/dy=f(y,p),再求解得p=φ(y,C 1) 即dy/dx=φ(y,C 1),即dy/φ(y,C 1)=dx,所以∫dy/φ(y,C 1)=x+C 2 5.二阶常系数齐次线性微分方程解法 一般形式:y ”+py ’+qy=0,特征方程r 2+pr+q=0 6.二阶常系数非齐次线性微分方程解法 一般形式: y ”+py ’+qy=f(x) 先求y ”+py ’+qy=0的通解y 0(x),再求y ”+py ’+qy=f(x)的一个特解y*(x) 则y(x)=y 0(x)+y*(x)即为微分方程y ”+py ’+qy=f(x)的通解 求y ”+py ’+qy=f(x)特解的方法: ① f(x)=P m (x)e λx 型 令y*=x k Q m (x)e λx ［k 按λ不是特征方程的根,是特征方程的单根或特征方程的重根依次取0,1或2］再代入原方程,确定Q m (x)的m+1个系数 ② f(x)=e λx ［P ｌ(x)cos ωx+P n (x)sin ωx ］型 令y*=x k e λx ［Q m (x)cos ωx+R m (x)sin ωx ］［m=max ﹛ｌ,n ﹜,k 按λ+i ω不是特征方程的根或是特征方程的单根依次取0或1］再代入原方程,分别确定Q m (x)和R m (x)的m+1个系数
第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+
相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ 第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果： (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理：若一步方法满足： (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1．局部截断误差 局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意：y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。 第八章 8.4讲 第四节 二阶常系数线性微分方程 一、二阶常系数线形微分方程的概念 形如 )(x f qy y p y =+'+'' (1) 的方程称为二阶常系数线性微分方程.其中p 、q 均为实数,)(x f 为已知的连续函数. 如果0)(≡x f ,则方程式 (1)变成 0=+'+''qy y p y (2) 我们把方程(2)叫做二阶常系数齐次线性方程,把方程式(1)叫做二阶常系数非齐次线性方程. 本节我们将讨论其解法. 二、二阶常系数齐次线性微分方程 1．解的叠加性 定理1 如果函数1y 与2y 是式(2)的两个解, 则2211y C y C y +=也是式(2)的解,其中21,C C 是任意常数. 证明 因为1y 与2y 是方程(2)的解,所以有 0111 =+'+''qy y p y 0222 =+'+''qy y p y 将2211y C y C y +=代入方程(2)的左边,得 )()()(22112211221 1y C y C q y C y C p y C y C ++'+'+''+'' =0)()(2222111 1=+'+''++'+''qy y p y C qy y p y C 所以2211y C y C y +=是方程(2)的解. 定理1说明齐次线性方程的解具有叠加性. 叠加起来的解从形式看含有21,C C 两个任意常数,但它不一定是方程式(2) 的通解. 2.线性相关、线性无关的概念 设,,,,21n y y y 为定义在区间I 内的n 个函数,若存在不全为零的常数 ,,,,21n k k k 使得当在该区间内有02211≡+++n n y k y k y k , 则称这n 个函数在区间I 内线性相关,否则称线性无关. 例如 x x 2 2 sin ,cos ,1在实数范围内是线性相关的,因为 0sin cos 12 2≡--x x 又如2,,1x x 在任何区间(a,b)内是线性无关的,因为在该区间内要使 02321≡++x k x k k 必须0321===k k k . 对两个函数的情形,若 =21y y 常数, 则1y ,2y 线性相关,若≠2 1y y 常数, 则1y ,2y 线性无关. 3.二阶常系数齐次微分方程的解法 定理 2 如果1y 与2y 是方程式(2)的两个线性无关的特解,则 212211,(C C y C y C y +=为任意常数)是方程式(2)的通解. 例如, 0=+''y y 是二阶齐次线性方程,x y x y cos ,sin 21==是它的两个解,且 ≠=x y y tan 2 1 常数,即1y ,2y 线性无关, 所以 x C x C y C y C y cos sin 212211+=+= ( 21,C C 是任意常数)是方程0=+''y y 的通解. 由于指数函数rx e y =(r 为常数)和它的各阶导数都只差一个常数因子, 常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步： 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。基本模型 1. 发射卫星为什么用三级火箭 2. 人口模型 3. 战争模型 4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1. 改进Euler 法： 2. 龙格—库塔( Runge—Kutta )方法： 【源程序】 1. 改进Euler 法： function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子 区间个数 if nargin<5,n=5O;end h=(x1-xO)/n; x(1)=xO;y(1)=yO; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command 窗口 f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O) 2 x +2x , (0 < x < , y(0) = 1 求解函数y'=-2y+2 2. 龙格—库塔( Runge—Kutta )方法： [t,y]=solver('F',tspan ，y0) 这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。 注：ode45和ode23变步长的，采用Runge-Kutta算法。 ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长(变步长)的常微分方程数值解法，其整体截断误差为(△ 口人5解 决的是Nonstiff(非刚性)常微分方程。 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 郑州航空工业管理学院 毕业论文（设计） 2015届数学与应用数学专业1111062班级 题目二阶常微分方程的降阶解法 姓名贾静静学号111106213 指导教师程春蕊职称讲师 2015年4月5号 二阶常微分方程的降阶解法 摘要 常微分方程是数学领域的一个非常重要的课题，并在实践中广泛于解决问题，分析模型。常微分方程在微分理论中占据首要位置,普遍应用在工程应用、科学研究以及物理学方面,不少应用范例都归结为二阶线性常微分方程的求解问题。而正常情况下，常系数微分方程依据线性常微分方程的日常理论是可以求解的.不过对于变系数二阶线性常微分方程的求解却有一定程度的困难,迄今为止还没有一个行之有效的普遍方法。 本文主要考虑了二阶常系数线性微分方程的降阶法。关于二阶常系数线性微分方程的求解问题,首先,我们给出二阶齐次常系数线性微分方程的特征方程，并求解出特征方程的两个特征根；其次，利用积分因子乘以微分方程和导数的运算，将二阶常系数线性微分方程化为一阶微分形式；最后，将一阶微分形式两边同时积分，求解一阶线性微分方程，可求得二阶常系数线性微分方程的一个特解或通解。关于二阶变系数齐次线性微分方程的求解问题，化为恰当方程通过降阶法求解二阶齐次变系数微分方程的通解。对于非齐次线性微分方程，只需再运用常数变易法求出它的一个特解，问题也就相应地解决了。 关键词 二阶常微分方程；降阶法；特征根；常数变易法；一阶微分形式 Order reduction method of second order ordinary differential equations Jingjing Jia Chunrui Cheng 111106213 Abstract Ordinary differential equation is a very important topic in the field of mathematics, it has been widely used in solving the problem and analyzing model in practice . Ordinary differential equations in the theory of differential occupied first place, it has been widely used in engineering application and scientific research as well as physics, many application examples are attributed to second order linear ordinary differential equation solving problem. And under normal circumstances,ordinary coefficient differential equation on the basis of the linear often daily theory of differential equations is can be solved. But for the solution for variable coefficient second order linear ordinary differential equations have a certain degree of difficulty, so far we haven't a well-established general method. This paper mainly introduces the method of reduction of order two order linear differential equation with constant coefficients.On the problem of solving the linear differential equation with two order constant coefficients,first, we give homogeneous ordinary coefficient linear differential equation of the characteristic equation and solve the two characteristic roots of characteristic equation;secondly,we should use the integral factor times differential equation and derivative operation and turn two order constant常微分方程数值解法
(整理)二阶常系数线性微分方程的解法word版.
常微分方程数值解法
常微分方程的数值解
二阶常微分方程的降阶解法
相关主题
文本预览