081数值计算方法—常微分方程(组)

科学计算—理论、方法

及其基于MATLAB 的程序实现与分析

微分方程（组）数值解法

§1 常微分方程初值问题的数值解法

微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton 第二定律：

()()()?????'='==0

00022x t x x t x t F dt x

d m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：

()

()0

0y y t f ay dt

dy

=+= （２） 的解为：

()()()τττd f e y e t y t

t a at

?-+=00 （３）

但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：

１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；

２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；

如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。

一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：

()()0

00,y t y t t t y t F dt

dy

f

=≤≤= （７）

其中

()()()()????

??

? ??=t y t y t y t y n 21 （８）

()()()()????

??

? ??=y t f y t f y t f y t F n ,,,,21

（９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。

§１.1 常微分方程（组）的Cauch 问题数值解法概论

假设要求在点（时刻）k t ，n k ,,2,1 =处初值问题（７）的解

()00,,y t t y y =的（近似）值,如果已求得k t 时刻的值()k t y 或它的近似值k y (如0t 时刻的值()0t y )，那么将式（７）的两端在区间[]1,+k k t t 上积分

()()()()dt t y t f t y t y dt dt

dy

k k

k k

t t k k t t ??

++=-=+11

,1 (10)

可得

()()()()dt t y t f t y t y k k

t t k k ?

++

=+1

,1 (11)

或

()()()dt t y t f y y t y k k

t t k k k ?

++

=≈++1

,11 (12)

显然,为了利用式（11）或（12）求得()1+k t y 的精确值（近似值1+k y ），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解, 可以说,不同的式值积分方法将给出不同的Cauch 问题的数值解法。 §1.2 最简单的数值解法——Euler 方法

假设要求在点（时刻）kh t t k +=0，n

t t h f 0

-=

，n k ,,2,1 =处初值问

题（７）的解()t y y =的近似值。首先对式（７）的两端积分，得

()()()()dt t y t f t y t y dt dt

dy

k k k k

t t k k t t ??

++=-=+11

,1 （13） 对于式（13）的右边，如果用积分下限k t t =处的函数值()()k k t y t f ,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有

()()()()()[]k k t t

k k t y t hf dt t y t f t y t y k k

,,1

1≈=-?++ （14）

进而得到下式给出的递推算法—Euler 方法

()()0

01,,2,1,y t y n

k y t hf y y k k k k ==+=+ （15）

例１ 用Euler 方法解如下初值问题，取3.0=h ，

()1

03

0sin 23=≤≤-=y t t y y dt

dy

解：由（15）得

()1

010

,,2,1sin 6.03.13

1==-=+y k t y y y k

k k k

结果如下：

如果取1.0=h ，其结果如下图所示：

Euler_Method

§1.3 改进的Euler 方法

对于（15）的右边，如果被积函数用积分限k t t =和1+=k t t 处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有

()()()()()()()()[]111,,2

,1++++≈

=-?

+k k k k t t k k t y t f t y t f h

dt t y t f t y t y k k

（16） 进而得到下式给出的递推算法：

()()[]()0

0111,,2,1,,2

y t y n

k y t f y t f h

y y k k k k k k ==++=+++ （17）

通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求1+k y 时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了

１） 尽可能地保持算法（17）精度高的优点； ２） 尽可能地利用Euler 方法计算简单的长处； 人们采取了如下的称之为改进的Euler 方法的折衷方案：

预测 ()k k k k y t hf y y ,01+=+ （18） 修正 ()()[]

n k y t f y t f h y y k k k k k k ,,2,1,,2

111 =++

=+++ （19）

()00y t y =

例２ Euler 方法与改进的Euler 方法的比较 下图是当3.0=h 时比较的结果：

open Improved_Euler_Method.m

§1.4 Euler 方法和改进的Euler 方法的误差分析

由Taylor 公式

()()()()()()()()()()

2

2

111,!

21h O h t y t f t y t t t y t t t y t y t y k k k k k k k k k k k ++=+-''+-'+=+++

（19） 说明Euler 方法的截断误差是()2h O ，类似地，由

()()()()()()

321!

21

,h O h t y h t y t f t y t y k k k k k +''+

+=+ （20） ()()()()()()

321111!

21

,h O h t y h t y t f t y t y k k k k k +''+

-=++++ （21） 以及

()()()()()()h O t y t y h t y t y t y k k k k k +''=''?+'''+''=''++11 (22)

让式（20）的两端减式（21）的两端，可得

()()()()()()[]()()

()()()()()[]()

3

1132111,,2

4

1,,2h O t y t f t y t f h

t y h O h O h t y t f t y t f h

t y t y k k k k k k k k k k k +++=++++

=+++++ （23）

从上述推导Euler 方法、改进的Euler 方法的过程以及例１、例２容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于： 1 在推导Euler 方法时，我们是用待求解函数()t y y =在一点处的变化 率()()k k t y t f ,代替()t y y =在区间],[1+k k t t 上的平均变化率：

()()()()t y t hf dt t y t f k k

t t ,,1

=?

+ （24）

2 而在推导改进的Euler 方法时，我们是用待求解函数()t y y =在两点处变化率的平均值()()()()[]11,,2

1+++k k k k t y t f t y t f 代替()t y y =在区间],[1+k k t t 上的平均变化率；

显然，通常()()()()[]11,,2

1+++k k k k t y t f t y t f 比()()k k t y t f ,更接近于()t y y =在区间],[1+k k t t 上的平均变化率。由此启发人们：适当地选取区间],[1+k k t t 上函数()t y y =若干点处的变化率，用它们加权平均值代替()t y y =在区间

],[1+k k t t 上的平均变化率，近似解的精度应更高。

下面将要介绍的Runge —Kutta 法就是基于上述想法得到的。

§2 Runge —Kutta 法

Runge —Kutta 法是按选取区间],[1+k k t t 上函数()t y y =变化率()y t f ,的个数的多少和截断误差的阶数()m h O 来区分的一系列方法，如 １ 二阶的Runge —Kutta 法（改进的Euler 方法）

()()()???

?

???

++=+==++211

11212,,K K h

y y hK y t f K y t f K k k k k k k （25） ２ 三阶的Runge —Kutta 法

()

()()????

????

?+++=++=??? ??

++==+32112

312146,2,2,K K K h y y hK y h t f K K h y h t f K y t f K k k k k k k k k （26） ３ 四阶的Runge —Kutta 法 1) 古典形式

()()()????

?

??

????++++=++=?

?? ??++=?

?? ??

++==+432113423121226,2,22,2,K K K K h

y y hK y h t f K K h y h t f K K h y h t f K y t f K k

k k k k k k k k k （27） 2) Gill 公式(具有减小舍入误差的优点)

()

()()()

??????

???

????+++-++=??

?

? ??++-+=???

? ??-+-++=?

?? ??

++==+4

32113242131212222622222,22

2212,22,2,K K K K h

y y hK hK y h t f K hK hK y h t f K K h y h t f K y t f K k k k k k k k k k k （28） 4 Runge —Kutta 法的一般形式

()∑=++=+=n

i i i k k k k k K a h y h y t h y y 1

1,,φ (29)

其中()h y t k k ,,φ称为增量函数(Increment Function)以及

()

()

()()

?????????++++=

+++=++==

-----11,111,1122212123111121,,,,n n n n k n k n k k k k k k hK q hK q y h p t f K hK q hK q y h p t f K hK q y h p t f K y t f K （30）

需要特别指出的是：在确定（29）、（30）中的参数i a ，i p 和ij q 时，应该使（29）的右端和适当阶的（如n 阶）Taylor 展式一致，这样，

至少对于底阶的R-K 法来说，参与加权的斜率个数n 与方法的阶数是一致的。

例如：一阶的R-K 法即Euler 法，局部的截断误差（Truncation Error ）是()2h O ，所以，当微分方程的阶是一次函数时是精确的；二阶R-K 法即改进的Euler 法，局部的截断误差是()3h O ，所以，当微分方程的阶是二次函数时是精确的；类似地，四阶R-K 法局部的截断误差是()5h O 。

但是，一般五阶以及五阶以上的R-K 法的局部截断误差不再具有上述的特点，因此用“性价比”来衡量，四阶R-K 得到了广泛的应用。

下面是Butcher's Fifth-Order R-K 方法（1964）：

()6543117321232790

K K K K K h

y y k k +++++

=+ （28） ()

??

???

??

?

??

???

???

?

?

??+-++-+=?

?

? ??+-+=?

?? ??+-+=?

?? ??+++=?

?? ?

?

++==543216415324213121787127127273,169163,4321,28181,44,4,hK hK hK hK hK y h t f K hK hK y h t f K hK hK y h t f K hK hK y h t f K K h y h t f K y t f K k k k k k k k k k k k k （29）

例三：疾病的传播与防疫模型 open Lunge_Kutta_Method.m

其中()t x 表示一种传染病传染者的人数，()t y 表示易感染者人数。

?????++--=-=e dx y cx axy dt

dy

bx

axy dt dx

2 (31) 该系统的Jacobi 矩阵为：

()??

?

?

??--+---=??=

22,cx ax d

cxy ay ax b ay y x f

J (32) 该系统的平衡点为：

()()()()?

???

????????? ??+---=???? ???

???

????????? ??+-+-=???? ??a b acbe d b a b d a bc y x a b acbe d b a b d a bc y x 42142122112200 (33)

根据问题的实际意义，有意义的平衡点为：

()()()???

????? ?

?+-+-=a b acbe d b a b d a bc y x ,421,2

200

(34)

§3 解刚性问题的数值方法

例四：单个微分方程的例子

考虑如下的一阶微分方程描述的系统

()?????=+--=-0

03000

20001000x x e x dt

dx

t (35) 容易求的该微分方程初值问题的通解：

()39992000999997,010000+-

??? ?

?-=--t t e x e x t x (36) f=dsolve('Dy = -1000*y+3000-2000*exp(-t)', 'y(0) = c')

并且有:

()3,lim 0=+∞

→x t x t (37)

当初值999

997

0≠

x 时，方程的解中包含“快变”和“慢变”两部分. 应用不同数值方法求解的结果：

１ 解非刚性问题的方法：ODE45: Sol_StiffODE01.m ２ 解刚性问题的方法： ODE15s: Sol_StiffODE01.m ３ 解非刚性问题的方法：Euler Method: Sol_StiffODE01_Euler.m

系统（35）所对应的自治系统为

()?????=-=0

01000x x x

dt

dx

(38) 利用Euler 法求解上述初值问题：

()()()k y h k y 100011-=+， ,2,1,0=k (39)

为保证时间序列()k y 的收敛性，步长应满足不等式：001.0

考虑如下的一阶微分方程组描述的系统

()()??

?==???? ????????--=????????? ???-=+-=0

0003011003530110035y y x x y x y x y x dt

dy y x dt dx

(40) ()

()

??

????+-++=--------t t t t

t

t t t At

e e e e

e e e e e

0101.3029899.30101.3029899.30101.3029899.30101.3029899.39966.00034.03355.00101.09966.00034.0 (41)

()()()()?

??

? ??+-+++++=?

??

? ??=???? ??----t t t

t

At e y x e y x e

y x e y x y x e y x 0101.302009899.3000101.302009899.300009966.03355.00034.03355.00101.09966.00101.00034.0 (42)

利用Euler 法求解上述初值问题的迭代公式为

()()k y h h h h

k y ??

????---+=+3011100131511， ,2,1,0=k (43) 由于系数矩阵的特征值为

()h h h 0101

.1491531555121531±+≈±+=λ (44) 所以Euler 法不适用于该初值问题的求解。

利用四阶R-K 法求解上述初值问题的迭代公式为

()()k My k y =+1， ,2,1,0=k (45)

其中系数矩阵为

4

4332232

3252A h A h A h hA E M ++++= (46)

其特征值为：

()4

4332232

32521A A A A M h h h h h λλλλλ++++= (47)

１ 解非刚性问题的方法：ODE45: Sol_StiffODE02.m ２ 解刚性问题的方法： ODE15s: Sol_StiffODE02.m ３ 解非刚性问题的方法：Euler Method: Sol_StiffODE02_Euler.m

４ 解非刚性问题的方法：Euler Method: Sol_StiffODE02_RK4.m

f=dsolve('Dx=-5*x+3*y','Dy=100*x-301*y','x(0)=1','y(0)=2') transient dominate motivation orientation scope

discretization

predictor – corrector approach increment function recurrence

adaptive steps - size control abrupt change

例2：二体问题

open Runge_Kutta_Method01

§4 线性多步法与常微分方程数值解法的分类

一般的常微分方程初值问题的数值解法都是以递推的形式给出的，即是递推算法，根据递推算法，在计算1+k y 时，已经得到了前面各时刻的近似值i y ，k i ,,1,0 =，前面介绍的各种数值解法都有一个共同点：在计算1+k y 时，只用到了前一个时刻（当前时刻k t ）的“信息”k

y

预测未来某时刻1+k t 系统的“状态”()1+k t y ，这样的数值解法称为单步法，对于一个动态过程()t y y =在1+k t 时刻的状态()1+k t y 而言，不仅前一个时刻的信息k y 对它有影响，前若干个时刻的信息通常对它也有影响，显然，单步法的缺点是没有充分利用已得到的信息。

基于上述考虑，适当取前若干个时刻的信息i y ，m k k k i --=,,1, ，并用()i i y t f ,的线性组合代替()t y y =在区间],[1+k k t t 上的平均变化率所给出的方法称为线性多步法。

§4.1 线性多步法: 从二步法谈起

在改进的Euler 方法中

预测 ()k k k k y t hf y y ,01+=+ （48） 较正 ()()[]

n k y t f y t f h y y k k k k k k ,,2,1,,2

111 =++

=+++ （49）

()00y t y =

如果改进预测的手段，如利用数值微分的研究成果：

()()()

()

()()()()()

3112

112,2h O h t y t f t y t y h O h

x f x f x f k k k k k k k ++=?+-=

'-+-+ (50) 那么，由于（48）的局部截断误差是()3h O ，精度高于Euler 公式的局部截断误差()2h O ，所以，下列预测-校正方法

Predictor: ()k k k k y t hf y y ,2101+=-+ （51） Corrector: ()()[]

n k y t f y t f h y y k k k k k k ,,2,1,,2

111 =++

=+++ （52）

()00y t y =

尽管步长增大了，精度却提高了。公式（49）、（50）给出的数值方法称为非自始的（non-self-starting ）二步法，其原因在于：利用该方法求解微分方程初值问题时，自身不能提供所需的初值1-y 。

§4.3 线性二步法的进一步改进

在上述讨论中，我们只考虑了预测手段方面的改进，从方法论的角度进行思考，也应该考虑校正方面的改进问题，并且，可以希望适当的改进会带来精度提高的回报。下面推荐的方案就是常用的一种： Predictor: ()m k k m k k y t hf y y ,2101+=-+ （53） Corrector:

()()[]

m

i y t f y t f h y y i k k m k k m k m k ,,2,1,,2

1

1

11 =++

=-+++ , （54） 其中m k y 1+是1+k t 时刻最终计算结果,校正的次数m 按如下公式确定

εε≤?-=+-++%1001

1

11i k i k i k a y y

y （55）

其中

ε是预先指定的误差精度（Prespecified Error Tolarence ）.

§4.4 线性二步法的误差分析与步长控制

1 预测误差：

()()()()()()p k p p k k k k k t f h E t f h h t y t f t y t y ξξ,3

1,3

12,3

3

11'''=?'''++=-+ (56) 其中下表表示预测误差，这表明

()()p k k k t f h y

t y ξ,3

13

1

1'''+=++ （57） 2 校正误差：类似地可得

()()()

c k c c k m k k t f h E t f h y

t y ξξ,12

1,12

13

3

1

1'''-=?'''-=++ （58）

3 步长控制准则

有误差估计（53）-（55）得到

()ξ,12

1531

01k m k k t f h y y E '''-=-=++ （59）

如果E 大于事先指定的误差精度，那么就减小步长。

§4.5 线性多步法的进一步扩展

在积分公式

()()()()()()dt

t y t f y y dt

t y t f t y t y dy y k n

k k n

k k n

k t t n k k t t n k k y y ?

?

?

+-+-+-+=?=-='-+-+11

1

,,11 （60）

中，用1+n 个节点()()i i i i y t f t ,,，1,,1,+---=k n k n k i ，的n 次插值多项式()t P n 代替被积函数()()t y t f ,，得到的就是一般的线性多步法。

如当1=n 时，

()k k k k y t hf y y ,211+=-+ （61）

如当2=n 时，

()()[]k k k k k k y t f y t f h

y y ,,2

31121++

=---+ （62） 特别地，用四个节点()()i i i i y t f t ,,，k k k k i ,1,2,3---=，的三次插值多项式()t P 3代替被积函数()()t y t f ,，得到

()3211937595524

---+-+-+

=k k k k k k f f f f h

y y （63） 四个节点()()i i i i y t f t ,,，1,,1,2+--=k k k k i ，的三次插值多项式()t P 3代替，又得到

()2111519924

--+++-++

=k k k k k k f f f f h

y y （64） 其中误差分别为

()()ξ55720251y h 、()

()ξ55720

19y h -，可见，公式（64）的精度比公式（63）的精度高。 Adams 预报-校正公式：

()()

??

???

+-++=-+-+=--++---+21113211

519924937595524n n n n n n n n n n n n f f f h

y y f f f f h y y （65） 其中，()n n n y t f ,，()111,+++n n n y t f ，

说明：公式（65）自身不能提供最初所需的三个初值，需要Runge-Kutta 法的配合。

§5 应用常微分方程组初值问题数值解法的基本要求 §5.1 常微分方程组初值问题数值解法的本质

以Runge-Kutta 为例，对于给定的常微分方程（组）的初值问题：

()()0

00,y t y t t t y t F dt

dy

f

=≤≤= （７）

其一般形式是

()∑=++=+=n

i i i k k k k k K a h y h y t h y y 1

1,,φ (29)

从一般性的系统角度看，（7）刻画的是连续时间动态（力）系统的演化规律，而（29）刻画的是离散时间动态系统的演化规律，由此可见，常微分方程（组）初值问题数值解法的本质是用离散时间动态系统模型代替连续时间动态系统模型，这意味着，应用数值方法求解对于给定的常微分方程（组）的初值问题时，所选择的离散时间动态系统模型应该是连续时间动态系统模型的“有意义的近似”。 为此，应用常微分方程组初值问题数值解法需要了解的三个基本概念： 1 相容性：

因为

()h y t h y y k k k k ,,1φ=-+ (66)

所以，应该保证

()()()()()k k k k k k h t y t f t y t h

t y y ,0,,lim

10

==-+→φ (67)

2 收敛性：

如果已知k t 时刻函数的精确值()k t y ，

()()()h t y t h t y y k k k k ,,1φ+=+ (68)

那么应该有

()k k h t y y =+→10

lim (69)

3 稳定性：

§5 常微分方程边界值问题的数值解法

对于给定的n 阶常微分方程组

()f t t t y t F dt

dy

≤≤=0, (70)

其定解问题要求给定n 个独立的定解条件,如果n 个独立的定解条件指定在求解区间的一端,那么这样的定解问题就是前面讨论过的初值问题,如果n 个独立的定解条件分别在求解区间的两端给定,这样的定解问题就是所谓的边界值问题(Two-Points Boundary Valve Problems ).常微分方程的边界值问题在工程实际中广泛存在，如两端固定的弹簧振动问题，与外界有热交换的杆状物体的热平衡问题等。

下面介绍两种常用的常微分方程边界值问题的数值解法。

§5.1 边界值问题的数值解法一: Shooting Method

采用 Shooting Method 解常微分方程边界值问题的主导思想是将边界值问题转化为等价的初值问题。 例6 考虑如下的边界值问题：

()L x T T h dx T

d a ≤≤-'=02

2 （71）

其中L 是细长杆状物体的长度，a T 时周围空气的温度，h '时热交换系数（单位：2-cm ），那么（71）就是描述该物体处于稳态时温度变化规律的数学模型，如果在杆的两端，外界的温度固定：

()()?

?

?==L T L T T T 0

0 （72） 那么式（71）和（72）就构成了简单的常微分方程边界值问题，设参数和边界条件如下：

081数值计算方法—常微分方程(组)

科学计算—理论、方法 及其基于MATLAB 的程序实现与分析 微分方程（组）数值解法 §1 常微分方程初值问题的数值解法 微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton 第二定律： ()()()?????'='==0 00022x t x x t x t F dt x d m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题： () ()0 0y y t f ay dt dy =+= （２） 的解为： ()()()τττd f e y e t y t t a at ?-+=00 （３） 但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：

１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解； ２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）； 如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。 一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题： ()()0 00,y t y t t t y t F dt dy f =≤≤= （７） 其中 ()()()()???? ?? ? ??=t y t y t y t y n 21 （８） ()()()()???? ?? ? ??=y t f y t f y t f y t F n ,,,,21 （９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。 §１.1 常微分方程（组）的Cauch 问题数值解法概论

高阶线性微分方程常用解法介绍

高阶线性微分方程常用解法简介 关键词：高阶线性微分方程 求解方法 在微分方程的理论中，线性微分方程是非常值得重视的一部分内容，这不仅 因为线性微分方程的一般理论已被研究的十分清楚，而且线性微分方程是研究非线性微分方程的基础，它在物理、力学和工程技术、自然科学中也有着广泛应用。下面对高阶线性微分方程解法做一些简单介绍. 讨论如下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 λαβ=-也是特征根，因而于这对共轭复根

常微分方程边值问题的数值解法

第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 αβ''-=<

常微分方程考研讲义 一阶微分方程解的存在定理

第三章一阶微分方程解的存在定理 [教学目标] 1.理解解的存在唯一性定理的条件、结论及证明思路，掌握逐次逼近法，熟练 近似解的误差估计式。 2.了解解的延拓定理及延拓条件。 3.理解解对初值的连续性、可微性定理的条件和结论。 [教学重难点] 解的存在唯一性定理的证明，解对初值的连续性、可微性定理的 证明。 [教学方法] 讲授，实践。 [教学时间] 12学时 [教学内容] 解的存在唯一性定理的条件、结论及证明思路，解的延拓概念及延 拓条件，解对初值的连续性、可微性定理及其证明。 [考核目标] 1.理解解的存在唯一性定理的条件、结论，能用逐次逼近法解简单的问题。 2.熟练近似解的误差估计式，解对初值的连续性及可微性公式。 3.利用解的存在唯一性定理、解的延拓定理及延拓条件能证明有关方程的某些性质。 §1 解的存在性唯一性定理和逐步逼近法 微分方程来源于生产实践际，研究微分方程的目的就在于掌握它所反映的客 观规律，能动解释所出现的各种现象并预测未来的可能情况。在第二章介绍了一 阶微分方程初等解法的几种类型，但是，大量的一阶方程一般是不能用初等解法 求出其通解。而实际问题中所需要的往往是要求满足某种初始条件的解。因此初 值问题的研究就显得十分重要，从前面我们也了解到初值问题的解不一定是唯一的。他必须满足一定的条件才能保证初值问题解的存在性与唯一性，而讨论初值 问题解的存在性与唯一性在常微分方程占有很重要的地位，是近代常微分方程定 性理论，稳定性理论以及其他理论的基础。 例如方程 过点(0,0)的解就是不唯一，易知0 y=是方程过(0,0)的解，此外，容易验证，2 =或更一般地，函数 y x 都是方程过点(0,0)而且定义在区间01 <<的任一数。 c ≤≤上的解，其中c是满足01 x

常微分方程组(边值)

常微分方程组边值问题解法 打靶法Shooting Method （shooting.m ） %打靶法求常微分方程的边值问题 function [x,a,b,n]=shooting(fun,x0,xn,eps) if nargin<3 eps=1e-3; end x1=x0+rand; [a,b]=ode45(fun,[0,10],[0,x0]'); c0=b(length(b),1); [a,b]=ode45(fun,[0,10],[0,x1]'); c1=b(length(b),1); x2=x1-(c1-xn)*(x1-x0)/(c1-c0); n=1; while (norm(c1-xn)>=eps & norm(x2-x1)>=eps) x0=x1;x1=x2; [a,b]=ode45(fun,[0,10],[0,x0]'); c0=b(length(b),1); [a,b]=ode45(fun,[0,10],[0,x1]'); c1=b(length(b),1) x2=x1-(c1-xn)*(x1-x0)/(c1-c0); n=n+1; end x=x2; 应用打靶法求解下列边值问题： ()()??? ????==- =010004822y y y dx y d 解：将其转化为常微分方程组的初值问题

()????? ? ?????==-==t dx dy y y y dx dy y dx dy x 0011221 048 命令： x0=[0:0.1:10]; y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解 plot(x0,y0,'r') hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1))

郑州大学研究生课程数值分析复习---第八章 常微分方程数值解法

郑州大学研究生课程（2012-2013学年第一学期）数值分析 Numerical Analysis 习题课 第八章常微分方程数值解法

待求解的问题：一阶常微分方程的初值问题/* Initial-Value Problem */: ?????=∈=0 )(] ,[),(y a y b a x y x f dx dy 解的存在唯一性（“常微分方程”理论）：只要f (x , y ) 在[a , b ] ×R 1 上连续，且关于y 满足Lipschitz 条件，即存在与x , y 无关的常数L 使 对任意定义在[a , b ] 上的y 1(x ) 和y 2(x ) 都成立，则上述IVP 存在唯一解。 1212|(,)(,)||| f x y f x y L y y ?≤?一、要点回顾

§8.2 欧拉(Euler)法 通常取(常数),则Euler 法的计算格式 h h x x i i i ==?+1?? ?=+=+) (),(001x y y y x hf y y i i i i i =0,1,…,n ( 8.2 )

§8.2 欧拉(Euler)法(1) 用差商近似导数 )) (,()()()()(1n n n n n n x y x hf x y x y h x y x y +=′+≈+?? ?=+=+) (),(01a y y y x hf y y n n n n 差分方程初值问题向前Euler 方法h x y x y x y n n n ) ()()(1?≈ ′+)) (,() ()(1n n n n x y x f h x y x y ≈?+))(,()(n n n x y x f x y =′

第十一章 常微分方程边值问题的数值解法汇总

第十一章 常微分方程边值问题的数值解法 工程技术与科学实验中提出的大量问题是常微分方程边值问题．本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法. 11.1 引言 在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程 ),,(y y x f y '='', b x a <<, (11.1.1) 在如下三种边界条件下的定解问题: 第一种边界条件: α=)(a y , β=)(b y (11.1.2) 第二种边界条件: α=')(a y , β=')(b y (11.1.2) 第三种边界条件: ? ? ?=-'=-'101 0)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a . 常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法. 11.2 打靶法 对于二阶非线性边值问题 ()()().,,βα==≤≤'=''b y a y b x a y y x f y ，，， (11.2.1) 打靶法近似于使用初值求解的情况． 我们需要利用一个如下形式问题初值解的序列： ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α， （11.2.2） 引进参数v 以近似原边界值问题的解．选择参数k v v =，以使： ()()β==∞ →b y v b w k k ,lim ， （11.2.3）

其中),(k v x w 定义为初值问题（11.2.2）在k v v =时的解，同时()x y 定义为边值问题（11.2.1）的解． 首先定义参数0v ，沿着如下初值问题解的曲线，可以求出点),(αa 对应的初始正视图 ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α． （11.2.4） 如果),(0v b w 不严格收敛于β，那么我们选择1v 等值以修正近似值，直到),(0v b w 严格逼近β． 为了取得合适的参数k v ，现在假定边值问题（11.2.1）有唯一解，如果),(v x w 定义为初始问题（11.2.2）的解，那么v 可由下式确定： 0),(=-βv b w ． （11.2.5） 由于这是一个非线性方程，我们可以利用Newton 法求解．首先选择初始值0v ，然后由下式生成序列 ),)(()),((111----- =k k k k v b dv dw v b w v v β，此处),(),)(( 11--=k k v b dv dw v b dv dw ， （11.2.6） 同时要求求得),)(( 1-k v b dv dw ，因为),(v b w 的表达式未知，所以求解这个有一点难度；我们只能得到这么一系列的值。 ，，，),(),(),(),(1210-??k v b w v b w v b w v b w 假如我们如下改写初值问题（11.2.2），使其强调解对x 和v 的依赖性 ()()v v a w v a w b x a v x w v x w x f w ='=≤≤'=''),(,),(),,(,，，，α，（11.2.7） 保留初始记号以显式与x 的微分相关．既然要求当k v v =时),)((v b dv dw 的值，那么我们需要求出表达式（11.2.7）关于v 的偏导数．过程如下： )),(),,(,(),(v x w v x w x v f v x v w '??=?''? ),()),(),,(,()),(),,(,(v x v w v x w v x w x w f v x v x w v x w x x f ??'??+??'??= ) ,()),(),,(,(v x v w v x w v x w x w f ?'?''??+ 又因为x 跟v 相互独立，所以当b x a ≤≤上式如下；

数值分析_第五章_常微分方程数值解法

图５畅２ 令珔h ＝h λ，则y n ＋１＝１＋珔 h ＋１２珔h ２ ＋１６珔h ３＋１２４ 珔 h ４y n ．由此可知，绝对稳定性区域在珔h ＝h λ复平面上满足 ｜１＋珔 h ＋１２珔h ２＋１６珔h ３＋１２４珔h ４ ｜≤１的区域，也就是由曲线 １＋珔h ＋ １２珔h ２＋１６珔h ３＋１２４ 珔h ４＝e i θ 所围成的区域．如图５畅２所示． 例22　用Euler 法求解 y ′＝－５y ＋x ，y （x ０）＝y ０， 　x ０≤x ≤X ． 从绝对稳定性考虑，对步长h 有何限制？ 解　对于模型方程y ′＝λy （λ＜０为实数）这里λ＝抄f 抄y ＝－５．由 ｜１＋h λ｜＝｜１－５h ｜＜１ 得到对h 的限制为：０＜h ＜０畅４． 四、习题 １畅取步长h ＝０畅２，用Euler 法解初值问题 y ′＝－y －x y ２ ， y （０）＝ １． 　（０≤x ≤０畅６）， ２畅用梯形公式解初值问题 y ′＝８－３y ， 　（１≤x ≤２），

取步长h＝０畅２，小数点后至少保留５位． ３畅用改进的Euler公式计算初值问题 y′＝１x y－１x y２， y（１）＝０畅５， 　１＜x＜１畅５， 取步长h＝０畅１，并与精确解y（x）＝ x １＋x比较． ４畅写出用梯形格式的迭代算法求解初值问题 y′＋y＝０， y（０）＝１ 的计算公式，取步长h＝０畅１，并求y（０畅２）的近似值，要求迭代误差不超过１０－５． ５畅写出用四阶经典Runge唱Kutta法求解初值问题 y′＝８－３y， y（０）＝２ 的计算公式，取步长h＝０畅２，并计算y（０畅４）的近似值，小数点后至少保留４位． ６畅证明公式 y n＋１＝y n＋h９（２K１＋３K２＋４K３）． K１＝f（x n，y n）， K２＝f x n＋h２，y n＋h２K１， K３＝f x n＋３４h，y n＋３４h K２， 至少是三阶方法． ７畅试构造形如 y n＋１＝α（y n＋y n－１）＋h（β０f n＋β１f n－１）

一阶常微分方程的解法

一阶常微分方程的解法 摘要：常微分方程是微积分学的重要组成部分，广泛用于具体问题的研究中，在整个数学中占有重要的地位。本文对一阶常微分方程的解法作了简要的总结，并举例加以分析了变量可分离方程，线性微分方程，积分因子，恰当微分方程，主要归纳了一阶微分方程的初等解法，并以典型例题加以说明。 关键词：变量分离;积分因子;非齐次微分方程;常数变易法 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 =

常微分方程组(边值)

常微分方程组边值问题解法 打靶法Shooti ng Method (shoot in g.m ) % 丁靶法求常微分方程的边值问题 function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3 eps=1e-3; end x1=x0+ra nd; [a,b]=ode45(fu n, [0,10],[0,x0]'); c0=b(le ngth(b),1); [a,b]=ode45(fu n, [0,10],[0,x1]'); c1=b(le ngth(b),1); x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=1; while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2; [a,b]=ode45(fu n,[ 0,10],[0,x0]'); cO=b(le ngth(b),1); [a,b]=ode45(fu n,[ 0,10],[0,x1]'); c1= b(le ngth(b),1) x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=n+1; end x=x2; 应用打靶法求解下列边值问题: y 10 0 解：将其转化为常微分方程组的初值问题

命令: xO=[O:O.1:1O]; y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'r') hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1)) dy i dx y 2 dy 2 dx y i 0 y 4 y o dy dx X0 真实解 30 ' 12^4567^9 10

偏微分方程数值解法

《偏微分方程数值解法》 课程设计 题目: 六点对称差分格式解热传导方程的初边 值问题 姓名: 王晓霜 学院: 理学院 专业: 信息与计算科学 班级: ０91１01２ 学号: ０９１10１218 指导老师：翟方曼 ２０12年12月14

日 一、题目 用六点对称差分格式计算如下热传导方程的初边值问题 222122,01,01(,0),01 (0,),(1,),01x t t u u x t t x u x e x u t e u t e t +???=<<<≤?????=≤≤??==≤≤??? 已知其精确解为 2(,)x t u x t e += 二、理论 1.考虑的问题 考虑一维模型热传导方程 （1.１) )(22x f x u a t u +??=??,T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.１)的定解问题分两类: 第一,初值问题（Cauch ｙ 问题)：求足够光滑的函数()t x u ,,满足方程（１.1）和初始条件： (1.２) ()()x x u ?=0,， ∞<<∞-x 第二,初边值问题（也称混合问题）:求足够光滑的函数()t x u ,，满足方程(１.１)和初始条件： ()13.1 ()()x x u ?=0,， l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑,并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 现在考虑边值问题（１.1),（1.３)的差分逼近 取 N l h = 为空间步长,M T =τ为时间步长,其中N ,M 是自然数,

matlab常微分方程和常微分方程组的求解

下载的，感觉不错，共享一下 常微分方程和常微分方程组的求解 一、实验目的： 熟悉Matlab 软件中关于求解常微分方程和常微分方程组的各种命令，掌握利用Matlab 软件进行常微分方程和常微分方程组的求解。 二、相关知识 在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 例1：求解常微分方程1dy dx x y = +的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')， 注意，系统缺省的自变量为t ，因此这里要把自变量写明。 结果为：-lambertw(-C1*exp(-x-1))-x-1 其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。 例2：求解常微分方程2 '''0yy y -=的MATLAB 程序为： Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为： Y2 =[ exp((x+C2)/C1)] [ C2] 我们看到有两个解，其中一个是常数。 例3：求常微分方程组253t t dx x y e dt dy x y e dt ?++=??? ?--=??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t') 例4：求常微分方程组020 210cos ,224,0 t t t dx dy x t x dt dt dx dy y e y dt dt =-=?+-==??? ?++==??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')

常微分方程组的MATLAB求解范例

微分方程求解是系统仿真、数学模型实现以及很多工程问题求解的核心部分，应用MATLAB可以方便地对一阶常微分方程组进行求解，这里将对其基本方法进行介绍。值得注意的是，高阶微分方程组可以通过引进参变量化为一阶常微分方程组，也可以同样方便解决。 若有一个微分方程（组）的参变量为列向量,即,且它参变量随时间变化的微分方程可以有以下方程描述： 这里的f函数是一个列向量，即, i=1,2,3…,n，它可以是任意非线性函数。 则一般微分方程可以如此求解： [t,x]=ode45(f,timespan,x0) 对于刚性方程，即一些解变化缓慢，一些解变化剧烈，且两者相差较为悬殊的这种方程，通常调用ode15s而非o de45进行求解。 例1： 解：编写function或者用匿名函数 表达f=y-2*x/y即可； function dy=f(t,y) dy=y-2*t/y; end 命令： t=[0,1];%y0=1; [x,y]=ode45('f',t,1);%注意 这里的x相当于自变量t plot(x,y,x,sqrt(1+2*x)),legend('数值解','解析解');

可见求解效果不错。 例2、 解：编写function function dx=f(t,x)%返回值是列向量 dx=[-x(2)-x(3); x(1)+0.2*x(2); 0.2+(x(1)-5.7)*x(3)]; end 命令： t=[0,100]; y0=[0 0 0]';%注意是列向量 [x,y]=ode45('f',t,y0); plot(x,y)； 例3、 这是一个二阶微分方程组，可以引进变量,由此ODE可以化成如下形式 可以采用和例2相同的方法求解： function dx=f(t,x) dx=[x(2); -(x(1)^2-1)*x(2)-x(1)]; End

一阶常微分方程解法归纳

第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(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 =+为变量可分离方程，得到)(0),,(为常数C C x u f =再把u 代入得到)(0),,(为常数C C x x y f =。 ②、形如)0(),(≠+=ab by ax G dx dy 解法：令by ax u +=，则b du adx dy +=，代入得到)(1u G b a dx du b =+为变量可分离方程， 得到)(0),,(为常数C C x u f =再把u 代入得到)(0),,(为常数C C x by ax f =+。 ③、形如 )(2 221 11c y b x a c y b x a f dx dy ++++= 解法：01、 02 2 11=b a b a ，转化为 )(by ax G dx dy +=，下同①； 02、 022 1 1≠b a b a ，???=++=++00 222111 c y b x a c y b x a 的解为),(00y x ，令???-=-=00y y v x x u 得到，)()( )(221 12211u v g u v b a u v b a f v b u a v b u a f du dv =++=++=，下同②； 还有几类：xy u dy xy xg dx xy yf ==+,0)()( xy v xy f dx dy x ==),(2 22),(x y w x y xf dx dy == θθsin ,cos ,0))(,())(,(r y r x ydx xdy y x N ydy xdx y x M ===-++ 以上都可以化为变量可分离方程。 例2.1、 2 5 --+-=y x y x dx dy 解：令2--=y x u ，则du dx dy -=，代入得到u u dx du 7 1+= - ，有dx udu 7-= 所以)(72 2 为常数C C x u +-=，把u 代入得到)(72 22 为常数） （C C x y x =+--。 例2.2、 1 212+-+-=y x y x dx dy

(整理)二阶常系数线性微分方程的解法word版.

第八章 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 为常数)和它的各阶导数都只差一个常数因子,

Matlab求解常微分方程边值问题的方法

Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2];

(完整版)二阶常微分方程边值问题的数值解法毕业论文

二阶常微分方程边值问题的数值解法 摘要 求解微分方程数值解的方法是多种多样的，它本身已形成一个独立的研究方向，其要点是对微分方程定解问题进行离散化．本文以研究二阶常微分方程边值问题的数值解法为目标，综合所学相关知识和二阶常微分方程的相关理论，通过对此类方程的数值解法的研究，系统的复习并进一步加深对二阶常微分方成的数值解法的理解，为下一步更加深入的学习和研究奠定基础． 对于二阶常微分方程的边值问题，我们总结了两种常用的数值方法：打靶法和有限差分法．在本文中我们主要探讨关于有限差分法的数值解法．构造差分格式主要有两种途径：基于数值积分的构造方法和基于Taylor展开的构造方法．后一种更为灵活，它在构造差分格式的同时还可以得到关于截断误差的估计．在本文中对差分方法列出了详细的计算步骤和Matlab

程序代码，通过具体的算例对这种方法的优缺点进行了细致的比较．在第一章中，本文将系统地介绍二阶常微分方程和差分法的一些背景材料．在第二章中，本文将通过Taylor展开分别求得二阶常微分方程边值问题数值解的差分格式．在第三章中，在第二章的基础上利用Matlab求解具体算例，并进行误差分析． 关键词：常微分方程,边值问题，差分法，Taylor展开，数值解

The Numerical Solutions of Second-Order Ordinary Differential Equations with the Boundary Value Problems ABSTRACT The numerical solutions for solving differential equations are various. It formed an independent research branch. The key point is the discretization of the definite solution problems of differential equations. The goal of this paper is the numerical methods for solving second-order ordinary differential equations with the boundary value problems. This paper introduces the mathematics knowledge with the theory of finite difference. Through solving the problems, reviewing what have been learned systematically and understanding the ideas and methods of the finite difference method in a deeper layer, we can establish a foundation for the future learning.