常微分方程地数值解法及其应用.docx
- 文档编号:11390514
- 上传时间:2023-02-28
- 格式:DOCX
- 页数:65
- 大小:303.39KB
常微分方程地数值解法及其应用.docx
《常微分方程地数值解法及其应用.docx》由会员分享,可在线阅读,更多相关《常微分方程地数值解法及其应用.docx(65页珍藏版)》请在冰豆网上搜索。
常微分方程地数值解法及其应用
引言
自然界中很多事物的运动规律可用微分方程来刻画。
常微分方程是研究自然科学和社会科学中的事物、物体和现象运动、演化和变化规律的最为基本的数学理论和方法。
物理、化学、生物、工程、航空航天、医学、经济和金融领域中的许多原理和规律都可以描述成适当的常微分方程,如牛顿的运动定律、万有引力定律、机械能守恒定律,能量守恒定律、人口发展规律、生态种群竞争、疾病传染、遗传基因变异、股票的涨幅趋势、利率的浮动、市场均衡价格的变化等,对这些规律的描述、认识和分析就归结为对相应的常微分方程描述的数学模型的研究。
因此,常微分方程的理论和
方法不仅广泛应用于自然科学,而且越来越多的应用于社会科学的各个领域。
它的学术价值是无价的,应用价值是立竿见影的。
求一阶常微分方程的解是数学工作者的一项基本的且重要的工作。
由于国内外众多数学家的努力,使此学科基本上形成了一套完美的学科体系;由于该问题比较复杂且涉及的面广,使得有些问题的解析解很难求出,而对于一些典型的微分方程(如线性方程、某些特殊的一阶非线性方程等)可以运用基本方法求出其解析解,并在理论上可以根据初值问题的条件把其中的任意常数完全确定下来。
然而,在生产实际和科学研究中所遇到的微分方程往往很复杂,在很多
情况下都不可能给出解的解析表达式,有时即使能求出形式的解,也往往因计算量太大而不实用,而且高次代数方程求根也并不容易,所以用求解析解的方法来计算微分方程的数值解往往是不适宜的。
实际上,对于解微分方程初值问题,一般只要求得到解在若干个点上的近似解或者解的便于计算的近似表达式(只要满足规定的精度就行)。
所以,研究数学建模中常微分方程模型理论性数值解法迫在眉睫。
本文研究的数值解法主要是针对常微分方程初值问题多种数值解法精度比较而言。
从而得到更常
用的数值解法在微分方程模型中的应用
在自然科学和经济的许多领域中。
常常会遇到一阶常微分方程的初值问题
黑f(X,y),a
y(x。
)y°.
这里f(x,y)是充分光滑,即关于x或y满足李普希茨条件的二元函数,yo是给定的初始值,y(xo)yo称为初始条件。
常微分方程初值问题的数值解法一般分为两大类:
单步法:
所谓单步法是指这类方法在计算yn1时,只用到前一步的值Xni,Xn,yn,然后逐步往下计算。
这个算法的代表是龙格库塔算法,简称R—K方法。
四阶显示
Runge——Kutta方法是求解普通常微分方程初值问题数值解法中的重要方法,而隐
式Runge——Kutta公式是求解刚性常微分方程初值问题的重要方法。
多步法:
这类方法在计算yn1时,除了用到前一步的值Xni,Xn,yn,之外,还要用到
Xnp,ynp(P1,2,,k;k0),
这前面k步的值,这个算法的代表就是阿达姆斯(Adams)方法。
一般用微分方程建立的动态模型来描述动态过程的变化规律,但对于某些实际问
题,建模的主要目的并不是要寻求动态过程每个瞬时的性态,而是研究某种意义下稳
定状态的特征,特别是当时间充分长以后动态过程的变化趋势,为了分析这种稳定与
不稳定的规律,常常不需要求解微分方程,而可以利用微分方程的稳定性理论,直接研究平衡状态的稳定性就行了。
微分方程的数值解法有显式解和隐式解法,一般来说,隐式解要优于显式解。
欧拉方法是一种最简单的单步法,计算量小,但精确度比较低。
一般的初值问题,多采用改进的欧拉方法,因为它的数值稳定性和计算精确度比一般的欧拉方法好。
龙格一
—库塔方法是一类应用较广的高精度单步法,当解充分光滑时的4阶龙格一一库塔方法一般可以达到很高的精确度。
常微分方程的初值对计算方法的收敛是有影响的。
为
了更好地比较这几种常用的方法,本文采用这几种数值方法对被积函数光滑连续,初
值精确的微分方程做了数值试验。
第一章微分方程的基本概念
自变量、未知函数均为实值的微分方程称为实值微分方程;未知函数取复值或自变量及未知函数均取复值时称为复值微分方程。
在本文中只讨论实值微分方程。
1.1微分方程和解
含有未知量的等式称为方程,它表达了未知量所必须满足的某些条件。
方程是根据对未知量所进行的运算来分类的,如代数方程、超越方程等。
微分方程与代数方程和超越方程不同,它的未知量是函数,对其所施加的运算涉及求导或微分。
1.1.1微分方程的概念
一般说来,微分方程就是联系自变量、未知函数以及未知函数的某些导数或微分的关系式。
如果其中未知函数是一元函数,则称为常微分方程;如果未知函数是多元函数,并且在方程中出现偏导数,则称为偏微分方程。
本论文主要介绍常微分方程,也简称微分方程或方程。
在一个微分方程中,出现的未知函数导数的最高阶数,称为
方程的阶。
以y为未知函数,x为自变量的一阶常微分方程的一般形式可表示为
F(x,y,yJ0,(1.1.1)
如果在(1.1.1)中能将y,解出,则得到方程
y,f(x,y),(1.1.2)
或
M(x,y)dxN(x,y)dy0,(1.1.3)
也称(1.1.1)为一阶隐式微分方程,(1.1.2)为一阶显式微分方程,(1.1.3)为一阶微分方程的微分形式。
n阶隐式方程的一般形式为
F(x,y,y,,,y(n))0,(1.1.4)
n阶显式方程的一般形式
(n)丄/,(n1)\/>>u\
yf(x,y,y,,y),(1.1.5)
方程(1.1.4)中,如果函数「对未知函数丁和它的各阶导数y,y;,y(n1)都是一次的,则称其为线性常微分方程,否则,称其为非线性微分方程。
以y为未知函数,x
为自变量的n阶线性微分方程具有如下形式:
y(n)Pi(x)y(n1)Pni(x)y‘Pn(x)yf(x),(1.1.6)
1.1.2微分方程的解一一通解与特解
定义1.1设函数y(x)在区间I上具有直到n阶的导数。
如果把y(x)代
入方程(1.1.4),有
F(x,(x),'(X),,(n)(x))0
在区间:
上关于x恒成立,则称y(x)为方程(1.1.4)在区间【上的一个解。
依据定义1.1可以直接验证:
:
2
(1)函数ysin(arcsinxC)是方程些f在区间(1,1)上的解,其中C是任
dxV1x2
意常数。
另外,该方程还有两个解y1(x(1,1)),它们不包含在前面解中。
、d2x、
(2)函数xC1costC2sint是方程一2xx0在区间(,)上的解,其中G
dt
和C2是独立的任意常数。
当然,xsint,xcost都是方程的解,它们包含在前面解中。
从上面的讨论中看到事实:
微分方程的解可以包含任意常数,其任意常数的个数可以多到与方程的阶数相等,也可以不含任意常数。
二阶常微分方程(1.1.4)的含有
n个独立的任意常数C1,C2,,Cn的解y(x,C1,C2,,5)称为该方程的通解,而方
程满足给定条件的解y(x)称为特解。
一般地,方程的特解可由其通解中任意常数
取确定的常数导出。
以隐函数形式表示的通解称为通积分,而以隐函数形式表示的特解称为特积分,对于通解或者通积分的说法或使用,通常是不加区分的。
另外方程的通解不一定表示方程的所有解。
为了便于研究方程解的性质,常常需要考虑解的图像,或者以图形方式表示微分方程的解。
一阶微分方程的特解y(x)的函数图像是xOy平面上的一条曲线,称为微分方程的积分曲线,而通解y(x,C)的函数图像是平面上的一族曲线,称为积分曲线族。
1.2常微分方程初值问题的一般提法
常微分方程初值问题的一般提法是求函数y(x),axb,满足
dy
f(x,y),axb,dx
y(a).
(1.2.1)
(1.2.2)
其中f(x,y)是已知函数,是已知值。
假设f(x,y)在区域D(x,y)axb,y上满足条件:
(1)f(x,y)在D上连续;
(2)f(x,y)在D上关于变量丁满足Lipschitz条件:
f(x,yi)f(x,y2)L%y2,axb,/y,(123)
其中常数L称为Lipschitz常数。
我们简称条件
(1)、
(2)的基本条件。
由常微分方程的基本理论,我们有:
定理1.1当f(x,y)在D上满足基本条件时,一阶常微分方程初值问题(1.2.1)、
(1.2.2)对任意给定存在唯一解y(x)在a,b上连续可微。
定义1.2方程(1.2.1)、(1.2.2)的解y(x)称为适定的,若存在常数0和
K0,对任意满足条件及||(x)的和(x),常微分方程初值问题
(1.2.4)
f(x,z)(x),axb,dx
z(a)a.
由常微分方程
存在唯一解z(x),且|y(x)z(x)|K||。
适定问题的解y(x)连续依赖于(1.2.1)右端的f(x,y)和初值
的基本理论,还有:
定理1.2当f(x,y)在D上满足基本条件时,微分方程(1.2.1)、1.2.2)的解y(x)
是适定的。
在本章中假设f(x,y)在D上满足基本条件,从而(1.2.1)、(1.2.2)的解y(x)存
在且适定。
般的一阶常微分方程组初值问题是求解
Yi(a)i,i1,2,,n.
(1.2.5)的向量形式是
dy
(1.2.6)
F(x,y),axb,dx
y(a).
其中
y(x)(y1,,yn(x))T,F(x,y)(f,x,y),,fn(x,y))T,(1,n)T,
记D(x,y1,,yn)axb,y,i1,2,,n。
类似于定理1.1和定理1.2,有:
定理1.3若映射F(x,y)满足条件
(1)F(x,y)在D上是从Rn1到Rn上的连续映射;
(2)F(x,y)在D上关于y满足Lipschits条件;
F(x,yi)F(x,y2)|Lyiy2〔ab,yi,y2任意。
则常微分方程组初值问题(1.2.5)存在的唯一的连续可微解y(x)而且解y(x)是适定的。
高阶常微分方程初值问题一般为
其中,an为给定值。
弓I进新的变量函数
则初值问题(1.2.7)化成了一阶常微分方程组初值问题
yi(a)ii,i1,2,,n.
通过求解(1.2.9)得到(1.2.7)的解y(x)y1(x)
1.3初值问题数值解基本概念
a
(1.2.1),(1.2.2)的准确解为y(x),记y(Xk)的近似值为Yk,记侬皿)为fk
求值冋题数值解的方法是逼近法,即在计算出y:
ik后计算yk1。
数值的方法有单步与多步法之分。
单步法在计算yki时只利用yk而多步法在计算yki时不仅要利用yk还有利用前面已算出的若干个ykj,j1,2,,l1。
我们称要用到
yk,yki,,ykii的多步法为】步法。
单步法可以看作多步法,但两者有很大差别。
L步法只能用于yk,k1的计算,yo,yi,,%i要用其它的方法计算;而且在稳定性上单步法比I1的多步法容易分析;此外单步法容易改变步长。
单步法和多步法又都有显式方法和稳式方法之分。
单步显式法的计算公式可写成
yk1ykh(Xk,yk,h),(1.3.1)
隐式单步法的计算公式可写成
yk1ykh(Xk,yk,yk1,h),(1.3.2)
在(1.3.2)中右端项显含yk1。
从而(1.3.2)是yk1的方程式,要通过解方程求出yk1。
显式多步法计算公式为
yk1ykh(Xk,yk,yk1,h),(1.3.3)
而隐式多步法计算公式为
yk1ykh(Xk,yk1,yk,,yki1,h),(1.3.4)
右端项含yk1。
多步法中一类常用方法是线性多步法
i1i1
yk1aiyk1hifki,kI1,(1.3.5)
i0k1
其中0,1,,i1,1,0,,i1是独立于咗和「的常数。
10时(1.3.5)是显式的,
10时是隐式的。
数值解法及方法构造、误差分析等内容。
一些概念和定义在后面的论述中逐步引入。
第二章常微分方程的数值解法
一般说来通过初等积分法求解的方程仅有很少的一些类型,绝大部分从实际问题
中提出的微分方程往往求不出其解析解。
而在实际问题中,对于复杂微分方程的求解,一般只需要得到解在若干个点上的近似值或者解的便于计算的近似表达式(只要满足
所规定的精度)即可。
本章将选讲关于微分方初值问题的基本数值解法及其matlab
程序设计。
2.1微分方程数值解法的一般步骤
一阶微分方程的初值问题为
(2.1)
Sf(X,y),
y(Xo)yo.
寻求微分方程初值问题(2.1)的数值解,就是求解函数yy(x)在一系列离散点Xi(i1,2,,n)(称为结点)上精确值y(xJ,y(X2),,y(Xn)的近似值y",,y。
在使用数值解法求解微分方程初值问题时,一般是按一下步骤:
(1)引入点列Xi,其中XiXi1hi,i1,2,,h,称为步长。
为了便于使用计算机
进行编程计算,一般取步长为定值,即hih,
XiXoih,i0,1,,(2.2)
(2)寻求数值解的方法,即寻求由y「计算出yi(i1,2,,n)的递推公式。
(3)利用
(2)中的格式逐步求出近似解y1,y2,,yn。
2.2Euler法
Euler法又称Euler折线法,它是解常微分方程的数值方法中最简单的一种方法。
Euler法的基本思想是:
在每一个小区间上,用一条切线来代替原函数曲线。
从整体上看就是用一条折线来代替解(或积分)曲线,并以此来求取一系列离散结点处函数近似值。
2.2.1Euler法的基本思想
对于一阶微分方程初值问题
y,f(x,y),(2.3)
y(x。
)y°.(24)
考虑在结点Xi的导数y,(Xi)用差商y(Xi1)y(xi)代替,则(2.3)可近似地写成
h
y(Xi1)y(Xi)hf(Xi,y(Xi)),i0,1,2,,(2.5)
从Xo出发,由初值(2.4),并利用(2.5)就可以得微分方程精确解的值y(Xi)的近似值yiyohf(xo,y°),再以yi作为y(xj的近似值代入(2.5)的右端,可得yg)的近似值yyihf(Xi,yi),继续做下去,一般可写成
yi1yihf(Xi,yi),i0,i,2,,(2.6)
这种方法称为解初值问题的Euler法。
公式(2.6)称为求解初值问题(2.3)和
(2.4)的Euler公式。
Euler法的几何意义如图3.1所示。
初值问题(2.3)和(2.4)的解曲线y(x)过点
P0(X0,y°),从P0点出发以f(X0,y°)为斜率作直线,与直线xXi交于点R(Xi,yi),显然yiy°hf(X0,y°),再从Pi点出发以f(Xi,yi)为斜率作直线段,与直线xx?
交于点P2(X2,y2),以此类推。
这样就得到解曲线的一条近似折线P0PE,所以,Euler
法又称Euler折线法
Euler折线法的递推结构:
(1)令xixoh,使用计算公式
yiyohf(Xo,yo)
计算yi;
(2)输出Xi和yi,并使XoXi,yoyi,重复以上过程。
2.2.2Euler法的局部截断误差进行估计。
在用Euler法求解一阶微分方程初值问题(2.3)和(2.4)时,实际上是利用Euler
公式yiyihf(Xi,yi)进行计算,计算出的近似值yii与精确值y(Xii)之间必然是
有误差的,所产生的误差为
Riiy(Xii)yii,(2.7)
将函数yy(x)在点Xi处进行Taylor展开,
yy(x)y‘(Xi)(xx)x)2■y^(xxj3
将xXi1,hXi1Xi代入上式,则
注意到
可得
观察(2.8),可以发现,当步长h趋于零时,Rii是h2的同阶无穷小,记为0(h2).
由此可知,我们用Euler折线法求解一阶常微分方程初值问题(2.3)和(2.4)时,在某一结点处的局部截断误差与选择的步长的平方成正比。
如果求解公式的局部截断误差是O(hp1),则称该求解公式是p阶的,或者具有p阶精度的,因此,Euler折线法是具有一阶精度的一种数值方法。
作为理论误差的应用,在用计算机编程进行数值解计算时,可以通过式(2.8)
和给定的误差要求初步确定步长h和选取的点个数n。
例如,如果给定的计算误差为
E,xa,b则可通过
初步确定出步长h,进行适当调整,使得等式nhba成立,并尽可能地取:
:
为正数,h为有限小数。
2.2.3Euler折线法的数值实例及matlab源程序
例2.1使用Euler折线法求解初值问题
dy
1y,
dx0x1,
y(o)1.
取步长为h0.10
解应用Euler折线法,在计算机中运用matlab软件,程序为:
clearelf
szy=[];
y=1;
szy=[szy;y];
forx=0.1:
0.1:
1
y=y+(1+y)*0.1;
szy=[szy;y];
end
szy
输出为:
szy=1.0000
1.2000
1.4200
1.6620
1.9282
2.2210
2.5431
2.8974
3.2872
3.7159
4.1875
用折线法算出的初值问题在点x0:
0.1:
1处的数值解(存放在数组中),而输入:
y=dsolve('Dy=1+y','y(0)=1')
可以算得此微分方程初值问题的精确解:
y=2*exp(x)-1
x0:
0.1:
1处的纵坐标(存
为了比较精确解和近似解的误差,编写程序算出精确解在点
放在数组jqy中)
jqy=[];
forx=0:
0.1:
1
y=-1+2.*exp(x);
jqy=[jqy;y];
end
jqy
现在将折线和曲线在同一坐标系中画出来,如图2.2所示.
x=0:
0.1:
1;
plot(x,szy,'o',x,jqy)
45
IL[I[L[[I
4-
3.5-
3—
2.5一
2」
1.5-
1[Ic[I[[l[[
00.10.20.30.40.50.60.70.80.9
图2.2折线法的数值解与函数真实值
在点x0:
0.1:
1处的精确解与近似解相除所得的差。
为了减少误差,
从输出容易观察到:
此方法在经过多步以后,其误差会积累起来.
一种方法是减少步长dx的值。
另外,在近似算法上也有很多改进办法。
224梯形法
可以注意到,在利用Euler法求解微分方程初值问题时计算简单、计算量小、程序简洁,但是,计算精度却很低。
为了提高计算精确度,同时又保持计算简单,程序简洁的优点,所以给出了梯形法。
其基本思想是:
在每个小区间上,仍然用一条直线来代替原函数曲线,但是这条直线较Euler折线法中的切线更近似于解函数曲线。
为寻求这条直线,考虑微分方程
y(x)f(x,y),
在区间Xi,Xi1上对方程两边关于x积分,则
于是有
(2.9)
xi1
y(XiJy(xjf(x,y)dx,
xi
现在,利用左矩形公式来计算(2.9)中的积分,有
X1
f(x,y)dxf(Xi,y(Xi))(Xi1Xi)hf(人,y(Xi)),
Xi
由此可得
y(Xi1)y(Xi)hf(Xi,y(Xi)).
如果用yi1和y分别代替y(XiJ和y(xj,则可看出这个迭代公式就是Euler公式。
如果改用梯形公式来计算(2.9)中的积分,即
x11h
f(x,y)dx-f(Xi,y(Xi))f(Xi1,y(Xi1))(Xi1Xj);f(人,y(Xj))f(Xi1,y(Xi1))
x22
代入(2.9)中,得到
h
y(Xii)y(Xi)二f(Xi,y(Xi))f(Xi1,y(Xi1))
2
在数值计算格式中,用yii和y分别代替y(Xi1)和y(Xi),则得到
h
y(Xi1)y-f(Xi,yi)f(Xii,yiJ,i0,1,2,(2.10)
2
从Xo出发,由初值y(Xo)yo,利用(2.10)可以依次得如宀,,yn。
公式(2.10)是由数值积分中的梯形公式得出来的,将它称为梯形公式。
和前面
的Euler公式相比。
梯形公式计算%1是也只用到了前一步的数值yi。
但是,如果y已
知,将yi代入(2.10)的右端,一般就不能直接得到yi1,而需要通过其他方法(如迭代法)来求解,称梯形公式为单步隐式公式。
单步隐式公式的一般形式为
yi1yih(xy,yi1,h),i0,1,,n1,
y。
其中(Xi,yi,~,h)为增函数。
梯形公式的增量函数为
1
(Xi,yi,~,h)-f(X,y)f(xh,~)。
2
下面不讨论梯形法的局部截断误差。
将y(x1)在Xi处用Taylor展开,得:
yihy'(Xi)》y”(Xi)h1O(h2)
22
因此,梯形法的误差为
Riiy(xii)yi
v"(x)23h213
y(Xi)y(Xi)h2f^h0(h)Vihy,(Xi)-yy,,(Xi)h-O(h)
1
-0(h)
2
由此可知,梯形公式的计算误差与步长的三次方称正比。
可见,梯形法具有二阶
精度。
2.2.5改进的Euler法
比较梯形公式(2.10)和Euler公式(2.6),梯形公式是隐式方程,且所耗费的计算量大于Euler公式求出
~iiyihf(Xi,yJ
用~i-近似准确值y(Xii),称~i-为预测值。
然后用梯形公式将它校正为较为准确的值yi-,称由此得到的y「为校正值,即
h〜
yiiyif(Xi,yi)f(Xi1,yi1).
2
这样建立起来的预测-校正系统通常称为改进的Euler公式,即
〜i1yhf(Xi,yi),
hf()f(〜、(2.11)
yi1y-f(^,yi)f(x1,yi1)•
2
为了便于编程计算,可将改进的Euler公式(2.11)改写为下列形式
〜pyhf(Xi,yi),
(2.12)
〜cyihf(Xi1,yp),
1
yi1尹卩〜c).
改进的Euler法的递推结构为
⑴令XiXoh,利用式(2.12)依次计算出Yp,y和yi;
⑵输出xi,yi,并使xoXi,yoyi,重复以上过程。
226改进的Euler法求解微分初值问题实例及matlab源程序
例2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法 及其 应用