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

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

工程技术与科学实验中提出的大量问题是常微分方程边值问题．本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法.

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 ≤≤上式如下；

),()),(),,(,(),()),(),,(,(),(v x v

w v x w v x w x w f v x v w v x w v x w x w f v x v w ?'

?''??+??'??=?''?（11.2.8）

初始条件为：

1),(0),(=?'

?=??v a v w v a v w ， 如果简单地用),(v x z 定义),)((v x v

w

??，并且假定x 和v 的微分阶翻转，（11.2.8）转化

为初值问题：

1

)(0)(),,(),,(....

==≤≤''??+'??=a z a z b x a z w w x w f

z w w x w f z ，，，，（11.2.9） 因此，对于（11.2.2）和（11.2.9）式的每次迭代需要求解两个初值问题．那么从（11.2.6）式可得：

),(),(111-----

=k k k k v b z v b w v v β

， （11.2.10）

事实上，这些初值问题很难精确求解，而将这些解近似为一个初值问题的解．

同样，我们可以按以上步骤考虑对于三阶非线性边值问题的打靶法算法． 对于三阶非线性边值问题

()()().)(,,,,γβα='='=≤≤'''='''b y a y a y b x a y y y x f y ，，，(11.2.11)

转变形式：

()()v a w a w a w b x a w w w x f w =''='=≤≤'''=''')()(,,,，，，，βα,(11.2.12)

选择参数k v v =，以使：

()()β=='∞

→b y v b w k k ,lim ， （11.2.13）

其中),(k v x w '定义为初值问题（11.2.12）在k v v =时的解，同时()x y 定义为边值问题（11.2.11）的解．

定义参数0v ，沿着如下初值问题解的曲线，可以求出点),(αa 对应的初始正视图

()()v a w a w a w b x a w w w x f w =''='=≤≤'''=''')()(,,,，，，，βα．（11.2.14）

如果),(0v b w '不严格收敛于β，那么我们选择1v 等值以修正近似值，直到),(0v b w '严

格逼近β．

为了取得合适的参数k v ，现在假定边值问题（11.2.11）有唯一解，如果),(v x w 定义为初始问题（11.2.12）的解，那么v 可由下式确定：

0),(=-'βv b w ． （11.2.15）

利用Newton 法求解．首先选择初始值0v ，然后由下式生成序列

),)(()),((111---'-'-

=k k k k v b dv

w d v b w v v β，此处),(),)((11--'

='k k v b dv w d v b dv w d ， （11.2.16）

同时要求求得),)((

1-'

k v b dv

w d ，因为),(v b w '的表达式未知，所以只能得到这么一系列的值。，，，),(),(),(),(1210-'??'''k v b w v b w v b w v b w

改写初值问题（11.2.12），使其强调解对x 和v 的依赖性

()(),(,),(,),(,),(,)(,)w f x w x v w x v w x v a x b w a v w a v w a v v

αβ''''''='''≤≤===，

，，， （11.2.17）

要求当k v v =时),)((

v b dv

w d '

的值，那么我们需要求出表达式（11.2.17）关于v 的偏导数．所以当b x a ≤≤上式如下；

v

w w f v w w f v w w f v x v w ?'

'?''??+

?'?'??+????=?'''?),( （11.2.18） 初始条件为：

1),(0),(0),(=?'

'?=?'?=??v a v

w v a v w v a v w ，， 对于第三种边界条件, 注意到 0

1

)(a a y'(a)a y -=

, 取 s a y =)('得到 ??

?

?

?=-==s a , y'a a a y'a y x,y,y'f y )()()()(01, b x a <<, (11.2.19) 设它的解为)(x,s y , 代入第二个边界条件得到

s βs F s b y βb,s y b y βb y ==-=-')() ,()(')()(00

即

s βs F =)(

当1s ββ= 时, 即为所求的解. 这样,同第一种边界条件一样, 可以利用打靶方法求解.

例11.1 利用打靶法求解两点边值问题

2()0

(0)0,(0.4) 1.8

y y x y y '''?-+=?

==? 解：(1) 为应用打靶方法，需要假定初值(0)y '来先求解初值问题，取初始参数

0 1.81

20.40t -==-， 解2()(0)0,(0)2y y x y y '''?=-?

'==?

令1y u =，2y u '=，则上述初值问题变为

1212

222,(0)1

,(0)2

u u u u u x u ?'==??

'=-=??， 由标准的R-K 方法，取0.2h =，可得0 2.540β=。 （2）再令11t =， 解

1212

222,(0)1

,(0)1

u u u u u x u ?'==??

'=-=??， 取0.2h =。 仍由标准的R-K 方法，可得1 1.497 1.8β=<。 （3）令212

1(1.8 1.497) 1.291.497 2.540

t -=+

-=-

再解1212

222,(0)1,(0) 1.29

u u u u u x u ?'==??'=-=??， 取0.11h =，得2 1.7094β=。

（4）令3(1.291)(1.8 1.7094)

1.29 1.4141.7094 1.497

t --=+=-，得3 1.81978β=

（5）类似有，4 1.38, 1.78879t β==，5 1.3960, 1.79995t β== 因为5 1.80.00005β-<， 可以认为 1.3960t =即为所求。

11.3 有限差分方法

有限差分方法常用来求解常微分方程边值问题, 下面就线性微分方程为例进行讨论. 设微分方程具有形式:

)()()()(')()(x f x y x q x y x p x 'y'=++, b x a <<, (11.3.1)

α=)(a y , β=)(b y . (11.3.2)

不失一般性, 将求解区间n 等份, 分点为a x =0, h a x +=1,…, ih a x i +=

),,1,0(n i =, 步长n

a

b h -=

, 这里 i x 称为结点. 在内部结点i x 上, 记 i i y x y =)(, i i p x p =)(, i i q x q =)(, i i f x f =)(, 并分别用一阶和二阶中心差分来代替一阶导数和二阶导数,

)(2)('1

12-i i i h O h

y y x y +-=

+,

)(2)(''2

1

12-i i i i h O h

y y y x y ++-=

+. 将上面各式代入方程(11.3.1), 略去)(2

h O 项, 得到关系式:

i i -i i i -i i i f h q h y y hp y y y 221111)(2

1

)2(=+-++-++,

将上式改写成

i i i i i -i i f h y c y a y b 2112

1

-=-+-+, )1,,1(-n i =, (11.3.3)

在上式中

i i q h a 212-=, )21(21i i p h b -=, )21(21i i p h

c +=,)1,,1(-n i =, (11.3.4)

利用边界条件 ，）α=a y (β=)(b y , 可以得到如下形式的方程组:

F AY = (11.3.5)

其中

????????

???

?

????????

??????---=---1123322211n n n a b c a b c a b c a A (11.3.6)

+?????????????

???=-1212 2n f f f h F ?????

??

?

????????-βα110 0 n c b (11.3.7) 系数矩阵(11.3.6)为n-1阶三对角阵, 当系数矩阵非奇异时,方程组(11.3.5) 有唯一解,可以利用追赶法求解.

对于第二种边界条件:

α=')(a y , β=')(b y (11.3.8)

和第三种边界条件:

?

?

?=-'=-'101

0)()()()(b b y b y a a y a y βα, (11.3.9) 其中0 0, ,00000>+≥≥b a b a .

由于边界条件中包含了导数, 故边界条件也必须用差商近似来表示. 因为不能利用区间][a,b 以外的点,所以在引进'

y 0(即 )(a y'或 '

n y (即 )(b y') 的近似式时,就不能利用中心差商. 如果只要求误差是)(h O , 则可以用最简单的近似公式

h y y y '

010-≈

, h

y y y n n '

n 1--≈ (11.3.10) 要使误差达到)(2

h O , 可以利用Newton 等距插值公式,得到如下近似公式:

h y y y y y h y '

2342110

120200-+-=??

????-≈

??, (11.3.11) h y y y y y h y n n n n n 'n 24321121221----+-=??

????+≈

??, (11.3.12) 将这些近似式代入边值条件中, 再和方程(11.3.4)联立, 就可以得到对应的差分方程组. 例11.2 利用差分法求解两点边值问题

2(1)1

(0)1,(1)3

y x y y y ''?-+=?

==? 解：将方程离散后写成矩阵形式

2.0101 1 1 2.0104 1 1 2.0164 1 ---12890.990.010.01 1 2.0181 2.99y y y y -????????????????????????????????????=??????????????????????????????--?????????

??? 取0.05h =， 可以得到129,,,y y y 三对角方程，下表11-1给出了计算结果。 表 11-1取0.05h =时计算结果

11.4 应用实例与Matlab

11.4.1 MATLAB 关于常微分方程边值问题数值解法的命令

常微分方程初、边值问题的符号解法 命令形式1：dsolve(，eqution ，) 功能：求常微分方程eqution 的解。

命令形式2：dsolve(，eqution ，，，condl ，cond2…，，，var ，) 功能：求常微分方程eqution 的满足初始条件的特解。其中eqution 是 求解的微分方程或微分方程组，condl ，cond2…是初始条件或边值，var 是自变

量.

例11.3求解两点边值问题：2

''3'(1)0(5)0

xy y x y y ?-=?==?.

解 Matlab 命令为

y=dsolve(…x ＊D2y-3＊Dy=x^2?，?y(1)=0，y(5)=0?,?x ?) 执行结果：

y ＝－1／3＊x^3+125／468+31／468＊x^4 11.4.2 应用实例

当流体以均匀的流速纵掠一平壁时，由于流体粘性的影响，在靠近壁面邻近会形成很大的速度梯度，这就是速度边界层．用来描述速度边界层的动量方程经过相似变换可以化为如下常微分方程两点边界值问题(著名的布拉修斯边界层方程)

1

02

f ff '''''+

= (11.3.13) 边界条件为

(0), (0)0, ()0f N f f ''=±=+∞= (11.3.14)

这里(0)f 表示壁面喷注(或吸入). 根据给定的喷注速度值，(0)f N =+(常数)表示吸入，

(0)f N =-(常数)表示喷注.

为了利用打靶法求解该问题的解，假设()0f β''=（这里β是一个待定常数，需要在计算中估计). 表11-1 给出了数值计算方法。当0η=时取()0f N =、()'

00f

=和

()0f β''=.这里β是一个常数(估计值)，这样就能根据式 (11.3.13) - (11.3.14), 求得()0f '''.接着，令η具有某一增量，分别计算f 、'f 、f ''和f '''，继续进行下去一直到边

界外边缘(在平壁的情况下它相当于 5.0η=). 此时，应该得到()'

1f ∞=. 如果这个要求

不能满足，就要修改原来的估计值，直至使()'

1f ∞=这个边界条件得到满足为止.，

表11-2 给出了数值计算方法。

表11-2 布拉修斯方程式的数值计算方法

在程序中为了能够方便的变换ff ''系数，我们把方程中ff ''系数设为epsilo ，从而把方程改写为：

epsilo**=0f f f '''''+ (11.3.13)

[Matlab 命令]：

?itBlasius(0,0,0,0.5,0.00001,0.1,5,0.001)

beta=0.435860

得出每一步的迭代结果： [Matlab 命令]：

?blasius(0,0,0.435860,0.1,5,0.001,[1,1,1,1,1]) beta=0.435860

表11-2给出了当0N =情况在区间[0，5]的部分计算结果.

Matlab 程序

function itaconv=Blasius(n,m,beta,dita,maxita,tol,options)

format long;

%%%%(f"'+epsilo*f*f"=0)%%%%%

%书上是epsilo=0.5

%epsilo=1;

epsilo=0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if options(5)==1

fprintf('\nbeta=%f\n\t\tita\t\tf\t\tf''\t\tf''''\t\tf''''''\n\n',beta); end

found=0;

x0=n;

y0=m;

z0=beta;

w0=-x0*z0*epsilo;

for ita=0 : dita : maxita

switch options(1)

case 1

multiplier=0.1;%推荐

otherwise

fprintf('\n出错了！\n\n');

break

end

switch options(2)

case 1

x=x0+multiplier*y0;%推荐

otherwise

fprintf('\n出错了！\n\n');

break

end

switch options(3)

case 1

y=y0+multiplier*z0;%推荐

otherwise

fprintf('\n出错了！\n\n');

break

end

switch options(4)

case 1

z=z0+multiplier*w0;%推荐

otherwise

fprintf('\n出错了！\n\n');

break

end

w=-x*z*epsilo;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if options(5)==1

fprintf('\t%2.5f\t%2.5f\t%2.5f\t%2.5f\t%2.5f\n',ita,x0,y0,z0,w0);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if abs(1-y)

found=1;

itaconv=ita;

break

end

x0=x;

y0=y;

z0=z;

w0=w;

end

if ~found

itaconv=-1;

end

function itBlasius(n,m,minbeta,maxbeta,dbeta,dita,maxita,tol)

found=0;

fprintf('\n白拉修斯方程数值计算迭代程序\n\n正在计算...\n\n');

for beta=minbeta:dbeta:maxbeta

ita=Blasius(n,m,beta,dita,maxita,tol,[1,1,1,1,0]);

if ita>-1

fprintf('\n满足白拉修斯收敛条件的beta已经找到，beta=%f\n\n',beta); found=found+1;

break

%若需找到所有beta, 请注释掉上一行的break

end

end

if found==0

fprintf('\n在给定的范围 (%f..%f) 内无法找到满足白拉修斯收敛条件的beta.\n\n',minbeta,maxbeta);

end

下面表11-3 给出了当0

N=情况在区间[0，5]的部分计算结果

表11-3：当0

N=部分计算结果.

小结

本章研究了常微分方程边值问题的数值求解方法. 主要介绍了三种边界条件下的定解问题, 通常主要有两大类求解方法, 即:打靶法算法和有限差分方法.

打靶算法对于线性或非线性边界值问题均适应. 边值问题的解可以视为一个 “弹道曲线”, 打靶算法是通过引进初值参数以近似原边界值问题, 即将原微分方程先看成一个初值问题求解, 需要利用一系列初值形式问题的解的序列来逼近, 沿着初值问题解的曲线积分，使该初值问题的解满足另一端点的边界条件, 每一步计算都要进行判断, 通过迭代得到新的初值参数, 以达到满足右端边界条件. 对于不同的边界条件, 需要引入不同形式的初值参数及建立不同代格式.

有限差分方法通常仅适应于求解常微分方程边值问题, 对于求解常微分方程初值问题有时则较困难. 通常通过将求解区间n 等份, 并分别用一阶和二阶中心差分来代替一阶导数和二阶导数, 略去)(2

h O 项, 结合定解条件, 得到相应于原微分方程的一个近似的线性差分方程组. 最后, 通过解线性方程组来得到原微分方程组的解.

本书给出了一个求解流体力学中动量边界层方程的例, 读者可根据实际问题进行相关练习.

习题11

一、 利用差分方法求解下列边值问题

(1)2(1)1, 11(1) (1)0, 0.5y x y x y y h ''?=+--<

(2)2(1), 11(1) (1)1, 0.5

y x y x y y h ''?=+-<

(0), (1)y q x y f x x y y αβ''-=<

采用如下差分方程

11111-1211

1

(2)(10)(10)

12

12

i i i i i i i i

i i

i i i L y y y y q y q y q y f f f h +-+++=

-+-++=++ βα==n 0y ,y , 试证明截断误差的阶为4()O h .

三、已知初值：

1

(0)(2)0

y xy y y y '''-+=??

==? （1） 将此边值问题化为两个初值问题.

（2） 取步长0.5h =写出用差分方法解此边值问题时的差分方程.

一、（1）1

012

2

0.2538,0.3369, 0.2538y

y y -===

（2）1230.7024,0.6244, 0.7024y y y === 二、略。 三、

（1）设()y x 1及()y x 2分别为以下两个初值问题的解

1(0)0,(0)1y xy y y y '''-+=??

'==? 及0

(0)0,(0)1

y xy y y y '''-+=??'==? 则11220(2)

()()()(2)

y y x y x y x y -=+

为边值问题的解。

（2）原微分方程边值问题化为如下三角方程

12314702101462011142y y y --????????????--=-???????????

?--??????

常微分方程边值问题与不动点定论文

目录 引言 (1) 1预备知识 (2) 定义1.1（奇异Sturm-Liouville边值问题的正解） (2) 引理1.1.1 (2) 定义1.2（凸集的概念) (3) 定义1.3锥的定义 (3) 定义1.4（全连续算子的概念） (3) 1.5 （常微分边值问题的定义） (4) 定义1.6混合单调算子得定义） (4) 2 常微分方程边值问题正解得存在性 (5) 2.1 奇异Sturm-Liouville常微分边值问题的正解存在在 (5) 子 (8) 2.2 一类二阶边值问题的存在性 (9) 3一类混合单调算子应用 (11) 3.1一类混合单调算子的存在唯一性?........................ 错误!未定义书签。 3.2 求常微分边值问题的例题 (13) 结束语 (15) 参考文献 (15) 致 (16)

常微分方程边值问题与不动点定 （数学与统计学院 11级数学与应用数学2班）指导教师：攀峰 引言 从历史上看在有了微积分这个概念以后，紧接着出现了常微分方程。发展初期是属于“求通解”得时代，当人们从初期的热潮中结束要从维尔证明了卡帝方程中是一定不会存在一般性的初等解的时候开始的，并且柯西紧接着又提出了初值问题，常微分方程开始从重视“求通解”转向重视“求定解”的历史时代。 大学我们都学习了常微分方程这门学科，如果要研究它的定解问题，我们首先就会知道是常微分方程的初值问题。然而，在科学技术、生产实际问题中，我们还是提出了另一类定解问题-边值问题。对于常微分方程边值问题，伟大的科学家最早在解决二阶线性微分方程时，提出了分离变量法。[]1.在牛顿时期，科学家们已经提出过常微分的边值问题，牛顿也对常微分边值问题进行过研究，并且在1666年10月牛顿已经在这个领域取得了很大的成就，但是由于种种原因当时并没有整理成论文，所以没有及时出版。但在1687年他终于把在常微分方程上研究的成果发表了，虽然不是在数学著作中，却是他的一本力学著作中（《自然哲学的数学原理》）。 在微积分刚创立时期，雅克.伯努利来自瑞士的科学家提出了远著文明的问题-悬链线问题，紧着的地二年著名数学家莱布尼兹就给出了正确的解答，通过对绳子上个点受力分析，建立了以下方程 这个方程满足的定解条件是y（a）=α；y(b)=β.这是一个典型的常微分方程的边值问题。从这开始，常微分边值问题已经是科学家研究微分方程是不可或缺的工具，我就简单列举几个例子：（比如种族的生态系统；梁的非线性震动）等。对于怎么研究它，

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

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

一阶常微分方程解法总结

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

常微分方程组(边值)

常微分方程组边值问题解法 打靶法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))

偏微分方程边值问题的数值解法论文

求解偏微分方程的边值问题 本实验学习使用MATLAB 的图形用户命令pdetool 来求解偏微分方程的边值问题。这个工具是用有限元方法来求解的，而且采用三角元。我们用个例题来说明它的用法。 一、MATLAB 支持的偏微分方程类型 考虑平面有界区域D 上的二阶椭圆型PDE 边值问题： ()c u u f α-??+=g (1.1) 其中 (1) , (2) a,f D c x y ?????=? ????? 是上的已知函数(3)是标量或22的函数方阵 未知函数为(,) (,)u x y x y D ∈。它的边界条件分为三类： （1）Direchlet 条件： hu f = (1.2) （2）Neumann 条件： ()n c u qu g ?+=g (1.3) （3）混合边界条件：在边界D ?上部分为Direchlet 条件，另外部分为Neumann 条件。 其中,,,,h r q g c 是定义在边界D ?的已知函数，另外c 也可以是一个2*2的函数矩阵，n 是沿边界的外法线的单位向量。 在使用pdetool 时要向它提供这些已知参数。 二、例题 例题1 用pdetool 求解 22D 1 D: 10u x y u ??-?=+≤??=?? (1.4)

解：首先在MATLAB 的工作命令行中键入pdetool ，按回牟键确定，于是出现PDE Toolbox 窗口，选Genenic Scalar模式． ( l ）画区域圆 单击椭圆工具按钮，大致在（0，0）位置单击鼠标右键，拖拉鼠标到适当位置松开。为了保证所绘制的圆是标准的单位园，在所绘园上双击，打开 Object Dialog 对话框，精确地输入

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

第十一章 常微分方程边值问题的数值解法 工程技术与科学实验中提出的大量问题是常微分方程边值问题．本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法. 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 ≤≤上式如下；

常微分方程数值解法

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 。

常微分方程组(边值)

常微分方程组边值问题解法 打靶法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

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

常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山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个系数

两点边值问题的两种数值解法

常微分方程组两点边值问题的数值解法 ----张亚苗2011年9月 3)1(1)0(04===-''y y y y 可化为微分方程组3 )1(1)0(41221==='='y y y y y y 方法一：配置法 Matlab 程序： function bvcollation clc solinit = bvpinit(linspace(0,1,20),[100 600]);% sol = bvp4c(@twoode,@twobc,solinit); x = linspace(0,1,20); y = deval(sol,x); y' plot(x,y(1,:),x,y(2,:)); end %微分方程组 function dydx = twoode(x,y) dydx = [ y(2) 4*y(1)]; end %边值条件 function res = twobc(ya,yb) res = [ ya(1)-1 yb(1)-3]; end 运行结果： 1.0000 -0.4203 0.9834 -0.2117 0.9777 -0.0055 0.9828 0.2007 0.9988 0.4091 1.0259 0.6220 1.0644 0.8419 1.1147 1.0710 1.1774 1.3121 1.2531 1.5677 1.3427 1.8407 1.4472 2.1341 1.5678 2.4512 1.7057 2.7954 1.8626 3.1707 2.0401 3.5811 2.2402 4.0313 2.4652 4.5261 2.7175 5.0712 3.0000 5.6724

常微分方程数值解法

第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 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 阶精度。

常微分方程数值解法

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

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

二阶常微分方程边值问题的数值解法 摘要 求解微分方程数值解的方法是多种多样的，它本身已形成一个独立的研究方向，其要点是对微分方程定解问题进行离散化．本文以研究二阶常微分方程边值问题的数值解法为目标，综合所学相关知识和二阶常微分方程的相关理论，通过对此类方程的数值解法的研究，系统的复习并进一步加深对二阶常微分方成的数值解法的理解，为下一步更加深入的学习和研究奠定基础． 对于二阶常微分方程的边值问题，我们总结了两种常用的数值方法：打靶法和有限差分法．在本文中我们主要探讨关于有限差分法的数值解法．构造差分格式主要有两种途径：基于数值积分的构造方法和基于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.

常微分方程的数值解

实验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