数值分析-第9章-常微分方程初值问题数值解法_第1页
数值分析-第9章-常微分方程初值问题数值解法_第2页
数值分析-第9章-常微分方程初值问题数值解法_第3页
数值分析-第9章-常微分方程初值问题数值解法_第4页
数值分析-第9章-常微分方程初值问题数值解法_第5页
已阅读5页,还剩102页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

9.1引言9.2简单的数值方法与基本概念9.3龙格-库塔方法9.4单步法的收敛性与稳定性9.5线性多步法9.6方程组和高阶方程,Ch9常微分方程初值问题数值解法,9.1引言,在科学技术中常常需要求解常微分方程的定解问题.主要有初值问题与边值问题,本章只考虑初值问题。,常微分方程初值问题中最简单的例子是人口模型。设某特定区域在t0时刻人口y(t0)=y0为已知的,该区域的人口自然增长与人口总数成正比,所以t时刻的人口总数y(t)满足以下微分方程:,只要f(x,y)在a,bR1上连续,且关于y满足Lipschitz条件,即存在与x,y无关的常数L使(1.3)对任意定义在a,b上的y1(x)和y2(x)都成立,则上述IVP存在唯一解。,一阶常微分方程的初值问题/*Initial-ValueProblem*/:,(1.1),(1.2),定理1设f在区域(x,y)|axb,yR上连续,关于y满足利普希茨条件,则对任意x0a,b,y0R,常微分方程初值问题(1.1)式和(1.)式当xa,b时存在唯一的连续可微解y(x).,定理设f在区域上连续,且关于y满足利普希茨条件,设初值问题,的解为y(x,s),则,求解常微分方程有多种解析方法.但解析方法只能用来求解一些特殊类型的问题.,例:,要计算出解函数y(x)在一系列节点a=x0x1xn=b处的近似值,节点间距为步长,通常采用等距节点,即取hi=h(常数),这时节点为xn=x0+nh,n=0,1,2,。,实际问题中归结出来的微分方程要靠数值解法.,返回主页,9.2简单的数值方法,9.2.1欧拉法与后退欧拉法,在xy平面上,微分方程(1.1)的解y=y(x)称作它的积分曲线.积分曲线上一点(x,y)的切线斜率等于函数f(x,y)的值.如果按函数f(x,y)在xy平面上建立一个方向场,那么,积分曲线上每一点的切线方向均与方向场在该点的方向相一致.,基于上述几何解释,我们从初始点P0(x0,y0)出发,先依方向场在该点的方向推进到x=x1上一点P1,然后再从P1依方向场的方向推进到x=x2上一点P2,循此前进做出一条折线P0P1P2(图9-1).,图9-1,一般地,设已作出该折线的顶点Pn,过Pn(xn,yn)依方向场的方向再推进到Pn+1(xn+1,yn+1),显然两个顶点Pn,Pn+1的坐标有关系,即yn+1=yn+hf(xn,yn).(2.1),这就是著名的欧拉(Euler)公式.,若初值y0已知,则依公式(2.1)可逐步算出,y1=y0+hf(x0,y0),y2=y1+hf(x1,y1),例1求解初值问题,(2.2),解为便于进行比较,本章将用多种数值方法求解上述初值问题.这里先用欧拉方法,欧拉公式的具体形式为,取步长y=0.1,计算结果见表9-1.,表9-1计算结果对比,初值问题(2.2)有解y=(1+2x)1/2,按这个解析式子算出的准确值y(xn)同近似值yn一起列在表9-1中,两者相比较可以看出欧拉方法的精度很差.,表9-1计算结果对比,functionx,y=eulerivp(fun,x0,xf,y0,h)n=fix(xf-x0)/h)+1;y(1)=y0;x(1)=x0;fori=1:(n-1)x(i+1)=x(1)+i*h;y(i+1)=y(i)+h*feval(fun,x(i),y(i);end,functionf=fun(x,y)f=y-2*x/y,x,y=eulerivp(fun,0,1,1,0.1);plot(x,y,-r)y1=sqrt(1+2*x)plot(x,y,-r)holdonplot(x,y1,-b),还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上,那么,按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线(图9-2).从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线,可见欧拉方法是相当粗糙的.,图9-2,为了分析计算公式的精度,通常可用泰勒展开将y(xn+1)在xn处展开,则有,在yn=y(xn)的前提下,f(xn,yn)=f(xn,y(xn)=y(xn).于是可得欧拉法(2.1)的公式误差,(2.3),称为此方法的局部截断误差.,如果对方程(1.1)从xn到xn+1积分,得,(2.4),如果在(2.4)中右端积分用右矩形公式hf(xn+1,y(xn+1)近似,则得另一个公式,yn+1=yn+hf(xn+1,yn+1),(2.5),称为后退的欧拉法.,右端积分用左矩形公式hf(xn,y(xn)近似,再以yn代替y(xn),yn+1代替y(xn+1)也得到(2.1),局部截断误差也是(2.3).,1.欧拉公式是关于yn+1的一个直接的计算公式,这类公式称作是显式的;2.后退的欧拉公式(2.5)的右端含有未知的yn+1,它实际上是关于yn+1的一个函数方程,这类公式称作是隐式的.,隐式方程(2.5)通常用迭代法求解,而迭代过程的实质是逐步显示化.,后退的欧拉公式与欧拉公式有着本质的区别:,设用欧拉公式,给出迭代初值,直接计算得,然后再用,代入(2.5)式,又有,如此反复进行,得,用它代入(2.5)式的右端,使之转化为显式,(2.6),由于f(x,y)对y满足利普希茨条件(1.3).由(2.6)减(2.5)得,由此可知,只要hL,我们反复将步长折半进行计算,直至为止,这时取最终得到的yn+1(h/2)作为结果;,2.如果为止,这时再将步长折半一次,就得到所要的结果.,这种通过加倍或折半处理步长的方法称为变步长方法.表面上看,为了选择步长,每一步的计算量增加了,但总体考虑往往是合算的.,返回主页,9.4单步法的收敛性与稳定性/*ConvergencyandStability*/,收敛性/*Convergency*/,例:就初值问题考察欧拉显式格式的收敛性。,解:该问题的精确解为,欧拉公式为,对任意固定的x=xi=ih,有,稳定性/*Stability*/,例:考察初值问题在区间0,0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.00002.00004.00008.00001.60001013.2000101,1.00002.50001016.25001021.56251023.90631039.7656104,1.00002.50006.25001.56261013.90631019.7656101,1.00004.97871022.47881031.23411046.14421063.0590107,Whatiswrong?!,一般分析时为简单起见,只考虑试验方程/*testequation*/,常数,可以是复数,例:考察隐式欧拉法,可见绝对稳定区域为:,注:一般来说,隐式欧拉法的绝对稳定性比同阶的显式法的好。,例:隐式龙格-库塔法,而显式14阶方法的绝对稳定区域为,其中2阶方法的绝对稳定区域为,无条件稳定,返回主页,9.5线性多步法/*MultistepMethod*/,用若干节点处的y及y值的线性组合来近似y(xi+1)。,其通式可写为:,当10时,为隐式公式;1=0则为显式公式。,基于数值积分的构造法,亚当姆斯显式公式/*Adamsexplicitformulae*/,Newton插值余项,/*显式计算公式*/,局部截断误差为:,例:k=1时有,注:一般有,其中Bk与yi+1计算公式中fi,fik各项的系数均可查表得到。,亚当姆斯隐式公式/*Adamsimplicitformulae*/,小于Bk,较同阶显式稳定,亚当姆斯预测-校正系统/*Adamspredictor-correctorsystem*/,Step1:用Runge-Kutta法计算前k个初值;,Step2:用Adams显式计算预测值;,Step3:用同阶Adams隐式计算校正值。,注意:三步所用公式的精度必须相同。通常用经典Runge-Kutta法配合4阶Adams公式。,4阶Adams隐式公式的截断误差为,Predictedvaluepi+1,Modifiedvaluemi+1,Correctedvalueci+1,Modifiedfinalvalueyi+1,外推技术/*extrapolation*/,基于泰勒展开的构造法,将通式中的右端各项yi1,yik;fi+1,fi1,fik分别在xi点作泰勒展开,与精确解y(xi+1)在xi点的泰勒展开作比较。通过令同类项系数相等,得到足以确定待定系数0,k;1,0,k的等式,则可构造出线性多步法的公式。,解:,/*y(xi)=yi*/,个未知数个方程,7,5,令1=2=0,以yi+1取代yi1,并取1=2=0,取1=1,2=0得到辛甫生/*Simpson*/公式与Milne公式匹配使用,辛普森/*Simpson*/公式,在区间xi1,xi+1上积分,并用Simpson数值积分公式来近似积分项,亦可得此Simpson公式。,Milne-Simpson系统的缺点是稳定性差,为改善稳定性,考虑另一种隐式校正公式:,要求公式具有4阶精度。通过泰勒展开,可得到个等式,从中解出个未知数,则有个自由度。,5,6,1,取1=1得Simpson公式,注:哈明公式不能用数值积分方法推导出来。,返回主页,9.5Adams方法/*AdamsMethod*/,前述Runge-Kutta方法是一类重要方法,但这类方法的每一步需要先预报几个点上的斜率值,计算量比较大.,考察到在计算yn+1之前已得出一系列节点上的斜率值,自然会问,能否利用这些”老信息”来减少计算量呢?这就是Adams方法的设计思想.,用若干节点处的y及y值的线性组合来近似y(xi+1)。,其通式可写为:,当10时,为隐式公式;1=0则为显式公式。,特别地,Euler格式,和隐式Euler格式,是一阶Adams方法.,9.5.1二阶Adams格式,设用xn,xn-1两点的斜率值加权平均生成区间xn,xn+1上的平均斜率,而设计如下形式的差分格式,适当选取参数,使上述格式具有二阶精度.为此,考察其对应的近似关系式,设xn=0,h=1,这里有,令对于,准确可定出,这样设计出的计算格式,二阶显式Adams格式,类似地,改用xn,xn+1两点的斜率值加权平均生成区间xn,xn+1上的平均斜率,而使格式,具有二阶精度,不难定出,从而有二阶隐式Adams格式,它是熟知的梯形格式.,9.5.2误差的事后估计,仿照改进的Euler格式的构造方式,可以将显式与隐式两种Adams格式匹配在一起,构成下列二阶Adams预报校正系统:,预报,校正,要用到更前一步信息,因此不能自行启动.,在实际计算时,可以先借助于某种单步法,譬如具有二阶精度的改进的Euler格式提供开始值y1,再启动上述预报校正系统计算下去.,上述预报校正技术不仅能设计出实用算法,而且还能用于误差的事后估计.为此再考察系统(20)中预报与校正两种格式,注意到它们均具有二阶精度,进一步将它们加工成具有三阶精度的计算格式,为此考察其对应的近似关系式:,不妨设,令对于y=x3准确成立,可定出,从而有,由于这里yn+1是具有三阶精度的”准确”值,因而可以用预报值与校正值两者的偏差来估计它们的误差:,(21),利用误差作为计算结果的一种补偿有可能改善精度,因而基于这种误差的事后估计可以进一步优化预报校正系统(20).,就是说,按(21),与,分别可以看作,与,的改进值.在校正值,求出之前,自然用上一步的偏差值,替代,进行计算,这样,系统(20)可修改为如下改进的二阶Adams预报校正系统:,预报,校正,改进,改进,因此在启动之前必须先提供开始值y1与p1-c1.同Adams预报校正系统(20)一样,开始值y1可用改进的Euler格式(9)提供,而p1-c1一般令其等于0.,9.5.3实用的四阶Adams预报校正系统,运用上述处理方法,不难导出如下显式与隐式四阶Adams格式,将两者匹配在一起,即可生成下列四阶Adams预报校正系统:,预报,校正,亚当姆斯预测-校正系统/*Adamspredictor-correctorsystem*/,Step1:用Runge-Kutta法计算前3个初值;,Step2:用Adams显式计算预测值;,Step3:用同阶Adams隐式计算校正值。,注意:三步所用公式的精度必须相同。通常用经典Runge-Kutta法配合4阶Adams公式。,例5用四阶Adams预报校正系统(24)求解初值问题(2.2),解:取步长h=0.1,用四阶Runge-Kutta格式(19)提供开始值,然后用四阶Adams系统(24)逐步计算.计算结果见表3.5.,functionA=Adams4PC(f,a,b,n,ya)ifn4break;endh=(b-a)/n;x=zeros(1,n+1);y=zeros(1,n+1);x=a:h:b;y(1)=ya;F=zeros(1,4);fori=1:nifi4k1=feval(f,x(i),y(i);k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1);k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2);k4=feval(f,x(i)+h,y(i)+h*k3);y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);elseF=feval(f,x(i-3:i),y(i-3:i);py=y(i)+(h/24)*(F*-9,37,-59,55);p=feval(f,x(i+1),py);F=F(2)F(3)F(4)p;y(i+1)=y(i)+(h/24)*(F*1,-5,19,9);endendA=x,y,A4=Adams4PC(fun,0,1,10,1),仿照二阶Adams格式的处理方法估计系统(24)中预报值pn+1与校正值cn+1的误差.为此,考察如下形式的5阶格式,不难定出,从而有误差估计式,(25),利用这一误差估计式修改四阶Adams预报校正系统(),即可导出下列改进的四阶Adams预报校正系统,预报,校正,改进,改进,9.6微分方程组与高阶方程/*SystemsofDifferentialEquationsandHigher-OrderEquations*/,一阶微分方程组,IVP的一般形式为:,前述所有公式皆适用于向量形式。,譬如,对于方程组,引进节

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论