常微分方程数值解法在动力天文中的应用.docx
- 文档编号:28400859
- 上传时间:2023-07-13
- 格式:DOCX
- 页数:23
- 大小:225.33KB
常微分方程数值解法在动力天文中的应用.docx
《常微分方程数值解法在动力天文中的应用.docx》由会员分享,可在线阅读,更多相关《常微分方程数值解法在动力天文中的应用.docx(23页珍藏版)》请在冰豆网上搜索。
常微分方程数值解法在动力天文中的应用
青岛科技大学
学士学位论文
(2010年〜2014年)
题目:
常微分方程数值解法在动力天文中的应用
学院:
专业班级:
学生姓名:
数理学院
信计101
潘力杨
学号:
1006010129
指导教师:
王斌
起讫日期:
2014年3月
—2014年6月
常微分方程数值解法在动力天文中的应用
、尸■、亠
前言常微分方程,作为数学分支的新发展在人类认识自然,改造自然方面发挥了巨大的力量。
牛顿在研究天体力学和机械动力学时,一微分方程为工具,从而在理论上得出了行星运动的规律。
后来,法国天文学家勒维烈和英国天文学家亚当斯都各自分别利用了常微分方程,演算出了当时为止的海王星位置及其运动规律。
如今,常微分方程更是在许多学科领域内发挥着重要作用。
小行星、人造天体的运动及其规律,导弹轨道的计算,飞机空气动力学的研究,化学反应过程问题等,这些问题都可以转化为求常微分方程的解,或研究解的性质的问题。
虽然常微分方程存在众多的研究和应用领域,但事实上,只有极少数诸如线性系数常微分方程和个别典型微分方程的解能用初等函数,特殊函数或它们的级数与积分表达。
而在变系数常微分方程的解析求解上已经比较困难了,更别说一般的非线性常微分方程了。
大多数情况下,常微分方程只能用近似方法求解,近似方法可大致分为两大类:
用级数解法,逐次逼近法等求展开式的近似解析法;和在方程一些离散点上给出近似解得数值解法。
在动力天文的具体应用中,事实上由于所处的角度关系,我们很难宏观直接地去观测天体的具体运动状态,并直接以经典物理的力学方法去计算它的准确运动轨道及规律。
由于涉及到太多的不确定因素以及运动状态的复杂性,我们只能尽量地以我们的实际观测值为参数,建立模型,运用数学的方法去研究,拟合它的运动轨迹,以此为基础才能进一步作出定性的分析。
在天体运动方程的求解问题上,基本上来说也只有二体问题能给出严格解。
因此为了得出满足具体观测精度的定量解,我们不得不求助于数值解法。
利用数值解法求解常微分方程初值问题,一般只需要在满足规定精度的条件下求得解在若干相应点上的近似解,或是求出解的近似表达式就行,而无需使用复杂的方法试图去拟合解的表达式,或花费大量的计算去得出形式解。
由于它的理论简单,可操作性强且能适应各种复杂的问题条件,因此在数学建模领域有着独一无二的重要地位。
在过去技术尚未得到发展的时代,在运用数值解法求解天体运动方程的过程中,处于计算精度,方程稳定性和收敛性的需要,我们往往得做出大量繁杂重复的运算从而必须花费大量的时间,可靠性也得不到保证,这也反过来制约了数值解法的应用。
但现如今,随着计算机的普及以及编程运算技术的成熟运用,关于数值解法的应用领域也出现了新的契机,设定步长,循环次数的程式化运算大大缩短了我们的运算时间同时也提高了运算的精度,例如只要通过Matlab就可以完成简单方程求解问题的建模、迭代、作图、误差分析等一系列的工作,大大减轻了我们的负担同时也提高了运算的效率。
使得我们可以更从容地将数值解法运用到天体运动求解问题的更深邻域中去。
第一章常微分方程数值解法介绍
常微分方程可分为线性、非线性、高阶方程与方程组等类,其中线性方程包含于非线性类中,高阶方
程可化为一阶方程组,若方程组中的所有未知量看作一个向量,则方程组可写成向量形式的单个方程。
因
此研究一阶常微分方程的初值问题
dy
f(x,y),^x
dx(9-1)
y(a)=y。
的数值解法具有典型性,其中方程的解为y:
R>R。
只有保证问题(9-1)的解存在惟一的前提下,其数值解法的研究才有意义。
由常微分方程的基本理论,我们有:
定理9.1如果(9-1)中的f(x,y)满足条件
(1)f(x,y)在区域D={(x,y)a兰x兰b,-旳vy£+旳}上连续;
(2)f(x,y)在D上关于y满足Lipschitz条件,即存在常数L,使得
f(x,y)—f(x,y)兰Ly—y|
则初值问题(9-1)在区间[a,b]上存在惟一的连续解y=y(x)。
所谓数值解法,是通过常微分方程离散化而给出解在某些节点上的近似值。
对于问题(9-1),在区间〔a,bl引入若干节点
a=x0:
为:
X2:
HI:
Xn=b
m二xn1-Xn,(n=0,1川l,N-1)称为由xn到Xn1的步长,当hn=h(常数)时称为定步长,否则称为变
b—a
步长。
多数情况下,采用等步长,即h,xn=a,nh(n=1,2,|I(,N)。
记问题(9-1)的精确解为y(x),
N
记y(Xn)的近似值为yn,记f(xn,yn)为fn,『ng=1,2/,N)称为问题(9-1)的数值解。
求初值问题数值解的方法一般采用步进法,即在计算出y,^n后计算yn,,方法分为单步法和多步法。
单步法是指在计算yn1时只利用yn,而多步法在计算yn1时不仅要利用y.,还要利用前面已经计算出来的若干个yn_j,j=1,2,川,I-1°我们称要用到Yn和yn_j,j=1,2,|11,I-1的多步方法为I步方法。
不论单步法和多步法,它们都有显式方法和隐式方法之分。
显式单步法的计算公式可以写为:
yn+1=ynh(xn,yn,h),此公式右端不含ynd。
隐式单步法的计算公式可以写为:
yn+1=yn•h'(Xn,yn,yn.1,h),此公式右端含有yn-1,从而是『n计的
方程式,需要求解方程或者采用迭代法。
显式多步法的计算公式可以写为yn+1=%•h(Xn,yn,yn-1,川,丫心・1,h)。
隐式多步法的计算公式可以写为yn+1=ynh(Xn,yn!
yn,yn-1,,yn41,h),k-1-1。
显式公式与隐式公式各有特点。
显式公式的优点是使用方便,计算简单,效率高,其缺点是计算精度
低,稳定性差。
隐式公式正好与它相反,它具有计算精度高,稳定性好等优点,但求解过程很复杂,一般采用迭代法。
为了结合各自的优点,通常将显式公式与隐式公式配合使用,由显式公式提供迭代初值,再经隐式公式迭代校正。
1.1欧拉方法(显)
欧拉方法的建立可以通过以下三种方法:
(一)用差商近似导数,用向前差商y(XnQ「y(Xn)代替y(xn),则得
(n=0,1,L,N-1)
(9-2)
这样,问题(9-1)的近似解可以表示为
Yn1=Ynhf(Xn,Yn)
Y0=Y(Xo)
按式此式由初值y0经过N步迭代,可逐次算出y1,y2/yN,此方程称为差分方程。
式(9-2)称为求解
一阶常微分方程初值问题(9-1)的欧拉公式,也称显式欧拉公式。
需要说明的是,用不同的差商近似导数,将得到不同的计算公式。
(2)用Taylor多项式近似,把y(Xn』在Xn点处Taylor展开,取一次多项式近似,则得:
h2„y(Xn1)=y(Xn)hy(Xn)万y()
h2
二y(Xn)hf(Xn,y(Xn))—y()[Xn,Xn1]
设步长h的值较小,略去余项,并以yn代替y(xn),便得:
%1二ynhf(Xn,yn)
(3)将问题(9-1)中的微分方程在区间[X^Xn」上两边积分,可得:
Xn+
y(Xn卅)—y(Xn)=Jf(x,y(x))dx(n=0,1,…,N-1)
Xn
用yn1,yn分别代替y(Xn1),y(xn),若对右端积分采用取左端点的矩形公式,即:
Xn1
Jf(x,y(x))d^hf(Xn,yn)
xn
同样可得到显式欧拉公式。
1.2欧拉方法(隐)
在微分方程离散化时,用向后差商代替导数,即y(xn1)也山心°,则得到如下差分方程
hyn1fn■hf(Xn1,yn1)(n=0,1,|11,N-1)
«(9-3)
y。
二Y(x。
)
此公式称为求解问题(9-1)所得的数值解称为隐式欧拉法。
隐式欧拉法与欧拉公式(9-2)式上相似,但实际计算时却复杂得多。
欧拉公式(9-2)计算yn★的公式中
不含yn1,但是隐式欧拉法计算yn1的公式中含有yn1。
在求解yn1时,yn,Xn-1为已知,『n1是方程
yn^ynhf(Xn!
yn1)的根。
一般说来,这是一个非线性方程,当不能精确求解得到yn1的表达式时,
我们需要运用简单迭代法来求解。
迭代格式为:
■'[0]
yn1ynhf(Xn,yn)ynv'ynhf(Xni,ynk]i)(k=0,1,2,川)
由于f(x,y)满足Lipschitz条件,所以
yn羊-yn+=hf(Xn^yn』—f(X.申,y.』EhLyn申—yn申
由此可知,只要0:
:
:
hL:
:
:
1,迭代法就收敛到解yn勺。
1.3梯形公式
利用数值积分方法将微分方程离散化时,若用梯形公式计算式中右端积分,即
冷ih
ff(x,y(x))dx止二[f(Xn,y(Xn))+f(Xn+,y(Xn^))]
Xn2
并用yn,yn1代替y(Xn),y(Xni),则得计算公式
这就是求解初值问题(9-1)的梯形公式。
梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为
yn°fr=yn+hf(Xn,yn)扌h-
[yn.=yn+2」(Xn,yn)+f(Xn和ynk;)](k=0,1川)
由于函数f(x,y)关于y满足Lipschitz条件,所以
[k出]
h
yn卅_yn+
=—
2
[k]
f(XnH1,yn^)-f(Xn卄y^)兰牛|yn^-Yn卅
2hL
其中L为Lipschitz常数。
因此,当0§:
:
:
1时,迭代法就收敛到解yn计。
1・4改进的欧拉法
相比于显式欧拉法和隐式欧拉法,梯形方法提高了精度,但是梯形方法,在应用迭代公式进行实际计算时,每迭代一次,都要重新计算函数f(x,y)的值,而迭代又要反复进行若干次,计算量很大。
当函数
f(x,y)比较复杂时,这个问题会变得更加突出。
为了控制计算量,我们先用欧拉公式求得一个初步的近似值yn1,称之为预测值,预测值yn1的精度可能很差,再用梯形公式将它校正一次得yn计,称为校正值。
这样的预测校正系统通常称为改进的欧拉公式。
即
预测:
ynynhf(Xn,yn)
校正:
yn1二yn■h〔f(Xn,yn)■f(Xn1,Vn1)1
2
为了便于编制程序上机,将上式改写成
yp=yn+hf(Xn,yn)
*yq=yn+hf(Xn+h,『p)(9-5)
1/、yn#=2(yp+yq)
算法实现如表
算法9.5
1输入a,b,f(x,y),h,初值y0b-a2N,n=0,x二a,y二y0h
3计算
y^yhf(x,y)y^yhf(xh,yp)
1-(ypyq)二y,xh-x输出(x,y)
4若n:
:
:
N-1,置n•1=•n,转③;否则停机。
1.5龙格—库塔(Runge-Kutta)法
Runge-Kutta法由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明。
设初值问题(9-1)的解y=y(x)・C1[a,b],由微分中值定理,我们知道,必存在[x,,xnd],使
yg)=y(xj+hy©)^Ynhf(,y())=Y^hK*
其中K*叫做y(x)在[xn,Xn1]上的平均斜率。
对于平均斜率K*,只要提供一种计算方法,该公式就给出一种数值解公式。
例如,用K1二f(Xn,yn)计算K*,就得到一阶精度的欧拉公式;用K2二f(xn1,yt).
替代K*,就得到隐式欧拉公式;如果用K1,K2的算术平均值计算K*,则可得到二阶精度的梯形公式。
由
此可以设想,如果在[x,,x,1]上能多预测几个点的斜率值,用它们的加权平均值来计算K*,就有望得到具
有较高精度的数值公式,这就是Runge-Kutta法的基本思想。
Runge-Kutta公式的一般形式是
s
Kj=f(人+c,hy瓦ha)K
j1
\s(i=1川S)(9-6)
Yn厂ynh'bKi
Li#
其中h为步长,s称为方法的级数,Ki是y二y(x)在xn^h(0乞q^1)点的斜率预测值,q,t|,aij均
为常数。
这些常数的选取原则是使公式(9-13)具有尽可能高的精度。
公式(9-13)叫做s级Runge-Kutta公式,
简称RK方法。
Runge-Kutta公式根据系数aj选择的不同,可以分为显式方法和隐式方法。
当aj=0,H;j时,方法为
显式方法,否则为隐式方法。
显式Runge-Kutta方法可用下面的形式表示
Ki二f(XnCih,ynh'ajKj)
s
yn1=ynh'biKi
i4
(i胡川,s)
(9-7)
2级(S=2)显式RK方法的公式为
K^f(Xn,yn)
“心=fX+C2h,yn+1^2!
®
=yn+h(bK+匕2心)
这里bi,b2,C2,a2i为待定常数。
(9-8)
常见而二阶公式,俗称中点公式为:
Ki=f(Xn,yn)
«K2=f(xn+0.5h,yn+0.5hKJ
必十=yn+h©
(9-9)
类似地,对s=3和S=4的显式RK公式,通过更复杂的计算,可以导出三阶和四阶
用的三阶和四阶RK公式分别为
K=f(人,%)
hh
K2=f(x^-,y^-Ki)
«22
K3=f(Xn+h,yn—hKi+2hK2)
yn十yn』(Ki+4K2+K3)
L6
和
K=f(Xn,yn)
hh
K2=f(Xn右,yn+-Ki)
hh
卞3二f(Xn右小+^心)
K4=f(Xn+h,yn+h&)
h
yn+=yn+—(心+2心+2心+心)
I.6
RK公式,其中常
(9-10)
(9-11)
式(9-22)称为经典的四阶RK公式,通常说四阶RK方法就是指用式(9-22)求解,具体算法实现如表9.3。
表9.3
算法9.2
②输入a,b,f(x,y),h,初值y02N=-—,n=O,x0=ah
3计算
心=f(xo,yo)hh
K2=f(xo■二,yo■二Ki)
22hhK3=f(xo■—,yo■—心)22K4=f(X。
h,yohK3)y二y。
-(Ki2K22K3K4)6X。
h=x,y=yo输出(x,y。
)4若n:
:
N-1,则置n•1=n,转③;否则停机。
1.6线性多步法
之前的各种数值解法,都是单步法,即在每一步计算时,只要前面一个值yn已知的条件下就可以计算
出yn1。
单步法的特点是,可以自成系统进行直接计算,因为初始条件只有一个已知y0,由y0可以计算y1,
不必借助于其它方法,这种单步法是自开始的。
如果考虑在计算值yn卅时,能够比较充分地利用前面的已知信息,如yn和yn斗,j=1,2,|||,1-1,那么
就可望使所得到的yn1更加精确。
这就是多步法的基本思想。
我们称要用到yn和yn_j,j=1,2,川,1-1的多步法为I步方法。
多步法中最常用的是线性多步法,其一般公式是
rr
yn1=7:
iyn「fn4(9-12)
i=0im
其中:
k,\均为常数,=y(XnJ,fk=f(Xk,yQ(k=n1JH,n-r)。
当:
4=0时,上式称为显式,否则称为隐式。
定义:
设y(x)是初值问题(9-1)的精确解,线性多步法(9-12)在Xn1上的局部截断误差定义为:
rr
Tn1二y(Xn1)'h':
kYn*
k=0k7
若Tn1=o(hp1),则称方法(9-12)是p阶方法。
线性多步法的构造原理:
利用Taylor展开导出线性多步公式(9-12)的基本方法是将线性多步公式(9-12)
在Xn处进行Taylor展开,然后与y(Xn1)在Xn处的Taylor展开式相比较,要求它们前面的项重合,由此确定参数:
k,'k。
设初值问题(9-1)的解y(x)充分光滑,记丫讣)曲(门(人)仃=0,1川),把它在焉处展开成
Taylor级数,则有:
(X-Xn)j
及
y(X」j^ynj)rSP。
川
j4(j-1)!
p!
用Xn±=x0•(n-k)h替代式以上两式中的x,则得:
yn「y(Xn^yn丁畔ynj)代SP1)川
j二j!
(P+1)!
把这两个式子代入式(9-12),得:
(j=1川,p)
广(-k)2k=1
k7
可见,只要式(9-12)的系数:
k,\满足上式,方法(9-12)就具有p阶精度。
由此可得局部截断误差为:
P1rr
T(1=厂1」(-k)p1:
k-(p1)'(-k)pry(p1)O(hp2)
(p+1)!
-心k=1」
式(9-13)共有p1个方程,2r3个待定常数:
rFjll,-,、「。
,^」l「r,只要2r•3—pT,
即r-子一1,就可以由式(9-13),定出具有p阶精度的公式(9-12)的系数:
\,:
k.
因为\=0,式(9-15)称为Adams隐式公式,它是四阶公式,局部截断误差为
Tn1「焉h5yn5)O(h6)
yn,yn4,川,yn_r之后,才能使用。
所以开头的%,丫2,川,齐
720
由于线性多步法公式只有在知道了前面的
的值必须先用同阶的单步法(例如Runge-Kutta法)求出,然后才能用线性多步法公式。
第二章常微分方程数值解法的程序实现
对于数值求解微分方程(组),MATLAB有专门的求解器可以调用。
MATLAB中求微分方程的数值解命令
的调用格式为:
[X,Y]=solver(odefun,tspan,y0)
其中solver为MATLAB^求解微分方程的求解器:
例如ode45、ode23、ode113、ode15s、ode23s、ode23t、
、dy
ode23tb等;odefun是要求解的显式常微分方程:
」=f(x,y),;tspa门=区兀血]为积分区间;y0为初
dx
始条件。
输出的X为tspanrXo’Xend]区间上点的列向量,丫为对应于X上各个点的数值解所组成的向量。
要
获得问题在其他指定时间点Xo,Xi,X2,|,Xend上的解,则令tSpa门=恢,花,X?
」]I,Xend](要求是单调的)。
因为没有一种算法可以有效地解决所有的求解常微分方程问题,为此,MATLAB提供了多种求解器
Solver,对于不同的常微分方程问题,采用不同的Solver,一些常用的Solver总结如下:
表9.8
求解器
Solver
求解方程类型
Solver特点
说明
ode45
非刚性
单步算法;4、5阶Runge-Kutta方程;
3
累计截断误差达(Ax)
大部分场合的首选算法
ode23
非刚性
单步算法;2、3阶Runge-Kutta方程;
3
累计截断误差达(占X)
使用于精度较低的情形
ode113
非刚性
多步法;Adams算法;高低精度均
可到103〜105
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法;Gear's反向数值微分;精度中等
若ode45失效时,可尝试使用
ode23s
刚性
单步法;2阶Rosebrock算法;低
精度
当精度较低时,计算时间
比ode15s短
ode23tb
刚性
梯形算法;低精度
当精度较低时,计算时间
比ode15s短
要特别的是:
ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题
的解的MATLAB的常用程序,其中:
ode23采用Runge-Kutta2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度;
ode45贝慄用Runge-Kutta4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度。
3.1欧拉法的Matlab程序实现
1
Plot(x,y,x1,y1)
Functiony=tier(dyfun,x,y,h)yO=y;e=1e_a;K=1e+4;y=y+h*feval(dyfun,x,y);y1=y+2*e;k=1;
Whileabs(y-y1)>e
y仁y;
向前欧拉法
Function[x,y]=naeuler(dyfun,xspan,y0,h)X=xspan
(1):
h:
xspan
(2);
y
(1)=yo;
Forn=1:
length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));
End
x=x';y=y'
2向后欧拉法(隐1=(1+2*心0.5;
Function[x,y]=naeulerb(dyfun,xspan,y0,h)x=xspan
(1):
h:
xspan
(2);
y(i)=y(0);
Forn=1:
length(x)-1
y(n+1)=iter(dyfun,x(n+1),y(n),h);
隐式方法的程序去掉,整篇论文只给显式方法的程序
3改进欧拉法
Function[x,y]=naeuler2(dyfun,xspan,y0,h)x=xspan
(1):
h:
xspan
(2);
y
(1)=y(0);
Forn=1:
length(x)-1
k1=feval(dyfun,x(n),y(n));y(n+1)=y(n)+h*k1;
k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;
End
x=x';y=y';
x1=0:
0.2:
1;y仁(1+2*x1)90.5;
Plot(x,y,x1,y1)
2.2龙格库塔法的Matlab程序实现
三阶龙格—库塔公式:
functiony=DELGKT3_kuta(f,h,a,b,y0,varvec)formatlong;
N=(b-a)/h;
y=zeros(N+1,1);
y
(1)=y0;
x=a:
h:
b;
var=findsym(f);
fori=2:
N+1
K1=Funval(f,varvec,[x(i-1)y(i-1)]);
K2=Funval(f,varvec,[x(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法 动力 天文 中的 应用