《最新常微分方程数值解法51252PPT课件.ppt》由会员分享，可在线阅读，更多相关《最新常微分方程数值解法51252PPT课件.ppt（65页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、常微分方程数值解法51252对象象一一阶常微分方程初常微分方程初值问题:一一阶常微分方程常微分方程组初初值问题:高高阶常微分方程初常微分方程初值问题:(4.1)整体整体误差差ek=y(xk)-yk，下面，下面对其加以分析其加以分析y1=y0+hf(x0,y0)=1+0.1(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1(1.1-20.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438其精确解其精确解为欧拉法（欧拉法（续）思想思想:用中心差商近似代替微商用中心差商近似代替微商.注注：计计算算算算时时，先用，先用，先用，先用欧拉法欧拉法求出求出y1，
2、以后再用，以后再用，以后再用，以后再用二步欧拉法二步欧拉法计计算。算。算。算。二步欧拉法二步欧拉法(4.4)欧拉法（欧拉法（续）公式公式单步否步否 显式否式否单步步显式式单步步隐式式二步二步显式式截断截断误差差y(xn+1)-yn+1 截断截断误差差Def4.1 设y(xn)是是(4.1)式的精确解，式的精确解，yn是是(4.2)式欧式欧拉法得到的近似解，称拉法得到的近似解，称y(xn)-yn为欧拉法的整体截断欧拉法的整体截断误差差.Def4.3 若某算法的局部截断若某算法的局部截断误差差为O(hp+1),称称该算法算法有有p阶精度精度.Def4.2 假假设yn=y(xn),即第即第n步步计算
3、是精确的前提下算是精确的前提下,称称Rn+1=y(xn+1)-yn+1为欧拉法的局部截断欧拉法的局部截断误差差.分析：分析：证明其局部截断明其局部截断误差差为O(h2),可通可通过Taylor展开展开式式分析。分析。证明明：Euler 公式公式为令令yn=y(xn)，下下证：y(xn+1)-yn+1=O(h2)由由 y(x)=f(x,y(x)定理定理4.4 欧拉法的精度是一欧拉法的精度是一阶。二步欧拉法的二步欧拉法的局部截断局部截断误差差Def4.5 假假设yn=y(xn),yn1=y(xn1),称Rn+1=y(xn+1)-yn+1为二步二步欧拉法的欧拉法的局部截断局部截断误差差.定理定理4.
4、6 隐式欧拉法的精度是一式欧拉法的精度是一阶,二步欧拉法的精度是二二步欧拉法的精度是二阶。证明明：对二步欧拉法二步欧拉法进行行证明，明，考考虑其局部截断其局部截断误差，差，令令yn=y(xn),yn1=y(xn1),将上两式左右两端同将上两式左右两端同时相减：相减：二步欧拉法的二步欧拉法的局部截断局部截断误差差为O(h3)，其其精度是二精度是二阶。数数值积分法分法对右端的定右端的定积分用数分用数值积分公式求近似分公式求近似值：（1）用左矩形数用左矩形数值积分公式：分公式：（2）用梯形公式：用梯形公式：梯形公式梯形公式 梯形公式梯形公式:将将显示欧拉公式示欧拉公式,隐式欧拉公式平均可得式欧拉公式
5、平均可得 梯形公式是梯形公式是隐式式、单步步公式公式,其精度其精度为二二阶证：令：令yn=y(xn),由由Talor公式有公式有分析分析：梯形公式是：梯形公式是隐式公式，式公式，证明其明其局部局部截断截断误差差为O(h3)要用到要用到 二元函数的二元函数的Taylor公式。公式。f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1)=f(xn+1,y(xn+1)+fy(xn+1,)(yn+1-y(xn+1)，(xn，xn+1)=y(xn+1)+fy(xn+1,)(yn+1-y(xn+1)=y(xn)+hy”(xn)+O(h2)+fy(xn+1,)(yn+1-y(xn
6、+1)=f(xn,yn)+hy”(xn)+fy(xn+1,)(yn+1-y(xn+1)+O(h2)又又y(xn+1)=y(xn+h)=y(xn)+hy(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+hf(xn,yn)+hy”(xn)/2+O(h3)定理定理4.7 梯形公式的精度是二梯形公式的精度是二阶。从而从而y(xn+1)=yn+1+h fy(xn+1,)y(xn+1)-yn+1/2+O(h3)y(xn+1)-yn+1=h fy(xn+1,)y(xn+1)-yn+1/2+O(h3)y(xn+1)-yn+1=
7、O(h3)/1-hfy(xn+1,)/2=O(h3)梯形公式的截断梯形公式的截断误差差为O(h3)，其精度是2阶。f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,)(yn+1-y(xn+1)+O(h2)y(xn+1)=yn+hf(xn,yn)/2+hf(xn,yn)+hy”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h f(xn+1,yn+1)-fy(xn+1,)(yn+1-y(xn+1)+O(h2)/2+O(h3)=yn+hf(xn,yn)+f(xn+1,yn+1)/2 +h fy(xn+1,)(y(xn+1)-yn+1)/2+O(h3)解解:取取h=
8、0.01 x0=0 y0=y(0)=1 则 y(0.01)y1 f(x,y)=y，由梯形公式由梯形公式,基于基于稳定性考定性考虑解析解解析解欧拉公式的比欧拉公式的比较欧拉法欧拉法最简单最简单,精度低精度低隐式欧拉法隐式欧拉法稳定性好稳定性好二步欧拉法二步欧拉法显式显式,但需要两步初值但需要两步初值,且第且第2个初值只能个初值只能由其他方法给出由其他方法给出,可能对后面的递推精度可能对后面的递推精度有影响有影响梯形公式法梯形公式法精度有所提高精度有所提高,但为隐式但为隐式,需要迭代求解需要迭代求解,计算量大计算量大4.2 改改进的的Euler法法 Euler公式公式 计算量小，精度低算量小，精度
9、低梯形公式梯形公式 计算量大，精度高算量大，精度高综合两个公式，提出合两个公式，提出预报校正公式校正公式：预报 校正校正改改进的的Euler法法写成嵌套公式：写成嵌套公式：平均化形式：平均化形式：例例4.4 用改用改进的的Euler法解初法解初值问题在区在区间0,0.4上，上，步步长h=0.1的解，并比的解，并比较与精确解与精确解(y=1/(1-x)的差异。的差异。解解：Euler法法的的具具体体形形式式为：yn+1=yn+hyn2，改改进的的Euler法法的具体形式的具体形式为：x0=0,h=0.1,则 x1=0.1,x2=0.2,x3=0.3,x4=0.4 计算算y1:yp=y0+0.1y
10、02=1+0.112=1.1 yc=1+0.1X1.12=1.121 y1=(1.1+1.121)/21.1118同同样可求可求y2,、y3、,y4注注：（1）令令 y(xn)=yn，可可 推推 导 改改 进 的的 Euler法法 的的 局局 部部 截截 断断 误 差差 为O(h3),具有二具有二阶精度。精度。（2）改改进的的Euler法也可写成如下平均化形式法也可写成如下平均化形式4.3龙格格库塔方法塔方法 如何构造高如何构造高阶的方法？的方法？(4.5)为了构造函数了构造函数使得使得(4.5)式成式成为高高阶方法，方法，Taylor对于一般的于一般的显式式单步法步法若若讨论其精度（其精度（
11、阶数），即考数），即考虑令令则有有考考虑到到上述方法求高上述方法求高阶微分，微分，实际上不上不实用用改改进的的Euler公式公式:Euler公式公式：yn+1=yn+hf(xn,yn)写成写成精度精度：一一阶 精度精度：二二阶 由由Lagrange中中值定理，定理，称称为y(x)在在xn,xn+1上的平均斜率上的平均斜率取取Euler公式公式 取取改改进Euler公式公式 Euler公式用一个点的公式用一个点的值k1作作为k*的近似的近似值，而改，而改进的的Euler公式用二个点的公式用二个点的值k1和和k2的平均的平均值作作为k*近似近似值，改，改进的的Euler法比法比Euler法精度高。
12、法精度高。(4.6)Runge-Kutta法法的的思思想想：在在xn,xn+1内内多多预报几几个个点点的的ki值并并用用其其加加权平平均均值作作为k*近近似似值。而而构构造造出出具具有有更更高高精精度度的的计算公式。算公式。多个斜率加多个斜率加权平均可提高精度平均可提高精度Runge-Kutta法一般形式法一般形式：以以k1与与k2的加的加权平均来近似平均来近似k*，设其其中中c1,c2,2,b21为待待定定参参数数。适适当当选取取参参数数，使使(*)式式的的精精度度为二二阶，即使其局部截断，即使其局部截断误差差为O(h3)令令y(xn)=yn，由泰勒公式：由泰勒公式：二二阶龙格格库塔方法塔方
13、法(*)由多元函数的泰勒公式和由多元函数的泰勒公式和(*)式式则有有 比比较(a)与与(b)要使要使 y(xn+1)-yn+1=O(h3)(a a)(b b)上述方程上述方程组有四个未知量，只有三个方程，有无有四个未知量，只有三个方程，有无穷多多组解。解。取任意一取任意一组解便得一种二解便得一种二阶龙库公式。公式。当当c1=c2=1/2,a2=b21=1时二二阶Runge-Kutta公式公式为yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1)此即改此即改进Euler法法 取取c2=0，c2=1，a2=1/2，b21=1/2yn+1=yn+k2 k1=h
14、f(xn,yn)k2=hf(xn+h/2,yn+k1/2)此此为中点法或中点法或变形的形的 Euler公式公式三三阶龙格格库塔法是用三个塔法是用三个值k1,k2,k3的加的加权平均来近似平均来近似k*，即有即有yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2)要使其具有三要使其具有三阶精度，必精度，必须使局部截断使局部截断误差差为O(h4)类似二似二阶龙格格库塔法的推塔法的推导，c1,c2,c3,a2,a3,b21,b31,b32应满足足c1+c2+c3=1a2=b21a3=b
15、31+b32c2a2+c3a3=1/2c2a22+c3a32=1/3c3a32=1/6由由该方方程程组任任意意解解可可得三得三阶龙格格-库塔公式塔公式例例：Kutta公式公式kn+1=yn+(k1+4k2+k3)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h,yn-k1+2k2)类似可推出四似可推出四阶龙格格-库塔公式，常用的有塔公式，常用的有 例：例：经典典Runge-Kutta法法 yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h/2,yn+k2/2)k4=h
16、f(xn+h,yn+k3)局部截断局部截断误差差 O(h5)还有有:Gill公式及公式及m(m4)阶龙格格库塔法。塔法。m4时：计算量太大，精确度不一定提高，有算量太大，精确度不一定提高，有时会降低。会降低。Gill公式公式节省存省存储单元元控制舍入控制舍入误差差 对于于经典的四典的四阶Runge-Kutta法法给出如下算法：出如下算法：算法算法4.2求解：求解：dy/dx=f(x,y)axb y(a)=y0 Step 1:输入入a,b,y0 及及NStep 2:(b-a)/N=h,a=x,y0=yStep 3:输出出(x,y)Step 4:For I=1 T0 Nhf(x,y)=k1hf(x
17、+h/2,y+k1/2)=k2hf(x+h/2,y+k2/2)=k3hf(x+h,yk3)=k4y+(k1+2k2+2k3+k4)/6=yx+h=x输出出(x,y)Step 5:停止停止例例4.3用四用四阶经典典RungeKutta方法解初方法解初值问题：(1)求求，(2)求求，自适自适应龙格格库塔法塔法用户提出问题I：问题：如何判断：如何判断|y(xn)-yn|精度精度值y(xn)未知。未知。：如何取如何取h=?解解：如用：如用p阶龙格格库塔法塔法计算，局部截断算，局部截断误差差为O(hp+1)y=f(x,y)y(x0)=y0 要求要求误差差10-8求数求数值解解 xn h/2 xn+1h如
18、如 xn-xn+1 令令 yn=y(xn)yn+1(h)则 y(xn+1)-yn+1(h)chp+1步步长折半折半xnxn+h/2xn+1分两步分两步计算算y(xn+1)的近似的近似值yn+1(h/2)。则 y(xn+1)-yn+1(h/2)2c(h/2)p+1 定理：对于问题I若用P阶龙格库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h),yn+1(h/2)则有误差公式 注：注：10 误差的事后估差的事后估计法法 20 停机准停机准则：(可保可保证|y(xn+1)-yn+1(h/2)|）将将步步长折折半半反反复复计算算，直直至至为止止，取取h为最后一次的步最后一次的步长，yn
19、+1为最后一次最后一次计算的算的结果。果。Else if(为止止，最最后后一次运算的前一次一次运算的前一次计算算结果即果即为所需。所需。4.5 收收敛与与稳定性定性 对于常微分方程初于常微分方程初值问题(4.1)求求y=y(x)-求求y(xn)-求求yn xn=x0+nh离散化离散化某种公式某种公式例：例：Euler法，改法，改进的的Euler法，法，龙格格库塔法都是塔法都是单步法。步法。数数值解法：解法：单步法：步法：计算算yn+1时只用到前一步的只用到前一步的结果果yn。显式式单步法：步法：yn+1=yn+h(xn,yn,h)(x,y,h)为增量函数，它依赖于f,仅是xn,yn,h的函数D
20、ef:若某数若某数值方法方法对于任意固定的于任意固定的xn=x0+nh,当当 h-0(n-)时，yn-y(xn)，则称称该法是法是收收敛的。的。即即 (xn=x0+nh为固定固定值)改改进Euler法的收法的收敛性性 判定改判定改进Euler法的收法的收敛性：性：yn+1=yn+hf(xn,yn)+f(xn+h,yn+hf(xn,yn)/2 则(x,y,h)=f(x,y)+f(x+h,y+hf(x,y)/2若：若：yo=y(x0)f(x,y)关于关于y满足李足李-条件条件:则：限定限定hh0，则即即(x,y,h)满足李普希足李普希兹条件条件,由由Th改改进的的Euler法收法收敛同同样可可验证
21、，若，若f(x,y)关于关于y满足李普希足李普希兹条件条件,且且 y0=y(x0)时，Euler法，法，标准四准四阶龙格格库法也收法也收敛。4.5.2 稳定性定性在在讨论收收敛性性时，一一般般认可可：数数值方方法法本本身身计算算过程是准确的。程是准确的。实际上，并非如此：上，并非如此：初始初始值y0有有误差差0y0-y(x0).计算的每一步有舍入算的每一步有舍入误差。差。这些初始些初始误差在差在计算算过程的程的传播中，是逐步衰减播中，是逐步衰减的，的，还是是恶性增性增长，这就是就是稳定性定性问题。Def4.1若若一一种种数数值方方法法在在节点点xn处的的数数值解解yn的的扰动n0，而在以后的各
22、而在以后的各节点点值ym(mn)上有上有扰动m.当当|m|n|，(m=n+1,n+2,)，则称称该数数值算算法法是是稳定定的。的。分析某算法的数分析某算法的数值稳定性，通常考察模型方程定性，通常考察模型方程 y=y ，(0)Def4.1：设在在节点点xn处用用数数值法法得得到到的的理理想想数数值解解为 yn，而而 实 际 计 算算 得得 到到 的的 近近 似似 值 为 ，称称 为第第n步数步数值解的解的扰动Euler法的法的稳定性定性Euler法：yn+1=yn+hf(xn,yn)考察模型方程考察模型方程 y=y，(0)即即yn+1=(1+h)yn假假设在在节点点值yn上上有有扰动n，在在yn
23、+1上上有有扰动n+1，且且n+1仅由由n引起引起(计算算过程不再引程不再引进新的新的误差差)Euler法法稳定定Euler法的法的稳定的条件定的条件为 0h-2/隐式式Euler法法稳定性定性隐式式Euler法法：yn+1=yn+hf(xn+1,yn+1)对于模型方程于模型方程 y=y(yn+1=yn/(1-h)yn+1的的扰动0 上式均成立上式均成立，隐式式Euler法无条件法无条件稳定定210ReImg梯形公式梯形公式稳定性定性梯形公式梯形公式 yn+1=yn+hf(xn,yn)+f(xn+1,yn+1)/2，设模型方程模型方程为y=y(0)则 yn+1=yn+hyn+yn+1/2 0时
24、恒成立恒成立 梯形公式恒梯形公式恒稳定。定。yn+1的的扰动 3 3阶阶RungeRungeKuttaKutta显式式 1 4 阶方法的方法的绝对稳定区域定区域为k=1k=2k=3k=4-1-2-3-123ReImgTaylor公式公式公式公式单步否步否显式否式否精度精度 单步步显式式一一阶单步步隐式式一一阶二步二步显式式二二阶数数值积分法分法单步步隐式式二二阶单步步显式式二二阶RungeKutta方法方法局部截断局部截断误差差整体截断整体截断误差差收收敛与与稳定性定性4.6 多步法多步法某一步解和前若干步解的某一步解和前若干步解的值有关有关m步步线性多步法的一般形式性多步法的一般形式其中其中
25、a0,am-1,b0,bm为和和h无关的常数无关的常数初初值为y0,y1,ym-1若若bm0，方法是方法是隐式的式的(implicit)否否则是是显式的式的(explicit)(4.7)AdamsBashforth四四阶方法方法AdamsMoulton四四阶方法方法多步法的截断多步法的截断误差差4.7 常微分方程常微分方程组和高和高阶微分方程的数微分方程的数值解法解法形如形如欧拉方法的欧拉方法的计算公式算公式隐式欧拉方法的式欧拉方法的计算公式算公式梯形公式梯形公式Adams公式公式R-K方法方法对于于高高阶微微分分方方程程，总可可以以化化成成方方程程组的的形形式式，例如例如可以化成可以化成例例
26、令令y(x)=z，则有有用用4阶RK方法，有方法，有4.8 常微分方程常微分方程边值问题的数的数值解法解法形如形如第一第一类边界条件界条件(4.8)第二第二类边界条件界条件第三第三类边界条件界条件定理定理 假定假定问题（4.8）式中的函数）式中的函数f在集合在集合上上连续，且偏，且偏导数数fy和和fy也在也在D上上连续。若有。若有则边值问题有唯一解。有唯一解。例例 边值问题有有因此因此边值问题有唯一解。有唯一解。定理定理 若若线性性边值问题满足足则边值问题有唯一解。有唯一解。打靶法打靶法（Shooting Method）其中，其中，y1(x)是初是初值问题的解的解的解，的解，y2(x)是初是初值问题有限差分法有限差分法(Finite-Difference Methods)以差商代替以差商代替导数，将微分方程离散成数，将微分方程离散成为差分方程差分方程用用Taylor展开方法展开方法两式相加得到两式相加得到由中由中值定理得定理得线性性边值问题可写成可写成然后可得到有限差分公式然后可得到有限差分公式上式可写成上式可写成写成矩写成矩阵形式形式其中其中