课程设计-常微分方程的求解

常微分方程的求解问题
摘要
题目中显然涉及到求解数值积分和求解常微分方程问题.
首先，需要对常微分方程中的数值积分进行求解.回顾所学过的求解数值积分的方法，有Newton-Cotes公式、复化求积公式、Romberg算法和Gauss-Legender公式等，我们选择其中的一种方法来进行求解.在求解时，选择了包含在Newton-Cotes公式中的梯形公式进行了求解.
然后，将积分求解的结果代入原方程中，进一步再求解常微分方程.回顾所学过的求解常微分方程的方法，有Euler格式、后退的Euler格式、梯形格式、改进的Euler格式、Euler两步格式、Runge-Kutta方法、Taylor级数法等，我们选择其中的一种方法进行求解.在求解时，选择了四阶Runge-Kutta方法进行了求解.
最后，应用MATLAB软件进行了求解计算了t=1,2,..,10时的值，并对模型进行了优缺点分析及模型的推广.
关键词：数值积分；常微分方程；Newton-Cotes公式；Runge-Kutta方法；
MATLAB软件
THE SOLUTION OF ORDINRY DIFFERENTIAL EQUATIONS
ABSTRACT
Problems in apparently involves numerical integration and solving problems in ordinary differential equations.
First of all, we need to give the solution of ordinary differential equations. Take a review of the numerical integration method, a Newton-Cotes formula, complex formula, Romberg algorithm and Gauss-Legender formula, we choose one of the methods to solve the problem. In this method, we choose the Newton-Cotes formula which includes the trapezoid formula to get the solution.
Then, the integral solution should be used in the original equation, and then get the solution of ordinary differential equation. Take a review of what has been learned for solving ordinary differential equation method, are Euler format, Euler format, the format back trapezoid, improved Euler format, Euler two, Runge-Kutta method, Taylor format series method, we choose one of the methods to solve the equation. In the solution, we choose four order Runge-Kutta method to solve.
Finally, we use the MATLAB software to get the solution, and analyze the advantages and the disadvantage of the model and give a promotion of the model.
Key words: Numerical integration of ordinary differential equations; formula;
Newton-Cotes; Runge-Kutta method；MATLAB
目录
1 问题的提出 (1)
2 问题的分析 (1)
3 问题假设 (1)
4 符号说明 (2)
5 模型的建立 (2)
5.1 模型的准备工作 (2)
5.2 建立模型 (5)
6 模型求解 (6)
6.1 MATLAB软件求解结果 (7)
7 模型验证及结果分析 (7)
7.1 模型验证 (7)
7.2 结果分析 (7)
8 模型的优缺点 (8)
8.1 模型的优点 (8)
8.2 模型的缺点 (8)
9 模型的推广 (8)
参考文献 (9)
附录 (10)
腹有诗书气自华
1 问题提出
用数值分析中求解常微分方程的知识，求解下列方程：
'
20
()()t
y t ay b y s ds =+⎰ (05t ≤≤ )
其中 0,0a b >< 初值为
(0)1
y = .
2 问题分析
从方程'
20
()()t
y t ay b
y s ds =+⎰
中可以看到，方程右端中含有积分求解.故
在求解此常微方程之前，必须先求解积分
20
()t
y s ds ⎰
.
求解数值积分的方法有Newton-Cotes 公式、复化求积公式、Romberg 算法和Gauss-Legender 公式等，我们可以选择其中的一种方法进行求解.
求解常微分方程的方法有Euler 格式、后退的Euler 格式、梯形格式、改进的Euler 格式、Euler 两步格式、Runge-Kutta 方法、Taylor 级数法等，我们可以选择其中常用的方法进行求解.
3 问题假设
(1) 设应用MATLAB 进行求解时，舍入误差可以忽略，不用考虑； (2) 设变量 t 的范围设定不会影响方程的求解；
(3) 设积分求解的精度和常微分方程的求解精度不同时，对方程的求解没有影响；
腹有诗书气自华
4 符号说明
T ：设用梯形积分公式求解的积分值； S ：设用Simpson 积分公式求解的积分值； C ：设用Cotes 公式求解的积分值；
n y : 设求解常微分方程时第n 次值；
5 模型的建立
5.1 模型的准备工作
1）对积分求解的准备
设()f x 是[,]a b 上连续可积的实函数，求()f x 在[,]a b 上的数值积分也就是求()f x 在[,]a b 上的定积分的数值解.对于这个问题来说，即使能找到()f x 的一个原函数()F x 的解析形式，并利用牛顿—莱布尼茨（Newton-Leibniz ）公式进行计算，在许多情况下这样做也是非常麻烦的.问题的关键在于，初等函数的原函数不一定就是初等函数，在这种情况下并不能利用牛顿—莱布尼茨公式进行计算.那么就想到绕过牛顿—莱布尼茨公式而直接求一般的定积分的数值解.
现在所要解决的常微分方程中的积分就是如此，我们可以用Newton-Cotes 公式、复化求积公式、Romberg 算法和Gauss-Legender 公式等公式进行近似计算定积分的数值解.
下面将其中的几种常用公式简单列出如下： ① Newton-Cotes 公式
当1n =时为[()()]2b a
T f a f b -=
+称为梯形公式； 当2n =时为[()4()()]62
b a b a
S f a f f b --=++称为Simpson 公式；
腹有诗书气自华 当4n =时为01234[7()32()12()32()7()]90
b a
C f x f x f x f x f x -=
++++称为Cotes 公式. ② 复化求积公式
复化梯形公式：1
1
101
[()()][()2()()]22n n n k k k k k h h
T f x f x f a f x f b --+===+=++∑∑
复化Simpson 公式：
1
11
11100122
[()4()()][()4()2()()]66n n n n k k k k k k k k h h
S f x f x f x f a f x f x f b ---+++====++=+++∑∑∑
复化Cotes 公式：
1111
1130001424
[7()32()12()32()14()7()]90n n n n n k k k k k k k k h
C f a f x f x f x f x f b ----+++=====+++++∑∑∑∑
③ Gauss-Legender 公式
令对()1f x =准确成立，定02A =则： 一点Gauss-Legender 公式为：1
1()2(0)f x dx f -≈⎰
两点Gauss-Legender
公式为：1
1
()(f x dx f f -≈+⎰ 三点Gauss-Legender
公式为：1
1
585()()(0)(95995
f x dx f f f -≈
-++⎰ 2）对常微分方程求解的准备
在生产实践和科学研究过程中，我们经常会遇到求解常微分方程或常微分方程组的定解问题.在许多情况下，首先找到问题的解析解，然后再进行相关的计算往往非常困难，有时甚至是行不通的.我们已经可以成功的绕过牛顿—莱布尼茨公式而直接求一般的定积分的数值解.基于同样的理由，对于常微分方程或常微分方程组来说，我们也可以避免求解析解而直接求相应的数值解.
一般的高阶常微分方程都可以通过相应的变量代换而化为一阶常微分方程组来求解.一阶常微分方程组的初值问题求数值解与一阶常微分方程的初值问题求数值解的方法基本相同.
求解一阶常微分方程的初值问题求数值解的公式有Euler 格式、后退的Euler 格式、梯形格式、改进的Euler 格式、Euler 两步格式、Taylor 级数法和Runge-Kutta 方法等.
腹有诗书气自华
下面将几种常用的公式简单列出如下： ① Euler 格式：1(,)n n n n y y hf x y +=+
② 后退的Euler 格式：111(,)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 +++=++
改进的Euler 格式：
预测 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 +++=++
④ Euler 两步格式：112(,)n n n n y y hf x y +-=+ 改进的Euler 两步格式： 预测 112(,)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 +++=++
⑤ Taylor 级数法：2'''
()12!
!
p p n n n n n h h y y hy y y p +=+++
+ ⑥ Runge-Kutta 方法：
二阶Runge-Kutta 方法：1112212
1()
(,)(,)n n n n n p n y y h K K K f x y K f x y phK λλ++⎧=++⎪=⎨⎪=+⎩ 其中满足：122112p λλλ+=⎧⎪
⎨=⎪⎩
三阶Runge-Kutta 方法：1112233121312()
(,)(,)(,())
n n n n n n n n y y h K K K K f x y K f x ph y phK K f x qh y qh rK sK λλλ+=+++⎧⎪=⎪⎨=++⎪⎪=+++⎩
腹有诗书气自华
其中满足：1232322
2
33111213
16r s p q p q pqs λλλλλλλλ⎧⎪
⎪+=⎪
++=⎪⎪⎪
+=⎨⎪
⎪+=⎪⎪
⎪=⎪⎩
四阶Runge-Kutta 方法：112341213
243(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 +⎧
=++++⎪⎪
=⎪⎪⎪=++⎨⎪⎪
=++⎪⎪=++⎪⎩
5.2 建立模型
首先，对常微分方程中的积分20
()t
y s ds ⎰
进行求解.
我们选用梯形公式[()()]2b a T f a f b -=
+来求解，则有200()[(0)()]2
t t y s ds f f t -=+⎰.
将其代入'
20
()()t
y t ay b
y s ds =+⎰
得：
'0(0)()(|)()222
t b t bt bt
y t ay y y a y =-=+
+=++ 然后，求解常微分方程时，用四阶Runge-Kutta 方法：
1
12341213
243(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 +⎧
=++++⎪⎪
=⎪⎪⎪
=++⎨⎪
⎪
=++⎪⎪=++⎪⎩ 则有：1
1234121
3
243(22)6
()22()()222()()222
()()22n n h y y K K K K bt bt K a y bt h bt K a y K bt h bt K a y K bt bt
K a y hK +⎧=++++⎪⎪
⎪=++⎪⎪⎪=+++⎨⎪
⎪
=+++⎪⎪
⎪=+++⎪⎩
模型建立的具体过程如下流程图：
腹有诗书气自华
6 模型求解
求解数值积分与常微分方程的方法有多种，同样应用数学软件进行求解时，也有多种软件可用，如常用的有MATLAB 软件、Lingo 软件等.这两种软件对目前我们面对的数学问题，几乎都可以用其来解决.Lingo 软件在求解线性规划问题时用的较多，一般的问题都习惯用MATLAB 软件进行求解.
6.1 MATLAB 软件求解结果
由于在用四阶Runge-Kutta 方法求解时，需要设初始值，所以设01y =，
1 1.123y =步长0.2h =.又由于
0,0a b ><，为了便于用软件进行求解，不妨设
2,2a b ==-
计算了1,2,
,10t =时，方程的值.所得结果如下表所示：
表1
7 模型验证及结果分析
7.1 模型验证
7.1.1 MATLAB软件求解结果验证
由于这只是一个方程的求解，没有实际背景，故没有真值可以用来进行比较.只能用不同的方法进行求解，来看不同方法求解结果相差多少来验证.
在实际生产实践中，我们总是对某一个实际问题数学化，然后建立数学模型来解决.并且在实际生产中，我们可以得到足够多的数据.再对模型求解后，就可以通过比较实际值与模型求解值来验证模型的合理性.但现在所做的题目没有实际背景，没有真值来作比较.故只能通过用不同的方法求解，比较所得结果，验证那种方法更好.
7.2 结果分析
在模型求解时，我们假设
01
y=，
11.123
y=.假设的值不是真值，故对模型的求解有一定的影响.并且只取了一个步长0.2
h=，不一定是最优的步长.
8 模型的优缺点
8.1 模型的优点
模型的优点是列出了多种求解数值积分和常微分方程的方法，有多种选择来进行求解.
8.2 模型的缺点
在模型求解时，我们假设
01
y=，
11.123
y=.假设的值不是真值，故对模型的求解有一定的影响.并且只取了一个步长0.2
h=，不一定是最优的步长.
9 模型的推广
腹有诗书气自华
在实际生产实践中，我们总是对某一个实际问题数学化，然后建立数学模型来解决.并且在实际生产中，我们可以得到足够多的数据.再对模型求解后，就可以通过比较实际值与模型求解值来验证模型的合理性.本问题中的常微分方程，可能适合与生活实践中的某个数学问题，那么可以用本文中的方法进行求解.
参考文献
[1] 李庆扬，王能超，易大义编著．数值分析（第四版）[M]．武汉：华中科技大学出
版社，2006.7.
[2]刘卫国编著.MATLAB程序设计教程（第二版）[M]. 北京：中国水利水电出版社，
2005.3.
[3] 姜启源，谢金星，叶俊编著．数学模型（第三版）[M]．北京：高等教育出版社，
2005：1-202.
[4] 张可村，赵英良编著. 数值计算的算法与分析[M]. 北京：科学出版社，2004.6.
[5] 甄西丰编著.实用数值计算方法 [M]. 北京：清华大学出版社，2007.2.
[6] 关治，陆金甫编著. 数值方法[M]. 北京：清华大学出版社，2006.12.
腹有诗书气自华
附录
1）程序1：
clear;
t=[];y=[];a=2;b=-2;
t(1)=0;y(1)=1;y(2)=1.123;
for j=1:10
for i=2:6
t(i)=t(i-1)+0.2;
i=i+1;
end
for i=2:6
K1=(a+b*t(i-1)/2)*y(i-1)+b*t(i-1)/2;
K2=(a+b*t(i-1)/2)*(y(i-1)+0.1*K1)+b*t(i-1)/2;
腹有诗书气自华
K3=(a+b*t(i-1)/2)*(y(i-1)+0.1*K2)+b*t(i-1)/2; K4=(a+b*t(i-1)/2)*(y(i-1)+0.2*K3)+b*t(i-1)/2; y(i)=y(i-1)+0.2/6*(K1+2*K2+2*K3+K4)+b*t(i-1)/2; i=i+1;
end
y(6)
y(1)=y(5);y(2)=y(6);t(1)=t(6);j=j+1;
end
2）求得的结果如下：
ans =
1.6266
ans =
-6.0344
ans =
-14.0483
ans =
-15.0195
ans =
腹有诗书气自华
-13.3441
ans =
-12.5391
ans =
-12.5193
ans =
-12.9037
ans =
-13.5898
ans =
-14.6088出师表
两汉：诸葛亮
腹有诗书气自华
先帝创业未半而中道崩殂，今天下三分，益州疲弊，此诚危急存亡之秋也。

然侍卫之臣不懈于内，忠志之士忘身于外者，盖追先帝之殊遇，欲报之于陛下也。

诚宜开张圣听，以光先帝遗德，恢弘志士之气，不宜妄自菲薄，引喻失义，以塞忠谏之路也。

宫中府中，俱为一体；陟罚臧否，不宜异同。

若有作奸犯科及为忠善者，宜付有司论其刑赏，以昭陛下平明之理；不宜偏私，使内外异法也。

侍中、侍郎郭攸之、费祎、董允等，此皆良实，志虑忠纯，是以先帝简拔以遗陛下：愚以为宫中之事，事无大小，悉以咨之，然后施行，必能裨补阙漏，有所广益。

将军向宠，性行淑均，晓畅军事，试用于昔日，先帝称之曰“能”，是以众议举宠为督：愚以为营中之事，悉以咨之，必能使行阵和睦，优劣得所。

亲贤臣，远小人，此先汉所以兴隆也；亲小人，远贤臣，此后汉所以倾颓也。

先帝在时，每与臣论此事，未尝不叹息痛恨于桓、灵也。

侍中、尚书、长史、参军，此悉贞良死节之臣，愿陛下亲之、信之，则汉室之隆，可计日而待也。

臣本布衣，躬耕于南阳，苟全性命于乱世，不求闻达于诸侯。

先帝不以臣卑鄙，猥自枉屈，三顾臣于草庐之中，咨臣以当世之事，由是感激，遂许先帝以驱驰。

后值倾覆，受任于败军之际，奉命于危难之间，尔来二十有一年矣。

先帝知臣谨慎，故临崩寄臣以大事也。

受命以来，夙夜忧叹，恐托付不效，以伤先帝之明；故五月渡泸，深入不毛。

今南方已定，兵甲已足，当奖率三军，北定中原，庶竭驽钝，攘除奸凶，兴复汉室，还于旧都。

此臣所以报先帝而忠陛下之职分也。

至于斟酌损益，进尽忠言，则攸之、祎、允之任也。

愿陛下托臣以讨贼兴复之效，不效，则治臣之罪，以告先帝之灵。

若无兴德之言，则责攸之、祎、允等之慢，以彰其咎；陛下亦宜自谋，以咨诹善道，察纳雅言，深追先帝遗诏。

臣不胜受恩感激。

今当远离，临表涕零，不知所言。

腹有诗书气自华。