常微分方程的解法

-293-

第十五章 常微分方程的解法

建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如

22x y dx

dy

+=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。

§1 常微分方程的离散化

下面主要讨论一阶常微分方程的初值问题，其一般形式是

?????=≤≤=0

)()

,(y a y b

x a y x f dx

dy

（1）

在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条

件，即存在常数L ，使得

|||),(),(|y y L y x f y x f ?≤?

这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。

所谓数值解法，就是求问题（1）的解)(x y 在若干点

b x x x x a N =<<<<=L 210

处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解，

n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数

若用向前差商h

x y x y n n )

()(1?+代替)('n x y 代入（1）中的微分方程，则得

)1,,1,0())(,()

()(1?=≈?+N n x y x f h

x y x y n n n n L

化简得

))(,()()(1n n n n x y x hf x y x y +≈+

如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ，

则有

)1,,1,0()

,(1?=+=+N n y x hf y y n n n n L （2）

这样，问题（1）的近似解可通过求解下述问题 ??

?=?=+=+)

()

1,,1,0(),(01a y y N n y x hf y y n n n n L （3）

得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。

-294-

需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 （ii ）用数值积分方法

将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得

∫

+?==?+1

)1,,1,0())(,()()(1n n

x x n n N n dx x y x f x y x y L （4）

右边的积分用矩形公式或梯形公式计算。

（iii ）Taylor 多项式近似

将函数)(x y 在n x 处展开，取一次Taylor 多项式近似，则得

))(,()()(')()(1n n n n n n x y x hf x y x hy x y x y +=+≈+

再将)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值1+n y ，得到离

散化的计算公式

),(1n n n n y x hf y y +=+

以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。

§2 欧拉（Euler ）方法

2.1 Euler 方法

Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出)(n x y 的近似值)1,,2,1(?=N n y n L 。这组公式求问题（1）的数值解称为向前Euler 公式。

如果在微分方程离散化时，用向后差商代替导数，即h

x y x y x y n n n )

()()('11?≈++，

则得计算公式

??

?=?=+=+++)()

1,,1,0(),(0

111a y y N n y x hf y y n n n n L （5）

用这组公式求问题（1）的数值解称为向后Euler 公式。

向后Euler 法与Euler 法形式上相似，但实际计算时却复杂得多。向前Euler 公式是显式的，可直接求解。向后Euler 公式的右端含有1+n y ，因此是隐式公式，一般要用迭代法求解，迭代公式通常为

?????=+=+=+++++)

,2,1,0()

,(),()(11)1(1)0(1L k y x hf y y y x hf y y k n n n k n n n n n （6）

2.2 Euler 方法的误差估计

对于向前Euler 公式（3）我们看到，当1,,2,1?=N n L 时公式右端的n y 都是近似的，所以用它计算的1+n y 会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。

假定用（3）式时右端的n y 没有误差，即)(n n x y y =，那么由此算出

))(,()(1n n n n x y x hf x y y +=+ （7）

-295-

局部截断误差指的是，按（7）式计算由n x 到1+n x 这一步的计算值1+n y 与精确值)(1+n x y 之差11)(++?n n y x y 。为了估计它，由Taylor 展开得到的精确值)(1+n x y 是

)()(''2

)(')()(32

1h O x y h x hy x y x y n n n n +++=+ （8）

（7）、（8）两式相减（注意到),('y x f y =）得

)()()(''2

)(232

1

1h O h O x y h y x y n n n ≈+=?++ （9） 即局部截断误差是2

h 阶的，而数值算法的精度定义为：

若一种算法的局部截断误差为)(1

+p h O ，则称该算法具有p 阶精度。

显然p 越大，方法的精度越高。式（9）说明，向前Euler 方法是一阶方法，因此它的精度不高。

§3 改进的Euler 方法

3.1 梯形公式

利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分，即

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

))(,(111

+++≈∫+n n n n x x x y x f x y x f h

dx x y x f n n

并用1,+n n y y 代替)(),(1+n n x y x y ，则得计算公式

)],(),([2

111+++++=n n n n n n y x f y x f h

y y

这就是求解初值问题（1）的梯形公式。

直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为

??

???=++=+=+++++)

,2,1,0( )],,(),([2),()

(11)1(1)0(1

L k y x f y x f h y y y x hf y y k n n n n n k n n n n n （10） 由于函数),(y x f 关于y 满足Lipschitz 条件，容易看出

||2

||)

1(1)(1)(1)1(1?+++++?≤?k n k n k n k n y y hL y y

其中L 为Lipschitz 常数。因此，当12

0<

时，迭代收敛。但这样做计算量较大。如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler 法。

3.2 改进Euler 法

按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将Euler 公式与梯形公式结合使用：先用Euler 公式求1+n y 的一个初步近似值1+n y ，称为预测值，然后用梯形公式校正求得近似值1+n y ，即

-296-

??

?

??++=+=++++校正预测 ] ),(),([2)

,(1111n n n n n n n n n n y x f y x f h

y y y x hf y y （11） 式（11）称为由Euler 公式和梯形公式得到的预测—校正系统，也叫改进Euler 法。

为便于编制程序上机，式（11）常改写成

??

??

???

+=++=+=+)(2

1),(),(1q p n p n n q n n n p y y y y h x hf y y y x hf y y （12）

改进Euler 法是二阶方法。

§4 龙格—库塔（Runge —Kutta ）方法

回到Euler 方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理应有

10),(')

()(1<<+=?+θθh x y h

x y x y n n n

注意到方程),('y x f y =就有

))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （13）

不妨记))(,(h x y h x f K n n θθ++=，称为区间],[1+n n x x 上的平均斜率。可见给出一种斜率K ，（13）式就对应地导出一种算法。

向前Euler 公式简单地取),(n n y x f 为K ，精度自然很低。改进的Euler 公式可理解为K 取),(n n y x f ，),(11++n n y x f 的平均值，其中),(1n n n n y x hf y y +=+，这种处理提高了精度。

如上分析启示我们，在区间],[1+n n x x 内多取几个点，将它们的斜率加权平均作为

K ，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。

4.1 首先不妨在区间],[1+n n x x 内仍取2个点，仿照（13）式用以下形式试一下

???

??<<++==++=+1

,0),,()

,()(12

122111βαβαλλhk y h x f k y x f k k k h y y n n n n n n （14） 其中βαλλ,,,21为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差11)(++?n n y x y ，因为)(n n x y y =，所以（14）可以化为

???

????+++=++===++=+)())(,())(,())(,( )

)(,()

('))(,()

()(2112

1

22111h O x y x f hk x y x hf x y x f hk x y h x f k x y x y x f k k k h x y y n n y n n x n n n n n n n n n βαβαλλ （15）

-297-

其中2k 在点))(,(n n x y x 作了Taylor 展开。（15）式又可表为

)()()(')()(322211h O ff f h x hy x y y y x n n n ++

+++=+α

β

αλλλ 注意到

)()(''2

)(')()(32

1h O x y h x hy x y x y n n n n +++=+

中f y ='，y x ff f y +=''，可见为使误差)()(3

11h O y x y n n =?++，只须令

121=+λλ，212=

αλ，1=α

β

（16） 待定系数满足（16）的（15）式称为2阶龙格—库塔公式。由于（16）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令2

1

21=

=λλ ，1==βα，即为改进的Euler 公式。可以证明，在],[1+n n x x 内只取2点的龙格—库塔公式精度最高为2阶。

4.2 4阶龙格—库塔公式

要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：

?????

????++++=+++=++==++++=+)

,()

,()

,(),()(362514342

3122311121443322111hk hk hk y h x f k hk hk y h x f k hk 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 βββαββαβαλλλλ （17） 其中待定系数m j i βαλ,,共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差)()(5

11h O y x y n n =?++的11个方程。取既满足这些方程、又较简单的一组m j i βαλ,,，可得

????

?

??

???

???++=++=++==+++=+)

,()2

,2()

2,2(),()22(63423

12143211hk y h x f k hk y h x f k hk y h x f k y x f k k k k k h y n n n n n n n n n （18） 这就是常用的4阶龙格—库塔方法（简称RK 方法）。

§5 线性多步法

以上所介绍的各种数值解法都是单步法，这是因为它们在计算1+n y 时，都只用到前一步的值n y ，单步法的一般形式是

-298-

)1,,1,0()

,,(1?=+=+N n h y x h y y n n n n L ? （19）

其中),,(h y x ?称为增量函数，例如Euler 方法的增量函数为),(y x f ，改进Euler 法的

增量函数为

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

1),,(y x hf y h x f y x f h y x +++=?

如何通过较多地利用前面的已知信息，如r n n n y y y ??,,,1L ，来构造高精度的算法计算1+n y ，这就是多步法的基本思想。经常使用的是线性多步法。

让我们试验一下1=r ，即利用1,?n n y y 计算1+n y 的情况。 从用数值积分方法离散化方程的（4）式 ∫

+=

?+1

))(,()()(1n n

x x n n dx x y x f x y x y

出发，记n n n f y x f =),(，111),(???=n n n f y x f ，式中被积函数))(,(x y x f 用二节点

),(11??n n f x ，),(n n f x 的插值公式得到（因)n x x ≥，所以是外插。

])()[(1

))(,(1111

11?????????=

??+??=n n n n n

n n

n n n n n f x x f x x h

x x x x f x x x x f x y x f （20）

此式在区间],[1+n n x x 上积分可得 12

23))(,(1

??=

∫

+n n x x f h

f h dx x y x f n n

于是得到

)3(2

11?+?+

=n n n n f f h

y y （21） 注意到插值公式（20）的误差项含因子))((1n n x x x x ???，在区间],[1+n n x x 上积分后

得出3

h ，故公式（21）的局部截断误差为)(3

h O ，精度比向前Euler 公式提高1阶。

若取L ,3,2=r ，可以用类似的方法推导公式，如对于3=r 有

)9375955(243211???+?+?+=n n n n n n f f f f h

y y （22） 其局部截断误差为)(5

h O 。

如果将上面代替被积函数))(,(x y x f 用的插值公式由外插改为内插，可进一步减小误差。内插法用的是11,,,+?+r n n n y y y L ，取1=r 时得到的是梯形公式，取3=r 时

可得

)5199(24

2111??+++?++=n n n n n n f f f f h

y y （23） 与（22）式相比，虽然其局部截断误差仍为)(5

h O ，但因它的各项系数（绝对值）大为减小，误差还是小了。当然，（23）式右端的1+n f 未知，需要如同向后Euler 公式一

样，用迭代或校正的办法处理。

-299-

§6 一阶微分方程组与高阶微分方程的数值解法

6.1 一阶微分方程组的数值解法 设有一阶微分方程组的初值问题

??

?==0

21)()

,,,,('i i m i i y a y y y y x f y L ),,2,1(m i L = （24）

若记T

m y y y y ),,,(21L =，T

m y y y y ),,,(020100L =，T

m f f f f ),,,(21L =，则初值问题（24）可写成如下向量形式

??

?==0

)()

,('y a y y x f y （25） 如果向量函数),(y x f 在区域D ：

m R y b x a ∈≤≤,

连续，且关于y 满足Lipschitz 条件，即存在0>L ，使得],[b a x ∈?，m

R y y ∈21,，都有

2121),(),(y y L y x f y x f ?≤?

那么问题（25）在],[b a 上存在唯一解)(x y y =。

问题（25）与（1）形式上完全相同，故对初值问题（1）所建立的各种数值解法可全部用于求解问题（25）。

6.2 高阶微分方程的数值解法

高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。 设有m 阶常微分方程初值问题

???===≤≤=???)

1(0)1()1(00)1()()(,,)(',)(),,',,(m m m m y a y

y a y y a y b

x a y y y x f y L L （26） 引入新变量)

1(21,,',?===m m y y y y y y L ，问题（26）就化为一阶微分方程组初值问

题

?????

????========????)1(0

1)

2(011

)1(02320

121)( ),,,(')( ' )( ')( 'm m m m m m m m y a y y y x f y y a y y y y a y y y y a y y y L M

M （27） 然后用6.1中的数值方法求解问题（27），就可以得到问题（26）的数值解。

最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组

)(x Ay dx

dy

Φ+= （28） 其中A R y m

,,∈Φ为m 阶方阵。若矩阵A 的特征值),,2,1(m i i L =λ满足关系 ),,2,1(0Re m i i L =<λ

-300-

|Re |min |Re |max 11i m

i i m

i λλ≤≤≤≤>>

则称方程组（28）为刚性方程组或Stiff 方程组，称数

|Re |min /|Re |max 11i m

i i m

i s λλ≤≤≤≤=

为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长h 不作任何限制。

§7 初值问题的Matlab 解法和符号解

7.1 Matlab 数值解

7.1.1 非刚性常微分方程的解法

Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113，其中ode45采用四五阶RK 方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK 方法，ode113采用的是多步法，效率一般比ode45高。

Matlab 的工具箱中没有Euler 方法的功能函数。 （I ）对简单的一阶方程的初值问题

???==0

0)()

,('y x y y x f y

改进的Euler 公式为

??

??

???

+=++=+=+)(2

1),(),(1q p n p n n q n n n p y y y y h x hf y y y x hf y y

我们自己编写改进的Euler 方法函数eulerpro.m 如下：

function [x,y]=eulerpro(fun,x0,xfinal,y0,n); if nargin<5,n=50;end h=(xfinal-x0)/n; x(1)=x0;y(1)=y0; 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

例1 用改进的Euler 方法求解 x x y y 222'2

++?=，)5.00(≤≤x ，1)0(=y 解 编写函数文件doty.m 如下： function f=doty(x,y); f=-2*y+2*x^2+2*x;

在Matlab命令窗口输入：

[x,y]=eulerpro('doty',0,0.5,1,10) 即可求得数值解。

-301-

（II）ode23，ode45，ode113的使用 Matlab的函数形式是

[t,y]=solver('F',tspan，y0)

这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程),('y x f y =右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。

例2 用RK方法求解

x x y y 222'2++?=，)5.00(≤≤x ，1)0(=y

解 同样地编写函数文件doty.m 如下： function f=doty(x,y); f=-2*y+2*x^2+2*x;

在Matlab命令窗口输入：

[x,y]=ode45('doty',[0,0.5],1) 即可求得数值解。

也可以利用Matlab的匿名函数，计算数值解，程序如下： f=@(x,y)-2*y+2*x^2+2*x; [x,y]=ode45(f,[0,0.5],1)

7.1.2 刚性常微分方程的解法

Matlab 的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s ，ode23s ，ode23t ，ode23tb ，这些函数的使用同上述非刚性微分方程的功能函数。 7.1.3 高阶微分方程),,',,()

1()

(?=n n y y y t f y L 解法 例3 考虑初值问题

1)0(''1)0('0)0(0'''3'''?====??y y y y y y y

解 （i ）如果设'',',321y y y y y y ===，那么

???

???=+=====1

)0( 3'1)0( '0)0( '31233

232121y y y y y y y y y y y

初值问题可以写成0)0(),,('Y Y Y t F Y ==的形式，其中T

y y y Y ][321=。

（ii ）把一阶方程组写成接受两个参数t 和y ，返回一个列向量的M 文件F.m ：

function dy=F(t,y);

dy=[y(2);y(3);3*y(3)+y(2)*y(1)];

注意：尽管不一定用到参数t 和y ，M —文件必须接受此两参数。这里向量dy 必须是列向量。

（iii）用Matlab 解决此问题的函数形式为

[T,Y]=solver('F',tspan,y0)

这里solver 为ode45、ode23、ode113，输入参数F 是用M 文件定义的常微分方程组，tspan=[t0 tfinal]是求解区间，y0是初值列向量。在Matlab 命令窗口输入

[T,Y]=ode45('F',[0 1],[0;1;-1])

就得到上述常微分方程的数值解。这里Y 和时刻T 是一一对应的，Y(:,1)是初值问题的解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。

也可以利用Matlab 的匿名函数，求数值解，编写如下程序： F=@(t,y)[y(2);y(3);3*y(3)+y(2)*y(1)]; [T,Y]=ode45(F,[0 1],[0;1;-1])

-302-

例4 求 van der Pol 方程 0')1(''2

=+??y y y y μ 的数值解，这里0>μ是一参数。

解 （i）化成常微分方程组。设',21y y y y ==，则有

?????==1

22

122

1)1(''y y y y y y μ （ii ）书写M 文件（对于1=μ）vdp1.m:

function dy=vdp1(t,y);

dy=[y(2);(1-y(1)^2)*y(2)-y(1)];

（iii ）调用Matlab 函数。对于初值0)0(',2)0(==y y ，解为 [T,Y]=ode45('vdp1',[0 20],[2;0]);

（iv ）观察结果。利用图形输出解的结果： plot(T,Y(:,1),'-',T,Y(:,2),'--')

title('Solution of van der Pol Equation,mu=1'); xlabel('time t');

ylabel('solution y');

例5 van der Pol 方程，1000=μ（刚性） 解 (i)书写M 文件vdp1000.m: function dy=vdp1000(t,y);

dy=[y(2);1000*(1-y(1)^2)*y(2)-y(1)];

（ii）观察结果

[t,y]=ode15s('vdp1000',[0 3000],[2;0]); plot(t,y(:,1),'o')

title('Solution of van der Pol Equation,mu=1000'); xlabel('time t');

ylabel('solution y(:,1)');

7.2 常微分方程的解析解

在Matlab 中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令dsolve 。常微分方程在Matlab 中按如下规定重新表达：

-303-

符号D 表示对变量的求导。Dy 表示对变量y 求一阶导数，当需要求变量的n 阶导数时，用Dn 表示，D4y 表示对变量y 求4阶导数。

由此，常微分方程y y y =+'2''在Matlab 中，将写成'D2y+2*Dy=y'。 7.2.1 求解常微分方程的通解

无初边值条件的常微分方程的解就是该方程的通解。其使用格式为： dsolve('diff_equation') dsolve(' diff_equation'，'var')

式中diff_equation 为待解的常微分方程，第1种格式将以变量t 为自变量进行求解，第2种格式则需定义自变量var 。

例6 试解常微分方程 0')2(2

=?++y y x y x 解 编写程序如下： syms x y

diff_equ='x^2+y+(x-2*y)*Dy=0'; dsolve(diff_equ,'x')

7.2.2 求解常微分方程的初边值问题

求解带有初边值条件的常微分方程的使用格式为：

dsolve('diff_equation'，'condition1，condition2，…'，'var') 其中condition1，condition2，… 即为微分方程的初边值条件。

例7 试求微分方程

x y y =?'''''，4)2('',7)1(',8)1(===y y y

的解。

解 编写程序如下：

y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x') 7.2.3 求解常微分方程组

求解常微分方程组的命令格式为：

dsolve('diff_equ1，diff_equ2，…'，'var')

dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')

第1种格式用于求解方程组的通解，第2种格式可以加上初边值条件，用于具体求解。

例8 试求常微分方程组：

?

?

?=+=+x f g x

g f cos ''sin 3'' 的通解和在初边值条件为1)5(,3)3(,0)2('===g f f 的解。

解 编写程序如下： clc,clear

equ1='D2f+3*g=sin(x)'; equ2='Dg+Df=cos(x)';

[general_f,general_g]=dsolve(equ1,equ2,'x')

[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

7.2.4 求解线性常微分方程组 （i）一阶齐次线性微分方程组

-304-

AX X ='，??????????=n x x X M 1，?????

?????=nn n n a a a a A L M M M L 1111

这里的’表示对t 求导数。At

e 是它的基解矩阵。AX X ='，00)(X t X =的解为

0)(0)(X e t X t t A ?=。

例9 试解初值问题

X X ???????????=200120312'，????

?

?????=121)0(X

解 编写程序如下：

syms t

a=[2,1,3;0,2,-1;0,0,2]; x0=[1;2;1]; x=expm(a*t)*x0

（ii）非齐次线性方程组

由参数变易法可求得初值问题

)('t f AX X +=，00)(X t X = 的解为

∫??+=t

t s t A t t A ds s f e X e

t X 00)()()(0)

(. 例10 试解初值问题

??????????+???????????=t e X X t

2cos 00123212001'，??

??

?

?????=110)0(X 。

解 编写程序如下：

clc,clear syms t s

a=[1,0,0;2,1,-2;3,2,1];ft=[0;0;exp(t)*cos(2*t)]; x0=[0;1;1];

x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t); x=simple(x)

§8 边值问题的Matlab 解法

Matlab 中用bvp4c 命令求解常微分方程的两点边值问题，微分方程的标准形式为 ),('y x f y =，0))(),((=b y a y bc 或者是

),,('p y x f y =，0)),(),((=p b y a y bc 其中p 是有关的参数。这里f y ,可以为向量函数，求解的区间为],[b a ，bc 为边界条件。

一般地说，边值问题在计算上比初值问题困难得多，特别地，由于边值问题的解可能是多值的，bvp4c 需要提供猜测的初始值。下面我们首先给出一个简单的例子。

-305-

例11 考察描述在水平面上一个小水滴横截面形状的标量方程

0)))((1))((1()(2

/3222=+?+x h dx d x h x h dx

d ，0)1()1(==?h h

这里)(x h 表示x 处水滴的高度。设)()(1x h x y =，dx

x dh x y )

()(2=，把上述微分方程写成两个一阶微分方程组

)()(21x y x y dx d

= 2/32212))(1)(1)(()(x y x y x y dx

d

+?= 上述微分方程组可以由如下函数表示 function yprime=drop(x,y);

yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];

边界条件通过残差函数指定，边界条件通过如下函数表示 function res=dropbc(ya,yb); res=[ya(1);yb(1)];

我们使用211)(x x y ?=和)11.0/()(2

2x x x y ?+?=作为初始猜测解，由如下函数定义

function yinit=dropinit(x);

yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];

利用如下的程序就可以求微分方程的边值问题并画出图2。 solinit=bvpinit(linspace(-1,1,20),@dropinit); sol=bvp4c(@drop,@dropbc,solinit); fill(sol.x,sol.y(1,:),[0.7,0.7,0.7]) axis([-1,1,0,1])

xlabel('x','FontSize',12)

ylabel('h','Rotation',0,'FontSize',12)

x

h

图2

这里调用函数bvpinit ，计算区间]1,1[?上等间距的20个点的数据，然后调用函数

bvp4c ，得到数值解的结构sol ，填充命令fill 填充1y x ?平面上的解曲线。

一般地，bvp4c 的调用格式如下：

-306-

sol=bvp4c(@odefun,@bcfun,solinit,options,p1,p2,…); 函数odefun 的格式为

yprime=odefun(x,y,p1,p2, …); 函数bcfun 的格式为

res=bcfun(ya,yb, p1,p2, …);

初始猜测结构solinit 有两个域：solinit.x 提供初始猜测的x 值，排列顺序从左到右排列，其中solinit.x(1)和solinit.x(end)分别为a 和b 。对应地，solinit.y(:,i)给出点solinit.x(i)处初始猜测解。

输出参数sol 是包含数值解的一个结构，其中sol.x 给出了计算数值解的x 点，sol.x(i)处的数值解由sol.y(:,i)给出，类似地，sol.x(i)处数值解的一阶导数值由sol.yp(:,i)给出。

我们可以把上面的所有函数都放在一个文件中，程序如下： function sol=example11;

solinit=bvpinit(linspace(-1,1,20),@dropinit); sol=bvp4c(@drop,@dropbc,solinit); fill(sol.x,sol.y(1,:),[0.7,0.7,0.7]) axis([-1,1,0,1])

xlabel('x','FontSize',12)

ylabel('h','Rotation',0,'FontSize',12)

function yprime=drop(x,y);

yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];

function res=dropbc(ya,yb); res=[ya(1);yb(1)];

function yinit=dropinit(x);

yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];

例12 描述0=x 处固定，1=x 处有弹性支持，沿着x 轴平衡位置以均匀角速度旋转的绳的位移方程

0)()(22

=+x y x y dx

d μ

具有边界条件 0)0(=y ，1))((

0==x x y dx d ，0))()((1

=+=x x y dx d x y 这个边值问题是一个特征问题，我们必须找到参数μ的值使得方程的解存在。如

果我们提供了参数μ的猜测值和对应解的猜测值，我们也可以利用函数bvp4c 求解特征问题。上述微分方程可以写成下面的微分方程组：

)()(21x y x y dx d

= )()(12x y x y dx

d

μ?= 编写程序如下（下面的所有程序放在一个文件中）： function sol=skiprun;

-307-

solinit=bvpinit(linspace(0,1,10),@skipinit,5); sol=bvp4c(@skip,@skipbc,solinit);

plot(sol.x,sol.y(1,:),'-',sol.x,sol.yp(1,:),'--','LineWidth',2) xlabel('x','FontSize',12) legend('y_1','y_2',0)

function yprime=skip(x,y,mu); yprime=[y(2);-mu*y(1)];

function res=skipbc(ya,yb,mu); res=[ya(1);ya(2)-1;yb(1)+yb(2)];

function yinit=skipinit(x); yinit=[sin(x);cos(x)];

图3给出上述边值问题解的图象。

x

图3

例13 微分方程组为

v u w u u /)(5.0'?=

)(5.0'u w v ??=

z u w w y w w /))(5.0)(10009.0('????= )(5.0'u w z ?= )(100'w y y ??=

边界条件为1)0()0()0(===w v u ，10)0(?=z ；)1()1(y w =。

我们使用如下猜测解 1)(=x u

1)(=x v

191.85.4)(2

++?=x x x w 10)(=x z

91.095.4)(2++?=x x x y

我们编写如下程序（所有程序放在一个文件中）： function sol=example13

solinit=bvpinit(linspace(0,1,5),@exinit); sol=bvp4c(@exode,@exbc,solinit); plot(sol.x,sol.y)

-308-

function yprime=exode(x,y)

yprime=[0.5*y(1)*(y(3)-y(1))/y(2) -0.5*(y(3)-y(1))

(0.9-1000*(y(3)-y(5))-0.5*y(3)*(y(3)-y(1)))/y(4) 0.5*(y(3)-y(1)) 100*(y(3)-y(5))];

function res=exbc(ya,yb)

res=[ya(1)-1;ya(2)-1;ya(3)-1;ya(4)+10;yb(3)-yb(5)];

function yinit=exinit(x)

yinit=[1;1;-4.5*x^2+8.91*x+1;-10;-4.5*x^2+9*x+0.91];

习 题 十 五

1. 用欧拉方法和龙格—库塔方法求微分方程数值解，画出解的图形，对结果进行分析比较。

（i ）πππ22',22,0)('''2

22?=??

?

???=????

??=?++y y y n x xy y x （Bessel 方程，令

21=n ，精确解x

x

y π

2sin =。 （ii ）0)0(',1)0(,0cos ''===+y y x y y ，幂级数解

L ?+?+?=8

642!

855!69!42!211x x x x y

2. 一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。已知河水流速1v 与船在静水中的速度2v 之比为k 。

（i ）建立小船航线的方程，求其解析解。

（ii ）设100=d m ，11=v m/s ，22=v m/s ，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。

3. Lorenz 方程是一个三阶的非线性系统，它是由描述大气动力系统的Navier-Stokes偏微分方程演化而来的。自由系统如下：

???

??+?=??=?=xy

z z

xz y x y x y x λβσ&&&)( 当系统参数λβσ,,

在一定范围内，系统就出现混沌，如10=σ，28=β，

3/8=λ时，出现混沌现象。求在初始条件]17,13,5[)]0(),0(),0([=z y x 时，方程组

的数值解，并画出解的图形。

一阶常微分方程解法总结

第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(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 =+为变量可分离方程，得到

常微分方程解题方法总结.doc

常微分方程解题方法总结 来源：文都教育 复习过半， 课本上的知识点相信大部分考生已经学习过一遍 . 接下来， 如何将零散的知 识点有机地结合起来， 而不容易遗忘是大多数考生面临的问题 . 为了加强记忆， 使知识自成 体系，建议将知识点进行分类系统总结 . 著名数学家华罗庚的读书方法值得借鉴， 他强调读 书要“由薄到厚、由厚到薄”，对同学们的复习尤为重要 . 以常微分方程为例， 本部分内容涉及可分离变量、 一阶齐次、 一阶非齐次、 全微分方程、 高阶线性微分方程等内容， 在看完这部分内容会发现要掌握的解题方法太多， 遇到具体的题 目不知该如何下手， 这种情况往往是因为没有很好地总结和归纳解题方法 . 下面以表格的形 式将常微分方程中的解题方法加以总结，一目了然，便于记忆和查询 . 常微分方程 通解公式或解法 ( 名称、形式 ) 当 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]

常微分方程数值解法的误差分析教材

淮北师范大学 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 。

各类微分方程的解法大全

各类微分方程的解法 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个系数

常微分方程数值解法

第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==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 阶精度。

常微分方程数值解法

常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步： 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(非刚性)常微分方程。

15第十五章 常微分方程的解法

-293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 22x y dx dy +=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ?????=≤≤=0 )() ,(y a y b x a y x f dx dy （1） 在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条 件，即存在常数L ，使得 |||),(),(|y y L y x f y x f ?≤? 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解)(x y 在若干点 b x x x x a N =<<<<=L 210 处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解， n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商h x y x y n n ) ()(1?+代替)('n x y 代入（1）中的微分方程，则得 )1,,1,0())(,() ()(1?=≈?+N n x y x f h x y x y n n n n L 化简得 ))(,()()(1n n n n x y x hf x y x y +≈+ 如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ， 则有 )1,,1,0() ,(1?=+=+N n y x hf y y n n n n L （2） 这样，问题（1）的近似解可通过求解下述问题 ?? ?=?=+=+) () 1,,1,0(),(01a y y N n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。

(整理)常微分方程数值解法

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 -=<<<<=L （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 +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。

(完整版)专题一(二阶常微分方程解法)

二阶微分方程： 时为非齐次 时为齐次，0)(0)()()()(22≠≡=++x f x f x f y x Q dx dy x P dx y d 二阶常系数齐次线性微分方程及其解法： 2 122,)(2,,(*)0)(1,0(*)r r y y y r r q pr r q p qy y p y 式的两个根、求出的系数； 式中的系数及常数项恰好是，，其中、写出特征方程：求解步骤： 为常数； ，其中?'''=++?=+'+''式的通解：出的不同情况，按下表写、根据(*),321r r 二阶常系数非齐次线性微分方程 型 为常数； 型，为常数 ，]sin )(cos )([)()()(,)(x x P x x P e x f x P e x f q p x f qy y p y n l x m x ωωλλλ+===+'+'' 二阶常系数非齐次线性微分方程的一般形式是 ''+'+=y py qy f x () （1） 其中p q ,是常数。 方程（1）的通解为对应的齐次方程 0=+'+''qy y p y （2） 的通解Y 和方程（1）的一个特解*y 之和。即 *y Y y +=.我们已解决了求二阶常系数齐 次线性方程通解的问题，所以，我们只需讨论求二阶常系数非齐次线性微分方程的特解* y 的方法。 下面我们只介绍当方程（1）中的)(x f 为如下两种常见形式时求其特解*y 的方法。 一、 f x e P x x m ()()=?λ型 由于方程（1）右端函数f x ()是指数函数e x λ?与m 次多项式P x m ()的乘积，而指数

函数与多项式的乘积的导数仍是这类函数，因此，我们推测： 方程（1）的特解应为 y e Q x x *?=λ()( Q x ()是某个次数待定的多项式 ) y e Q x e Q x x x *??'=+'λλλ()() y e Q x Q x Q x x *?"=?+'+''λλλ[()()()]22 代入方程（1），得 e Q x p Q x p q Q x e P x x x m λλλλλ???''++'+++≡?[()()()()()]()22 消去e x λ?，得 ''++'+++≡Q x p Q x p q Q x P x m ()()()()()()22λλλ （3） 讨论 01、如果λ不是特征方程 r pr q 20++=的根。 即 02≠++q p λλ 由于P x m ()是一个m 次的多项式，欲使（3）的两端恒等，那未Q x ()必为一个m 次多项式，设为 Q x b x b x b x b m m m m m ()=++++--0111Λ 将之代入（3），比较恒等式两端x 的同次幂的系数，就得到以b b b b m m 01 1,,,,Λ-为未知数的m +1个线性方程的联立方程组，解此方程组可得到这m +1个待定的系数，并得到特解 y e Q x x m *?=λ() 02、如果λ是特征方程 r pr q 20++=的单根。 即 λλ20++=p q ，但 20λ+≠p 欲使（3）式的两端恒等，那么'Q x ()必是一个m 次多项式。 因此，可令 Q x x Q x m ()()=? 并且用同样的方法来确定)(x Q 的系数b b b b m m 0 11,,,,Λ-。 03、如果λ是特征方程 r pr q 20++=的二重根。 即 λλ20++=p q ，且 20λ+=p 。 欲使（3）式的两端恒等，那么''Q x ()必是一个m 次多项式 因此， 可令 Q x x Q x m ()()=?2 并且用同样的方法来确定)(x Q 的系数b b b b m m 011,,,,Λ-。

二阶常微分方程的解法及其应用

目录 1 引言 (1) 2 二阶常系数常微分方程的几种解法 (1) 2.1 特征方程法 (1) 2.1.1 特征根是两个实根的情形 (2) 2.1.2 特征根有重根的情形 (2) 2.2 常数变异法 (4) 2.3 拉普拉斯变化法 (5) 3 常微分方程的简单应用 (6) 3.1 特征方程法 (7) 3.2 常数变异法 (9) 3.3 拉普拉斯变化法 (10) 4 总结及意义 (11) 参考文献 (12)

二阶常微分方程的解法及其应用 摘要：本文通过对特征方程法、常数变易法、拉普拉斯变换法这三种二阶常系数常微分方程解法进行介绍，特别是其中的特征方程法分为特征根是两个实根的情形和特征根有重根的情形这两种情况，分别使用特征值法、常数变异法以及拉普拉斯变换法来求动力学方程，现今对于二阶常微分方程解法的研究已经取得了不少成就，尤其在二阶常系数线性微分方程的求解问题方面卓有成效。应用常微分方程理论已经取得了很大的成就，但是，它的现有理论也还远远不能满足需要，还有待于进一步的发展，使这门学科的理论更加完善。 关键词：二阶常微分方程;特征分析法；常数变异法；拉普拉斯变换

METHODS FOR TWO ORDER ORDINARY DIFFERENTIAL EQUATION AND ITS APPLICATION Abstract:This paper introduces the solution of the characteristic equation method, the method of variation of parameters, the Laplasse transform method the three kind of two order ordinary differential equations with constant coefficients, especially the characteristic equation method which is characteristic of the root is the two of two real roots and characteristics of root root, branch and don't use eigenvalue method, method of variation of constants and Laplasse transform method to obtain the dynamic equation, the current studies on solution of ordinary differential equations of order two has made many achievements, especially in the aspect of solving the problem of two order linear differential equation with constant coefficients very fruitful. Application of the theory of ordinary differential equations has made great achievements, however, the existing theory it is still far from meeting the need, needs further development, to make the discipline theory more perfect. Keywords:second ord er ordinary differential equation; Characteristic analysis; constant variation method; Laplasse transform 1 引言 数学发展的历史告诉我们，300年来数学分析是数学的首要分支，而微分方程

常微分方程常用数值解法.

第一章绪论 1.1 引言 常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。 研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。 由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。 本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。

常微分方程的初等解法与求解技巧

师大学本科毕业论文（设计） 常微分方程的初等解法与求解技巧 姓名娟 院系数学与计算机科学学院 专业信息与计算科学 班级12510201 学号1251020126 指导教师王晓锋 答辩日期 成绩

常微分方程的初等解法与求解技巧 容摘要 常微分方程在数学中发挥着举足轻重的作用，同时它的应用在日常生活里随处可见，因此掌握常微分方程的初等解法与求解技巧是非常必要的.本论文主要论述了其发展、初等解法与求解技巧，前者主要有变量分离、积分因子、一阶隐式微分方程的参数表示，通过举例从中总结出其求解技巧，目的是掌握其求解技巧. 【关键词】变量分离一阶隐式微分方程积分因子求解技巧

Elementary Solution and Solving Skills of Ordinary Differential Equation Abstract Ordinary differential equations take up significant position in mathematics, and at the same time, the application of it can be seen everywhere in our daily life, therefore, it’s necessary to grasp the elementary solution of ordinary differential equations and solving skills. This paper mainly introduced the definition of ordinary differential equations, elementary solution method and solving skills, the former mainly included the separation of variables, integral factor, a parameter-order differential equations implicit representation, by way of examples to sum up their solving skills, the purpose is to master the skills to solve. 【Key Words】the separation of variables the first order implicit differential equation integrating factor solution techniques

各类微分方程的解法大全

各类微分方程的解法大 全 TTA standardization office【TTA 5AB- TTAK 08- TTA 2C】

各类微分方程的解法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,C1) 即dy/dx=φ(y,C1),即dy/φ(y,C1)=dx,所以∫dy/φ(y,C1)=x+C2 5.二阶常系数齐次线性微分方程解法 一般形式:y”+py’+qy=0,特征方程r2+pr+q=0 6.二阶常系数非齐次线性微分方程解法 一般形式: y”+py’+qy=f(x) 先求y”+py’+qy=0的通解y0(x),再求y”+py’+qy=f(x)的一个特解y*(x) 则y(x)=y0(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个系数